c2-7 (779465), страница 2

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 2)

Then the matrix A is the tridiagonal part of thematrix in (2.7.9), with two terms modified:b01 = b1 − γ,b0N = bN − αβ/γ(2.7.11)#include "nrutil.h"void cyclic(float a[], float b[], float c[], float alpha, float beta,float r[], float x[], unsigned long n)Solves for a vector x[1..n] the “cyclic” set of linear equations given by equation (2.7.9). a,b, c, and r are input vectors, all dimensioned as [1..n], while alpha and beta are the cornerentries in the matrix. The input is not modified.{void tridag(float a[], float b[], float c[], float r[], float u[],unsigned long n);unsigned long i;float fact,gamma,*bb,*u,*z;if (n <= 2) nrerror("n too small in cyclic");bb=vector(1,n);u=vector(1,n);z=vector(1,n);gamma = -b[1];Avoid subtraction error in forming bb[1].bb[1]=b[1]-gamma;Set up the diagonal of the modified tridibb[n]=b[n]-alpha*beta/gamma;agonal system.for (i=2;i<n;i++) bb[i]=b[i];tridag(a,bb,c,r,x,n);Solve A · x = r.u[1]=gamma;Set up the vector u.u[n]=alpha;for (i=2;i<n;i++) u[i]=0.0;tridag(a,bb,c,u,z,n);Solve A · z = u.fact=(x[1]+beta*x[n]/gamma)/Form v · x/(1 + v · z).(1.0+z[1]+beta*z[n]/gamma);for (i=1;i<=n;i++) x[i] -= fact*z[i];Now get the solution vector x.free_vector(z,1,n);free_vector(u,1,n);free_vector(bb,1,n);}Woodbury FormulaIf you want to add more than a single correction term, then you cannot use (2.7.8)repeatedly, since without storing a new A−1 you will not be able to solve the auxiliaryproblems (2.7.7) efficiently after the first step.

Instead, you need the Woodbury formula,which is the block-matrix version of the Sherman-Morrison formula,(A + U · VT )−1ih= A−1 − A−1 · U · (1 + VT · A−1 · U)−1 · VT · A−1(2.7.12)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).We now solve equations (2.7.7) with the standard tridiagonal algorithm, and thenget the solution from equation (2.7.8).The routine cyclic below implements this algorithm.

We choose the arbitraryparameter γ = −b1 to avoid loss of precision by subtraction in the first of equations(2.7.11). In the unlikely event that this causes loss of precision in the second ofthese equations, you can make a different choice.76Chapter 2.Solution of Linear Algebraic EquationsHere A is, as usual, an N × N matrix, while U and V are N × P matrices with P < Nand usually P N . The inner piece of the correction term may become clearer if writtenas the tableau,U −1   · 1 + VT · A−1 · U ·  VT(2.7.13)where you can see that the matrix whose inverse is needed is only P × P rather than N × N .The relation between the Woodbury formula and successive applications of the ShermanMorrison formula is now clarified by noting that, if U is the matrix formed by columns out of theP vectors u1 , .

. . , uP , and V is the matrix formed by columns out of the P vectors v1 , . . . , vP ,        U ≡ u 1  · · · u P       V ≡ v1  · · · vP   then two ways of expressing the same correction to A are!PXA+uk ⊗ vk = (A + U · VT )(2.7.14)(2.7.15)k=1(Note that the subscripts on u and v do not denote components, but rather distinguish thedifferent column vectors.)Equation (2.7.15) reveals that, if you have A−1 in storage, then you can either make theP corrections in one fell swoop by using (2.7.12), inverting a P × P matrix, or else makethem by applying (2.7.5) P successive times.If you don’t have storage for A−1 , then you must use (2.7.12) in the following way:To solve the linear equation!PXA+(2.7.16)uk ⊗ vk · x = bk=1first solve the P auxiliary problemsA · z 1 = u1A · z 2 = u2···(2.7.17)A · z P = uPand construct the matrix Z by columns from the z’s obtained,      Z ≡ z1  · · · zP   (2.7.18)Next, do the P × P matrix inversionH ≡ (1 + VT · Z)−1(2.7.19)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).2.7 Sparse Linear Systems77Finally, solve the one further auxiliary problemA·y =b(2.7.20)In terms of these quantities, the solution is given byhix = y − Z · H · (VT · y)(2.7.21)Once in a while, you will encounter a matrix (not even necessarily sparse)that can be inverted efficiently by partitioning.

