Re: How to create a 3D matrix?

["Dr. David Kirkby" <drkirkby at ntlworld dot com>]

"Robert G. Brown" wrote:

> Numerical Recipes has a nasty license requiring one to pay, pay and pay
> to use their code in any sort of real project, especially one in a
> cluster environment and (as you note) has code that ranges from decent
> to really poor.  There is a famous online rant against the code quality
> of NR that should show up in google.  One reason I was ecstatic when the
> GSL project came out was so that I could get away from using NR for
> anything, and have just about fully excised NR from all my projects.

The memory allocation routines were put into the public domain by the
NR authors, so that is not an issue. They are in a separate appendix,
with their own 'public domain license' - don't quote me literally on
the term, but it is less restrictive than the rest of the book. 
 
> The idea of using malloc's wrapped inside array creation routines to
> allocate space in rectangular matrices is sound enough and I'm sure was
> common practice long before NR was written.  It is very simple to write
> a routine that allocates space in a contiguous nrows X ncols block and
> repacks its addresses into a **array (or more, for higher order tensors)
> and not TOO difficult to write routines that allocate space efficiently
> for e.g. spherical harmonic (l,m) indexed arrays.  If anyone is
> interested in seeing an example I can post one.

I'd be interested, as it is something I have thought about. It would
seem more logical and efficient to me - or at least not less
efficient.

> Numerical efficiency is likely to be better for a contiguous array
> allocation, but a lot still depends on stride and access topology.  

I've been wondering what the effect is when the code is multi-threaded
on a symmetric multi processor (SMP) machine, with a different thread
accessing different parts of the array. That's not so clear in that
case. I've seen arguments to suggest that each thread should access a
different area of memory. 

> one is using lots of linear algebra, one should likely try to use the
> ATLAS optimized BLAS libraries, which automagically block out matrices
> and switch algorithms to get maximal performance from the particular
> cache and memory architecture one is running on.
 
That's something I should look at myself. I've certainly had some odd
results when accessing arrays using the code from the NR book in a way
that is suitable for use on a SMP machine. A 600 MHz Dec Alpha I have
seems to be particularly slow in this case, being only twice the speed
of a 75 MHz Sun and about 6x slower than a 400 MHz HP C3000
workstation. I think these results
- see http://www.medphys.ucl.ac.uk/~davek/
are due to problems with the accessing memory in an efficient way. 


-- 
Dr. David Kirkby,
Senior Research Fellow,
Department of Medical Physics,
University College London,
11-20 Capper St, London, WC1E 6JA.
Tel: 020 7679 6408 Fax: 020 7679 6269
Internal telephone: ext 46408
e-mail davek@medphys.ucl.ac.uk