c2-7 (779465), страница 5
Текст из файла (страница 5)
As long as therecurrence does not break down earlier because one of the denominators is zero, it mustterminate after m ≤ N steps with rm+1 = rm+1 = 0. This is basically because after at mostN steps you run out of new orthogonal directions to the vectors you’ve already constructed.To use the algorithm to solve the system (2.7.29), make an initial guess x1 for thesolution. Choose r1 to be the residualr1 = b − A · x1(2.7.36)and choose r1 = r1 .
Then form the sequence of improved estimatesxk+1 = xk + αk pk(2.7.37)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).1x·A·x−b·x2This function is minimized when its gradientf (x) =2.7 Sparse Linear Systems85Φ(x) =11r · r = |A · x − b|222(2.7.38)where the successiveiterates xk minimize Φ over the same set of search directions pk generatedin the conjugate gradient method. This algorithm has been generalized in various ways forunsymmetric matrices.
The generalized minimum residual method (GMRES; see [9,15] ) isprobably the most robust of these methods.Note that equation (2.7.38) gives∇Φ(x) = AT · (A · x − b)(2.7.39)For any nonsingular matrix A, AT · A is symmetric and positive definite.
You might thereforebe tempted to solve equation (2.7.29) by applying the ordinary conjugate gradient algorithmto the problem(AT · A) · x = AT · b(2.7.40)Don’t! The condition number of the matrix AT · A is the square of the condition number ofA (see §2.6 for definition of condition number). A large condition number both increases thenumber of iterations required, and limits the accuracy to which a solution can be obtained. Itis almost always better to apply the biconjugate gradient method to the original matrix A.So far we have said nothing about the rate of convergence of these methods.
Theordinary conjugate gradient method works well for matrices that are well-conditioned, i.e.,“close” to the identity matrix. This suggests applying these methods to the preconditionedform of equation (2.7.29),e(A−1e −1 · b· A) · x = A(2.7.41)e closeThe idea is that you might already be able to solve your linear system easily for some Ae −1 · A ≈ 1, allowing the algorithm to converge in fewer steps. Theto A, in which case Ae is called a preconditioner [11], and the overall scheme given here is known as thematrix Apreconditioned biconjugate gradient method or PBCG.For efficient implementation, the PBCG algorithm introduces an additional set of vectorszk and zk defined bye · zk = rkAande T · zk = rkA(2.7.42)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).while carrying out the recurrence (2.7.32). Equation (2.7.37) guarantees that rk+1 from therecurrence is in fact the residual b − A · xk+1 corresponding to xk+1 .
Since rm+1 = 0,xm+1 is the solution to equation (2.7.29).While there is no guarantee that this whole procedure will not break down or becomeunstable for general A, in practice this is rare. More importantly, the exact termination in atmost N iterations occurs only with exact arithmetic. Roundoff error means that you shouldregard the process as a genuinely iterative procedure, to be halted when some appropriateerror criterion is met.The ordinary conjugate gradient algorithm is the special case of the biconjugate gradientalgorithm when A is symmetric, and we choose r1 = r1 . Then rk = rk and pk = pk for allk; you can omit computing them and halve the work of the algorithm.
This conjugate gradientversion has the interpretation of minimizing equation (2.7.30). If A is positive definite aswell as symmetric, the algorithm cannot break down (in theory!). The routine linbcg belowindeed reduces to the ordinary conjugate gradient method if you input a symmetric A, butit does all the redundant computations.Another variant of the general algorithm corresponds to a symmetric but non-positivedefinite A, with the choice r1 = A · r1 instead of r1 = r1 . In this case rk = A · rk andpk = A · pk for all k. This algorithm is thus equivalent to the ordinary conjugate gradientalgorithm, but with all dot products a · b replaced by a · A · b.
It is called the minimum residualalgorithm, because it corresponds to successive minimizations of the function86Chapter 2.Solution of Linear Algebraic Equationsand modifies the definitions of αk , βk , pk , and pk in equation (2.7.32):rk · zkpk · A · pkrk+1 · zk+1βk =rk · zkpk+1 = zk+1 + βk pkαk =(2.7.43)For linbcg, below, we will ask you to supply routines that solve the auxiliary linear systemse then use the diagonal part(2.7.42). If you have no idea what to use for the preconditioner A,of A, or even the identity matrix, in which case the burden of convergence will be entirelyon the biconjugate gradient method itself.The routine linbcg, below, is based on a program originally written by Anne Greenbaum.(See [13] for a different, less sophisticated, implementation.) There are a few wrinkles youshould know about.What constitutes “good” convergence is rather application dependent.
The routinelinbcg therefore provides for four possibilities, selected by setting the flag itol on input.If itol=1, iteration stops when the quantity |A · x − b|/|b| is less than the input quantitytol. If itol=2, the required criterion is−1e|Ae −1 · b| < tol· (A · x − b)|/|A(2.7.44)If itol=3, the routine uses its own estimate of the error in x, and requires its magnitude,divided by the magnitude of x, to be less than tol.
The setting itol=4 is the same as itol=3,except that the largest (in absolute value) component of the error and largest component of xare used instead of the vector magnitude (that is, the L∞ norm instead of the L2 norm). Youmay need to experiment to find which of these convergence criteria is best for your problem.On output, err is the tolerance actually achieved.
If the returned count iter doesnot indicate that the maximum number of allowed iterations itmax was exceeded, then errshould be less than tol. If you want to do further iterations, leave all returned quantities asthey are and call the routine again. 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 solvese (possibly the trivial diagonal part of A).e · x = b or Ae T · x = b for some preconditioner matrix AA{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.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).pk+1 = zk+1 + βk pk2.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);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).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];e (not}the final 0 indicates that the matrix Aelse 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);eT .asolve(n,rr,zz,1);Final 1 indicates use of transpose matrix Afor (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];}e · z = r and check stopping criterion.asolve(n,r,z,0);Solve Aif (itol == 1)*err=snrm(n,r,itol)/bnrm;else if (itol == 2)*err=snrm(n,z,itol)/bnrm;88Chapter 2.Solution of Linear Algebraic Equationsfree_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.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