Suppose that the N × N matrixA is partitioned intoPA=RQS(2.7.22)where P and S are square matrices of size p × p and s × s respectively (p + s = N ).The matrices Q and R are not necessarily square, and have sizes p × s and s × p,respectively.If the inverse of A is partitioned in the same manner,"−1A=ePeReQeS#(2.7.23)e R,e Q,e ethen P,S, which have the same sizes as P, Q, R, S, respectively, can befound by either the formulase = (P − Q · S−1 · R)−1Pe = −(P − Q · S−1 · R)−1 · (Q · S−1 )Qe = −(S−1 · R) · (P − Q · S−1 · R)−1R(2.7.24)eS = S−1 + (S−1 · R) · (P − Q · S−1 · R)−1 · (Q · S−1 )or else by the equivalent formulase = P−1 + (P−1 · Q) · (S − R · P−1 · Q)−1 · (R · P−1 )Pe = −(P−1 · Q) · (S − R · P−1 · Q)−1Qe = −(S − R · P−1 · Q)−1 · (R · P−1 )R(2.7.25)eS = (S − R · P−1 · Q)−1The parentheses in equations (2.7.24) and (2.7.25) highlight repeated factors thatyou may wish to compute only once.

(Of course, by associativity, you can insteaddo the matrix multiplications in any order you like.) The choice between usinge or eequation (2.7.24) and (2.7.25) depends on whether you want PS to have theSample 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).Inversion by Partitioning78Chapter 2.Solution of Linear Algebraic Equationssimpler formula; or on whether the repeated expression (S − R · P−1 · Q)−1 is easierto calculate than the expression (P − Q · S−1 · R)−1 ; or on the relative sizes of Pand S; or on whether P−1 or S−1 is already known.Another sometimes useful formula is for the determinant of the partitionedmatrix,(2.7.26)Indexed Storage of Sparse MatricesWe have already seen (§2.4) that tri- or band-diagonal matrices can be stored in a compactformat that allocates storage only to elements which can be nonzero, plus perhaps a few wastedlocations to make the bookkeeping easier.

What about more general sparse matrices? When asparse matrix of dimension N × N contains only a few times N nonzero elements (a typicalcase), it is surely inefficient — and often physically impossible — to allocate storage for allN 2 elements. Even if one did allocate such storage, it would be inefficient or prohibitive inmachine time to loop over all of it in search of nonzero elements.Obviously some kind of indexed storage scheme is required, one that stores only nonzeromatrix elements, along with sufficient auxiliary information to determine where an elementlogically belongs and how the various elements can be looped over in common matrixoperations. Unfortunately, there is no one standard scheme in general use. Knuth [7] describesone method. The Yale Sparse Matrix Package [6] and ITPACK [8] describe several othermethods. For most applications, we favor the storage scheme used by PCGPACK [9], whichis almost the same as that described by Bentley [10], and also similar to one of the Yale SparseMatrix Package methods.

The advantage of this scheme, which can be called row-indexedsparse storage mode, is that it requires storage of only about two times the number of nonzeromatrix elements. (Other methods can require as much as three or five times.) For simplicity,we will treat only the case of square matrices, which occurs most frequently in practice.To represent a matrix A of dimension N × N , the row-indexed scheme sets up twoone-dimensional arrays, call them sa and ija. The first of these stores matrix element valuesin single or double precision as desired; the second stores integer values. The storage rules are:• The first N locations of sa store A’s diagonal matrix elements, in order. (Note thatdiagonal elements are stored even if they are zero; this is at most a slight storageinefficiency, since diagonal elements are nonzero in most realistic applications.)• Each of the first N locations of ija stores the index of the array sa that containsthe first off-diagonal element of the corresponding row of the matrix.

Характеристики

Список файлов книги