c2-7 (779465), страница 4
Текст из файла (страница 4)
These are useful for techniques to be described in §13.10.In general, the product of two sparse matrices is not itself sparse. One therefore wantsto limit the size of the product matrix in one of two ways: either compute only those elementsof the product that are specified in advance by a known pattern of sparsity, or else compute allnonzero elements, but store only those whose magnitude exceeds some threshold value. Theformer technique, when it can be used, is quite efficient.
The pattern of sparsity is specifiedby furnishing an index array in row-index sparse storage format (e.g., ija). The programthen constructs a corresponding value array (e.g., sa). The latter technique runs the danger ofexcessive compute times and unknown output sizes, so it must be used cautiously.With row-index storage, it is much more natural to multiply a matrix (on the left) bythe transpose of a matrix (on the right), so that one is crunching rows on rows, rather thanrows on columns. Our routines therefore calculate A · BT , rather than A · B. This meansthat you have to run your right-hand matrix through the transpose routine sprstp beforesending it to the matrix multiply 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).jp=ija[m];Use bisection to find which row elementjl=1;m is in and put that into ijb[k].ju=n2-1;while (ju-jl > 1) {jm=(ju+jl)/2;if (ija[jm] > m) ju=jm; else jl=jm;}ijb[k]=jl;82Chapter 2.Solution of Linear Algebraic EquationsThe two implementing routines,sprspm for “pattern multiply” and sprstmfor “thresholdmultiply” are quite similar in structure.
Both are complicated by the logic of the variouscombinations of diagonal or off-diagonal elements for the two input streams and output stream.if (ija[1] != ijb[1] || ija[1] != ijc[1])nrerror("sprspm: sizes do not match");for (i=1;i<=ijc[1]-2;i++) {Loop over rows.j=m=i;Set up so that first pass through loop does themn=ijc[i];diagonal component.sum=sa[i]*sb[i];for (;;) {Main loop over each component to be output.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];}sc[m]=sum;sum=0.0;Reset indices for next pass through loop.if (mn >= ijc[i+1]) break;j=ijc[m=mn++];}}}#include <math.h>void sprstm(float sa[], unsigned long ija[], float sb[], unsigned long ijb[],float thresh, unsigned long nmax, float sc[], unsigned long ijc[])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 sprspm(float sa[], unsigned long ija[], float sb[], unsigned long ijb[],float sc[], unsigned long ijc[])Matrix multiply A · BT where A and B are two sparse matrices in row-index storage mode, andBT is the transpose of B.
Here, sa and ija store the matrix A; sb and ijb store the matrix B.This routine computes only those components of the matrix product that are pre-specified by theinput index array ijc, which is not modified. On output, the arrays sc and ijc give the productmatrix in row-index storage mode. For sparse matrix multiplication, this routine will often bepreceded by a call to sprstp, so as to construct the transpose of a known matrix into sb, ijb.{void nrerror(char error_text[]);unsigned long i,ijma,ijmb,j,m,ma,mb,mbb,mn;float sum;2.7 Sparse Linear Systems83if (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.
AsSample 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).Matrix multiply A · BT where A and B are two sparse matrices in row-index storage mode, andBT is the transpose of B. 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;84Chapter 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 function(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+1βk =rk · rkpk+1 = rk+1 + βk pkpk+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].
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
