c1-2 (779456), страница 3
Текст из файла (страница 3)
Let us, however, defer it fora moment in favor of an even more fundamental matter, that of variable dimensionarrays (FORTRAN terminology) or conformant arrays (Pascal terminology). Theseare arrays that need to be passed to a function along with real-time informationabout their two-dimensional size. The systems programmer rarely deals with twodimensional arrays, and almost never deals with two-dimensional arrays whose sizeis variable and known only at run time. Such arrays are, however, the bread andbutter of scientific computing.
Imagine trying to live with a matrix inversion routinethat could work with only one size of matrix!There is no technical reason that a C compiler could not allow a syntax like211.2 Some C Conventions for Scientific Computing**m[0][1][0][2][0][3][0][4][1][0][1][1][1][2][1][3][1][4][2][0][2][1][2][2][2][3][2][4]*m[0][0][0][0][1][0][2][0][3][0][4]*m[1][1][0][1][1][1][2][1][3][1][4]*m[2][2][0][2][1][2][2][2][3][2][4](a)**m(b)Figure 1.2.1.
Two storage schemes for a matrix m. Dotted lines denote address reference, while solidlines connect sequential memory locations. (a) Pointer to a fixed size two-dimensional array. (b) Pointerto an array of pointers to rows; this is the scheme adopted in this book.float a[13][9],**aa;int i;aa=(float **) malloc((unsigned) 13*sizeof(float*));for(i=0;i<=12;i++) aa[i]=a[i];a[i] is a pointer to a[i][0]The identifier aa is now a matrix with index range aa[0..12][0..8].
You can useor modify its elements ad lib, and more importantly you can pass it as an argumentto any function by its name aa. 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.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).[0][0]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.float **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);Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.
To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).A typical use is1.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 functionfrees 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:Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).void free_submatrix(float **b, long nrl, long nrh, long ncl, long nch)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 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.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
