Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C, страница 11

That function, which declares the correspondingdummy argument as float **aa;, can address its elements as aa[i][j] withoutknowing its physical size.You may rightly not wish to clutter your programs with code like the abovefragment. Also, there is still the outstanding problem of how to treat unit-offsetindices, so that (for example) the above matrix aa could be addressed with the rangea[1..13][1..9]. Both of these problems are solved by additional utility routinesin nrutil.c (Appendix B) which allocate and deallocate matrices of arbitraryrange. The synopses arefloat **matrix(long nrl, long nrh, long ncl, long nch)Allocates a float matrix with range [nrl..nrh][ncl..nch].double **dmatrix(long nrl, long nrh, long ncl, long nch)Allocates a double matrix with range [nrl..nrh][ncl..nch].int **imatrix(long nrl, long nrh, long ncl, long nch)Allocates an int matrix with range [nrl..nrh][ncl..nch].void free_matrix(float **m, long nrl, long nrh, long ncl, long nch)Frees a matrix allocated with matrix.22Chapter 1.Preliminariesvoid free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch)Frees a matrix allocated with dmatrix.void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch)Frees a matrix allocated with imatrix.A typical use isfloat **a;a=matrix(1,13,1,9);...a[3][5]=......+a[2][9]/3.0...someroutine(a,...);...free_matrix(a,1,13,1,9);All matrices in Numerical Recipes are handled with the above paradigm, and wecommend it to you.Some further utilities for handling matrices are also included in nrutil.c.The first is a function submatrix() that sets up a new pointer reference to analready-existing matrix (or sub-block thereof), along with new offsets if desired.Its synopsis isfloat **submatrix(float **a, long oldrl, long oldrh, long oldcl,long oldch, long newrl, long newcl)Point a submatrix [newrl..newrl+(oldrh-oldrl)][newcl..newcl+(oldch-oldcl)] tothe existing matrix range a[oldrl..oldrh][oldcl..oldch].Here oldrl and oldrh are respectively the lower and upper row indices of theoriginal matrix that are to be represented by the new matrix, oldcl and oldch arethe corresponding column indices, and newrl and newcl are the lower row andcolumn indices for the new matrix.

(We don’t need upper row and column indices,since they are implied by the quantities already given.)Two sample uses might be, first, to select as a 2 × 2 submatrix b[1..2][1..2] some interior range of an existing matrix, say a[4..5][2..3],float **a,**b;a=matrix(1,13,1,9);...b=submatrix(a,4,5,2,3,1,1);and second, to map an existing matrix a[1..13][1..9] into a new matrixb[0..12][0..8],float **a,**b;a=matrix(1,13,1,9);...b=submatrix(a,1,13,1,9,0,0);1.2 Some C Conventions for Scientific Computing23Incidentally, you can use submatrix() for matrices of any type whose sizeof()is the same as sizeof(float) (often true for int, e.g.); just cast the first argumentto type float ** and cast the result to the desired type, e.g., int **.The functionvoid free_submatrix(float **b, long nrl, long nrh, long ncl, long nch)frees the array of row-pointers allocated by submatrix(). Note that it does not freethe memory allocated to the data in the submatrix, since that space still lies withinthe memory allocation of some original matrix.Finally, if you have a standard C matrix declared as a[nrow][ncol], and youwant to convert it into a matrix declared in our pointer-to-row-of-pointers manner,the following function does the trick:float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch)Allocate a float matrix m[nrl..nrh][ncl..nch] that points to the matrix declared in thestandard C manner as a[nrow][ncol], where nrow=nrh-nrl+1 and ncol=nch-ncl+1.

Theroutine should be called with the address &a[0][0] as the first argument.(You can use this function when you want to make use of C’s initializer syntaxto set values for a matrix, but then be able to pass the matrix to programs in thisbook.) The functionvoid free_convert_matrix(float **b, long nrl, long nrh, long ncl, long nch)Free a matrix allocated by convert_matrix().frees the allocation, without affecting the original matrix a.The only examples of allocating a three-dimensional array as a pointer-topointer-to-pointer structure in this book are found in the routines rlft3 in §12.5 andsfroid in §17.4.

