c2-7 (779465), страница 3
Текст из файла (страница 3)
(If there areno off-diagonal elements for that row, it is one greater than the index in sa of themost recently stored element of a previous row.)• Location 1 of ija is always equal to N + 2. (It can be read to determine N .)• Location N + 1 of ija is one greater than the index in sa of the last off-diagonalelement of the last row. (It can be read to determine the number of nonzeroelements in the matrix, or the number of elements in the arrays sa and ija.)Location N + 1 of sa is not used and can be set arbitrarily.• Entries in sa at locations ≥ N + 2 contain A’s off-diagonal values, ordered byrows and, within each row, ordered by columns.• Entries in ija at locations ≥ N +2 contain the column number of the correspondingelement in sa.While these rules seem arbitrary at first sight, they result in a rather elegant storagescheme.
As an example, consider the matrix3. 0. 1. 0. 0. 0. 4. 0. 0. 0. (2.7.27) 0. 7. 5. 9. 0. 0. 0. 0. 0. 2. 0. 0. 0. 6. 5.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).det A = det P det(S − R · P−1 · Q) = det S det(P − Q · S−1 · R)792.7 Sparse Linear SystemsIn row-indexed compact storage, matrix (2.7.27) is represented by the two arrays of length11, as follows1234567891011ija[k]78810111232454sa[k]3.4.5.0.5.x1.7.9.2.6.(2.7.28)Here x is an arbitrary value. Notice that, according to the storage rules, the value of N(namely 5) is ija[1]-2, and the length of each array is ija[ija[1]-1]-1, namely 11.The diagonal element in row i is sa[i], and the off-diagonal elements in that row are insa[k] where k loops from ija[i] to ija[i+1]-1, if the upper limit is greater or equal tothe lower one (as in C’s for loops).Here is a routine, sprsin, that converts a matrix from full storage mode into row-indexedsparse storage mode, throwing away any elements that are less than a specified threshold.Of course, the principal use of sparse storage mode is for matrices whose full storage modewon’t fit into your machine at all; then you have to generate them directly into sparse format.Nevertheless sprsin is useful as a precise algorithmic definition of the storage scheme, forsubscale testing of large problems, and for the case where execution time, rather than storage,furnishes the impetus to sparse storage.#include <math.h>void sprsin(float **a, int n, float thresh, unsigned long nmax, float sa[],unsigned long ija[])Converts a square matrix a[1..n][1..n] into row-indexed sparse storage mode.
Only elements of a with magnitude ≥thresh are retained. Output is in two linear arrays with dimension nmax (an input parameter): sa[1..] contains array values, indexed by ija[1..]. Thenumber of elements filled of sa and ija on output are both ija[ija[1]-1]-1 (see text).{void nrerror(char error_text[]);int i,j;unsigned long k;for (j=1;j<=n;j++) sa[j]=a[j][j];Store diagonal elements.ija[1]=n+2;Index to 1st row off-diagonal element, if any.k=n+1;for (i=1;i<=n;i++) {Loop over rows.for (j=1;j<=n;j++) {Loop over columns.if (fabs(a[i][j]) >= thresh && i != j) {if (++k > nmax) nrerror("sprsin: nmax too small");sa[k]=a[i][j];Store off-diagonal elements and their columns.ija[k]=j;}}ija[i+1]=k+1;As each row is completed, store index to}next.}The single most important use of a matrix in row-indexed sparse storage mode is tomultiply a vector to its right. In fact, the storage mode is optimized for just this purpose.The following routine is thus very simple.void sprsax(float sa[], unsigned long ija[], float x[], float b[],unsigned long n)Multiply a matrix in row-index sparse storage arrays sa and ija by a vector x[1..n], givinga vector b[1..n].{void nrerror(char error_text[]);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).index k80Chapter 2.Solution of Linear Algebraic Equationsunsigned long i,k;if (ija[1] != n+2) nrerror("sprsax: mismatched vector and matrix");for (i=1;i<=n;i++) {b[i]=sa[i]*x[i];Start with diagonal term.for (k=ija[i];k<=ija[i+1]-1;k++)Loop over off-diagonal terms.b[i] += sa[k]*x[ija[k]];It is also simple to multiply the transpose of a matrix by a vector to its right.
(We will usethis operation later in this section.) Note that the transpose matrix is not actually constructed.void sprstx(float sa[], unsigned long ija[], float x[], float b[],unsigned long n)Multiply the transpose of a matrix in row-index sparse storage arrays sa and ija by a vectorx[1..n], giving a vector b[1..n].{void nrerror(char error_text[]);unsigned long i,j,k;if (ija[1] != n+2) nrerror("mismatched vector and matrix in sprstx");for (i=1;i<=n;i++) b[i]=sa[i]*x[i];Start with diagonal terms.for (i=1;i<=n;i++) {Loop over off-diagonal terms.for (k=ija[i];k<=ija[i+1]-1;k++) {j=ija[k];b[j] += sa[k]*x[i];}}}(Double precision versions of sprsax and sprstx, named dsprsax and dsprstx, are usedby the routine atimes later in this section.
You can easily make the conversion, or else getthe converted routines from the Numerical Recipes diskettes.)In fact, because the choice of row-indexed storage treats rows and columns quitedifferently, it is quite an involved operation to construct the transpose of a matrix, given thematrix itself in row-indexed sparse storage mode.
When the operation cannot be avoided,it is done as follows: An index of all off-diagonal elements by their columns is constructed(see §8.4). The elements are then written to the output array in column order. As eachelement is written, its row is determined and stored. Finally, the elements in each columnare sorted by row.void sprstp(float sa[], unsigned long ija[], float sb[], unsigned long ijb[])Construct the transpose of a sparse square matrix, from row-index sparse storage arrays sa andija into arrays sb and ijb.{void iindexx(unsigned long n, long arr[], unsigned long indx[]);Version of indexx with all float variables changed to long.unsigned long j,jl,jm,jp,ju,k,m,n2,noff,inc,iv;float v;n2=ija[1];Linear size of matrix plus 2.for (j=1;j<=n2-2;j++) sb[j]=sa[j];Diagonal elements.iindexx(ija[n2-1]-ija[1],(long *)&ija[n2-1],&ijb[n2-1]);Index all off-diagonal elements by their columns.jp=0;for (k=ija[1];k<=ija[n2-1]-1;k++) {Loop over output off-diagonal elements.m=ijb[k]+n2-1;Use index table to store by (former) columns.sb[k]=sa[m];for (j=jp+1;j<=ija[m];j++) ijb[j]=k;Fill in the index to any omitted rows.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 Systems81}for (j=jp+1;j<n2;j++) ijb[j]=ija[n2-1];for (j=1;j<=n2-2;j++) {Make a final pass to sort each row byjl=ijb[j+1]-ijb[j];Shell sort algorithm.noff=ijb[j]-1;inc=1;do {inc *= 3;inc++;} while (inc <= jl);do {inc /= 3;for (k=noff+inc+1;k<=noff+jl;k++) {iv=ijb[k];v=sb[k];m=k;while (ijb[m-inc] > iv) {ijb[m]=ijb[m-inc];sb[m]=sb[m-inc];m -= inc;if (m-noff <= inc) break;}ijb[m]=iv;sb[m]=v;}} while (inc > 1);}}The above routine embeds internally a sorting algorithm from §8.1, but calls the externalroutine iindexx to construct the initial column index.
This routine is identical to indexx,as listed in §8.4, except that the latter’s two float declarations should be changed to long.(The Numerical Recipes diskettes include both indexx and iindexx.) In fact, you canoften use indexx without making these changes, since many computers have the propertythat numerical values will sort correctly independently of whether they are interpreted asfloating or integer values.As final examples of the manipulation of sparse matrices, we give two routines for themultiplication of two sparse matrices.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
