Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 26
Текст из файла (страница 26)
The routine loses its memory of the spanned conjugategradient subspace between calls, however, so you should not force it to return more oftenthan about every N iterations.Finally, note that linbcg is furnished in double precision, since it will be usually beused when N is quite large.#include <stdio.h>#include <math.h>#include "nrutil.h"#define EPS 1.0e-14void linbcg(unsigned long n, double b[], double x[], int itol, double tol,int itmax, int *iter, double *err)Solves A · x = b for x[1..n], given b[1..n], by the iterative biconjugate gradient method.On input x[1..n] should be set to an initial guess of the solution (or all zeros); itol is 1,2,3,or 4, specifying which convergence test is applied (see text); itmax is the maximum numberof allowed iterations; and tol is the desired convergence tolerance. On output, x[1..n] isreset to the improved solution, iter is the number of iterations actually taken, and err is theestimated error.
The matrix A is referenced only through the user-supplied routines atimes,which computes the product of either A or its transpose on a vector; and asolve, which solvesTA · x = b or A · x = b for some preconditioner matrix A (possibly the trivial diagonal part of A).{void asolve(unsigned long n, double b[], double x[], int itrnsp);void atimes(unsigned long n, double x[], double r[], int itrnsp);double snrm(unsigned long n, double sx[], int itol);unsigned long j;double ak,akden,bk,bkden,bknum,bnrm,dxnrm,xnrm,zm1nrm,znrm;double *p,*pp,*r,*rr,*z,*zz;Double precision is a good idea in this routine.2.7 Sparse Linear Systems87p=dvector(1,n);pp=dvector(1,n);r=dvector(1,n);rr=dvector(1,n);z=dvector(1,n);zz=dvector(1,n);Calculate initial residual.*iter=0;atimes(n,x,r,0);Input to atimes is x[1..n], output is r[1..n];for (j=1;j<=n;j++) {the final 0 indicates that the matrix (not itsr[j]=b[j]-r[j];transpose) is to be used.rr[j]=r[j];}/*atimes(n,r,rr,0); */Uncomment this line to get the “minimum residif (itol == 1) {ual” variant of the algorithm.bnrm=snrm(n,b,itol);asolve(n,r,z,0);Input to asolve is r[1..n], output is z[1..n];}the final 0 indicates that the matrix A (notelse if (itol == 2) {its transpose) is to be used.asolve(n,b,z,0);bnrm=snrm(n,z,itol);asolve(n,r,z,0);}else if (itol == 3 || itol == 4) {asolve(n,b,z,0);bnrm=snrm(n,z,itol);asolve(n,r,z,0);znrm=snrm(n,z,itol);} else nrerror("illegal itol in linbcg");while (*iter <= itmax) {Main loop.++(*iter);Tasolve(n,rr,zz,1);Final 1 indicates use of transpose matrix A .for (bknum=0.0,j=1;j<=n;j++) bknum += z[j]*rr[j];Calculate coefficient bk and direction vectors p and pp.if (*iter == 1) {for (j=1;j<=n;j++) {p[j]=z[j];pp[j]=zz[j];}}else {bk=bknum/bkden;for (j=1;j<=n;j++) {p[j]=bk*p[j]+z[j];pp[j]=bk*pp[j]+zz[j];}}bkden=bknum;Calculate coefficient ak, new iterate x, and newatimes(n,p,z,0);residuals r and rr.for (akden=0.0,j=1;j<=n;j++) akden += z[j]*pp[j];ak=bknum/akden;atimes(n,pp,zz,1);for (j=1;j<=n;j++) {x[j] += ak*p[j];r[j] -= ak*z[j];rr[j] -= ak*zz[j];}asolve(n,r,z,0);Solve A · z = r and check stopping criterion.if (itol == 1)*err=snrm(n,r,itol)/bnrm;else if (itol == 2)*err=snrm(n,z,itol)/bnrm;88Chapter 2.Solution of Linear Algebraic Equationselse if (itol == 3 || itol == 4) {zm1nrm=znrm;znrm=snrm(n,z,itol);if (fabs(zm1nrm-znrm) > EPS*znrm) {dxnrm=fabs(ak)*snrm(n,p,itol);*err=znrm/fabs(zm1nrm-znrm)*dxnrm;} else {*err=znrm/bnrm;Error may not be accurate, so loop again.continue;}xnrm=snrm(n,x,itol);if (*err <= 0.5*xnrm) *err /= xnrm;else {*err=znrm/bnrm;Error may not be accurate, so loop again.continue;}}printf("iter=%4d err=%12.6f\n",*iter,*err);if (*err <= tol) break;}free_dvector(p,1,n);free_dvector(pp,1,n);free_dvector(r,1,n);free_dvector(rr,1,n);free_dvector(z,1,n);free_dvector(zz,1,n);}The routine linbcg uses this short utility for computing vector norms:#include <math.h>double snrm(unsigned long n, double sx[], int itol)Compute one of two norms for a vector sx[1..n], as signaled by itol.
Used by linbcg.{unsigned long i,isamax;double ans;if (itol <= 3) {ans = 0.0;for (i=1;i<=n;i++) ans += sx[i]*sx[i];Vector magnitude norm.return sqrt(ans);} else {isamax=1;for (i=1;i<=n;i++) {Largest component norm.if (fabs(sx[i]) > fabs(sx[isamax])) isamax=i;}return fabs(sx[isamax]);}}So that the specifications for the routines atimes and asolve are clear, we list heresimple versions that assume a matrix A stored somewhere in row-index sparse format.extern unsigned long ija[];extern double sa[];The matrix is stored somewhere.void atimes(unsigned long n, double x[], double r[], int itrnsp){void dsprsax(double sa[], unsigned long ija[], double x[], double b[],unsigned long n);2.7 Sparse Linear Systems89void dsprstx(double sa[], unsigned long ija[], double x[], double b[],unsigned long n);These are double versions of sprsax and sprstx.if (itrnsp) dsprstx(sa,ija,x,r,n);else dsprsax(sa,ija,x,r,n);}extern unsigned long ija[];extern double sa[];The matrix is stored somewhere.void asolve(unsigned long n, double b[], double x[], int itrnsp){unsigned long i;for(i=1;i<=n;i++) x[i]=(sa[i] != 0.0 ? b[i]/sa[i] : b[i]);The matrix A is the diagonal part of A, stored in the first n elements of sa.
Since thetranspose matrix has the same diagonal, the flag itrnsp is not used.}CITED REFERENCES AND FURTHER READING:Tewarson, R.P. 1973, Sparse Matrices (New York: Academic Press). [1]Jacobs, D.A.H. (ed.) 1977, The State of the Art in Numerical Analysis (London: AcademicPress), Chapter I.3 (by J.K. Reid). [2]George, A., and Liu, J.W.H. 1981, Computer Solution of Large Sparse Positive Definite Systems(Englewood Cliffs, NJ: Prentice-Hall). [3]NAG Fortran Library (Numerical Algorithms Group, 256 Banbury Road, Oxford OX27DE, U.K.).[4]IMSL Math/Library Users Manual (IMSL Inc., 2500 CityWest Boulevard, Houston TX 77042).
[5]Eisenstat, S.C., Gursky, M.C., Schultz, M.H., and Sherman, A.H. 1977, Yale Sparse Matrix Package, Technical Reports 112 and 114 (Yale University Department of Computer Science). [6]Knuth, D.E. 1968, Fundamental Algorithms, vol. 1 of The Art of Computer Programming (Reading,MA: Addison-Wesley), §2.2.6. [7]Kincaid, D.R., Respess, J.R., Young, D.M., and Grimes, R.G. 1982, ACM Transactions on Mathematical Software, vol. 8, pp.
302–322. [8]PCGPAK User’s Guide (New Haven: Scientific Computing Associates, Inc.). [9]Bentley, J. 1986, Programming Pearls (Reading, MA: Addison-Wesley), §9. [10]Golub, G.H., and Van Loan, C.F. 1989, Matrix Computations, 2nd ed. (Baltimore: Johns HopkinsUniversity Press), Chapters 4 and 10, particularly §§10.2–10.3. [11]Stoer, J., and Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: Springer-Verlag),Chapter 8.
[12]Baker, L. 1991, More C Tools for Scientists and Engineers (New York: McGraw-Hill). [13]Fletcher, R. 1976, in Numerical Analysis Dundee 1975, Lecture Notes in Mathematics, vol. 506,A. Dold and B Eckmann, eds. (Berlin: Springer-Verlag), pp. 73–89. [14]Saad, Y., and Schulz, M.
1986, SIAM Journal on Scientific and Statistical Computing, vol. 7,pp. 856–869. [15]Bunch, J.R., and Rose, D.J. (eds.) 1976, Sparse Matrix Computations (New York: AcademicPress).Duff, I.S., and Stewart, G.W. (eds.) 1979, Sparse Matrix Proceedings 1978 (Philadelphia:S.I.A.M.).90Chapter 2.Solution of Linear Algebraic Equations2.8 Vandermonde Matrices and ToeplitzMatricesIn §2.4 the case of a tridiagonal matrix was treated specially, because thatparticular type of linear system admits a solution in only of order N operations,rather than of order N 3 for the general linear problem.
When such particular typesexist, it is important to know about them. Your computational savings, should youever happen to be working on a problem that involves the right kind of particulartype, can be enormous.This section treats two special types of matrices that can be solved in of orderN 2 operations, not as good as tridiagonal, but a lot better than the general case.(Other than the operations count, these two types having nothing in common.)Matrices of the first type, termed Vandermonde matrices, occur in some problemshaving to do with the fitting of polynomials, the reconstruction of distributions fromtheir moments, and also other contexts.
In this book, for example, a Vandermondeproblem crops up in §3.5. Matrices of the second type, termed Toeplitz matrices,tend to occur in problems involving deconvolution and signal processing. In thisbook, a Toeplitz problem is encountered in §13.7.These are not the only special types of matrices worth knowing about. TheHilbert matrices, whose components are of the form aij = 1/(i + j − 1), i, j =1, .
. . , N can be inverted by an exact integer algorithm, and are very difficult toinvert in any other way, since they are notoriously ill-conditioned (see [1] for details).The Sherman-Morrison and Woodbury formulas, discussed in §2.7, can sometimesbe used to convert new special forms into old ones. Reference [2] gives some otherspecial forms. We have not found these additional forms to arise as frequently asthe two that we now discuss.Vandermonde MatricesA Vandermonde matrix of size N × N is completely determined by N arbitrarynumbers x1 , x2 , . . . , xN , in terms of which its N 2 components are the integer powersxj−1, i, j = 1, .
. . , N . Evidently there are two possible such forms, depending on whetheriwe view the i’s as rows, j’s as columns, or vice versa. In the former case, we get a linearsystem of equations that looks like this,1x1x21···x22...x2N···1x2...xN1. ..···−1xN1 c1y1 −1 xN2 · c 2 = y2 . . ..
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