The necessary allocation and deallocation functions arefloat ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)Allocate a float 3-dimensional array with subscript range [nrl..nrh][ncl..nch][ndl..ndh].void free_f3tensor(float ***t, long nrl, long nrh, long ncl, long nch,long ndl, long ndh)Free a float 3-dimensional array allocated by f3tensor().Complex ArithmeticC does not have complex data types, or predefined arithmetic operations oncomplex numbers. That omission is easily remedied with the set of functions inthe file complex.c which is printed in full in Appendix C at the back of the book.A synopsis is as follows:24Chapter 1.Preliminariestypedef struct FCOMPLEX {float r,i;} fcomplex;fcomplex Cadd(fcomplex a, fcomplex b)Returns the complex sum of two complex numbers.fcomplex Csub(fcomplex a, fcomplex b)Returns the complex difference of two complex numbers.fcomplex Cmul(fcomplex a, fcomplex b)Returns the complex product of two complex numbers.fcomplex Cdiv(fcomplex a, fcomplex b)Returns the complex quotient of two complex numbers.fcomplex Csqrt(fcomplex z)Returns the complex square root of a complex number.fcomplex Conjg(fcomplex z)Returns the complex conjugate of a complex number.float Cabs(fcomplex z)Returns the absolute value (modulus) of a complex number.fcomplex Complex(float re, float im)Returns a complex number with specified real and imaginary parts.fcomplex RCmul(float x, fcomplex a)Returns the complex product of a real number and a complex number.The implementation of several of these complex operations in floating-pointarithmetic is not entirely trivial; see §5.4.Only about half a dozen routines in this book make explicit use of these complexarithmetic functions.

The resulting code is not as readable as one would like, becausethe familiar operations +-*/ are replaced by function calls. The C++ extension tothe C language allows operators to be redefined. That would allow more readablecode. However, in this book we are committed to standard C.We should mention that the above functions assume the ability to pass, return,and assign structures like FCOMPLEX (or types such as fcomplex that are definedto be structures) by value. All recent C compilers have this ability, but it is not inthe original K&R C definition. If you are missing it, you will have to rewrite thefunctions in complex.c, making them pass and return pointers to variables of typefcomplex instead of the variables themselves.

Likewise, you will need to modifythe recipes that use the functions.Several other routines (e.g., the Fourier transforms four1 and fourn) docomplex arithmetic “by hand,” that is, they carry around real and imaginary parts asfloat variables. This results in more efficient code than would be obtained by usingthe functions in complex.c.

But the code is even less readable. There is simply noideal solution to the complex arithmetic problem in C.Implicit Conversion of Float to DoubleIn traditional (K&R) C, float variables are automatically converted to doublebefore any operation is attempted, including both arithmetic operations and passingas arguments to functions. All arithmetic is then done in double precision.

If afloat variable receives the result of such an arithmetic operation, the high precision1.2 Some C Conventions for Scientific Computing25is immediately thrown away. A corollary of these rules is that all the real-numberstandard C library functions are of type double and compute to double precision.The justification for these conversion rules is, “well, there’s nothing wrong witha little extra precision,” and “this way the libraries need only one version of eachfunction.” One does not need much experience in scientific computing to recognizethat the implicit conversion rules are, in fact, sheer madness! In effect, they make itimpossible to write efficient numerical programs. One of the cultural barriers thatseparates computer scientists from “regular” scientists and engineers is a differingpoint of view on whether a 30% or 50% loss of speed is worth worrying about. Inmany real-time or state-of-the-art scientific applications, such a loss is catastrophic.The practical scientist is trying to solve tomorrow’s problem with yesterday’scomputer; the computer scientist, we think, often has it the other way around.The ANSI C standard happily does not allow implicit conversion for arithmeticoperations, but it does require it for function arguments, unless the function is fullyprototyped by an ANSI declaration as described earlier in this section.