Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 25
Текст из файла (страница 25)
Here, sa and ija store the matrix A; sb and ijb store the matrixB. This routine computes all components of the matrix product (which may be non-sparse!),but stores only those whose magnitude exceeds thresh. On output, the arrays sc and ijc(whose maximum size is input as nmax) give the product matrix in row-index storage mode.For sparse matrix multiplication, this routine will often be preceded by a call to sprstp, so asto construct the transpose of a known matrix into sb, ijb.{void nrerror(char error_text[]);unsigned long i,ijma,ijmb,j,k,ma,mb,mbb;float sum;if (ija[1] != ijb[1]) nrerror("sprstm: sizes do not match");ijc[1]=k=ija[1];for (i=1;i<=ija[1]-2;i++) {Loop over rows of A,for (j=1;j<=ijb[1]-2;j++) {and rows of B.if (i == j) sum=sa[i]*sb[j]; else sum=0.0e0;mb=ijb[j];for (ma=ija[i];ma<=ija[i+1]-1;ma++) {Loop through elements in A’s row.
Convoluted logic, following, accounts for thevarious combinations of diagonal and off-diagonal elements.ijma=ija[ma];if (ijma == j) sum += sa[ma]*sb[j];else {while (mb < ijb[j+1]) {ijmb=ijb[mb];if (ijmb == i) {sum += sa[i]*sb[mb++];continue;} else if (ijmb < ijma) {mb++;continue;} else if (ijmb == ijma) {sum += sa[ma]*sb[mb++];continue;}break;}}}for (mbb=mb;mbb<=ijb[j+1]-1;mbb++) {Exhaust the remainder of B’s row.if (ijb[mbb] == i) sum += sa[i]*sb[mbb];}if (i == j) sc[i]=sum;Where to put the answer...else if (fabs(sum) > thresh) {if (k > nmax) nrerror("sprstm: nmax too small");sc[k]=sum;ijc[k++]=j;}}ijc[i+1]=k;}}Conjugate Gradient Method for a Sparse SystemSo-called conjugate gradient methods provide a quite general means for solving theN × N linear systemA·x =b(2.7.29)The attractiveness of these methods for large sparse systems is that they reference A onlythrough its multiplication of a vector, or the multiplication of its transpose and a vector.
As84Chapter 2.Solution of Linear Algebraic Equationswe have seen, these operations can be very efficient for a properly stored sparse matrix. You,the “owner” of the matrix A, can be asked to provide functions that perform these sparsematrix multiplications as efficiently as possible. We, the “grand strategists” supply the generalroutine, linbcg below, that solves the set of linear equations, (2.7.29), using your functions.The simplest, “ordinary” conjugate gradient algorithm [11-13] solves (2.7.29) only in thecase that A is symmetric and positive definite. It is based on the idea of minimizing the function1x·A·x−b·x2This function is minimized when its gradientf (x) =(2.7.30)∇f = A · x − b(2.7.31)is zero, which is equivalent to (2.7.29).
The minimization is carried out by generating asuccession of search directions pk and improved minimizers xk . At each stage a quantity αkis found that minimizes f (xk + αk pk ), and xk+1 is set equal to the new point xk + αk pk .The pk and xk are built up in such a way that xk+1 is also the minimizer of f over the wholevector space of directions already taken, {p1 , p2 , . .
. , pk }. After N iterations you arrive atthe minimizer over the entire vector space, i.e., the solution to (2.7.29).Later, in §10.6, we will generalize this “ordinary” conjugate gradient algorithm to theminimization of arbitrary nonlinear functions. Here, where our interest is in solving linear,but not necessarily positive definite or symmetric, equations, a different generalization isimportant, the biconjugate gradient method. This method does not, in general, have a simpleconnection with function minimization.
It constructs four sequences of vectors, rk , rk , pk ,pk , k = 1, 2, . . . . You supply the initial vectors r1 and r1 , and set p1 = r1 , p1 = r1 . Thenyou carry out the following recurrence:rk · rkpk · A · pk= rk − αk A · pkαk =rk+1rk+1 = rk − αk AT · pk(2.7.32)rk+1 · rk+1rk · rk= rk+1 + βk pkβk =pk+1pk+1 = rk+1 + βk pkThis sequence of vectors satisfies the biorthogonality conditionri · rj = ri · rj = 0,j <i(2.7.33)and the biconjugacy conditionpi · A · pj = pi · AT · pj = 0,j<i(2.7.34)There is also a mutual orthogonality,ri · pj = ri · pj = 0,j<i(2.7.35)The proof of these properties proceeds by straightforward induction [14]. 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)2.7 Sparse Linear Systems85while 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 functionΦ(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),(A−1−1· A) · x = A·b(2.7.41)The idea is that you might already be able to solve your linear system easily for some A close−1to A, in which case A · A ≈ 1, allowing the algorithm to converge in fewer steps. Thematrix A is called a preconditioner [11], and the overall scheme given here is known as thepreconditioned biconjugate gradient method or PBCG.For efficient implementation, the PBCG algorithm introduces an additional set of vectorszk and zk defined byA · zk = rkandTA · zk = rk(2.7.42)86Chapter 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)pk+1 = zk+1 + βk pkFor linbcg, below, we will ask you to supply routines that solve the auxiliary linear systems(2.7.42).
If you have no idea what to use for the preconditioner A, then use the diagonal partof 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−1|A−1· (A · x − b)|/|A· b| < tol(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.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
