1 /***************************************************************************
5 ----------------------------------------------------------------------------
7 FUNCTION: general real matrix routines; includes
8 gaussj() for matrix inversion using
9 Gauss-Jordan method with full pivoting.
11 The routines in this module have come more or less from ref [1].
12 Note that, probably due to the heritage of ref [1] (which has a
13 FORTRAN version that was probably written first), the use of 1 as
14 the first element of an array (or vector) is used. This is accomplished
15 in memory by allocating, but not using, the 0 elements in each dimension.
16 While this wastes some memory, it allows the routines to be ported more
17 easily from FORTRAN (I suspect) as well as adhering to conventional
18 matrix notation. As a result, however, traditional ANSI C convention
19 (0-base indexing) is not followed; as the authors of ref [1] point out,
20 there is some question of the portability of the resulting routines
21 which sometimes access negative indexes. See ref [1] for more details.
23 ----------------------------------------------------------------------------
25 MODULE STATUS: developmental
27 ----------------------------------------------------------------------------
29 GENEALOGY: Created 950222 E. B. Jackson
31 ----------------------------------------------------------------------------
33 DESIGNED BY: from Numerical Recipes in C, by Press, et. al.
35 CODED BY: Bruce Jackson
39 ----------------------------------------------------------------------------
49 Revision 1.1 1999/04/05 21:32:45 curt
52 Revision 1.1 1998/06/27 22:34:57 curt
55 * Revision 1.1 1995/02/27 19:55:44 bjax
59 ----------------------------------------------------------------------------
61 REFERENCES: [1] Press, William H., et. al, Numerical Recipes in
62 C, 2nd edition, Cambridge University Press, 1992
64 ----------------------------------------------------------------------------
68 ----------------------------------------------------------------------------
72 ----------------------------------------------------------------------------
76 ----------------------------------------------------------------------------
80 --------------------------------------------------------------------------*/
94 #include "ls_matrix.h"
97 #define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}
99 static char rcsid[] = "$Id$";
102 int *nr_ivector(long nl, long nh)
106 v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
112 double **nr_matrix(long nrl, long nrh, long ncl, long nch)
113 /* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
115 long i, nrow=nrh-nrl+1, ncol=nch-ncl+1;
118 /* allocate pointers to rows */
119 m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
123 fprintf(stderr, "Memory failure in routine 'nr_matrix'.\n");
130 /* allocate rows and set pointers to them */
131 m[nrl] = (double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
134 fprintf(stderr, "Memory failure in routine 'matrix'\n");
141 for (i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
143 /* return pointer to array of pointers to rows */
148 void nr_free_ivector(int *v, long nl /* , long nh */)
150 free( (char *) (v+nl-NR_END));
154 void nr_free_matrix(double **m, long nrl, long nrh, long ncl, long nch)
155 /* free a double matrix allocated by nr_matrix() */
157 free((char *) (m[nrl]+ncl-NR_END));
158 free((char *) (m+nrl-NR_END));
162 int nr_gaussj(double **a, int n, double **b, int m)
164 /* Linear equation solution by Gauss-Jordan elimination. a[1..n][1..n] is */
165 /* the input matrix. b[1..n][1..m] is input containing the m right-hand */
166 /* side vectors. On output, a is replaced by its matrix invers, and b is */
167 /* replaced by the corresponding set of solution vectors. */
169 /* Note: this routine modified by EBJ to make b optional, if m == 0 */
172 int *indxc, *indxr, *ipiv;
173 int i, icol, irow, j, k, l, ll;
174 double big, dum, pivinv, temp;
176 int bexists = ((m != 0) || (b == 0));
178 indxc = nr_ivector(1,n); /* The integer arrays ipiv, indxr, and */
179 indxr = nr_ivector(1,n); /* indxc are used for pivot bookkeeping */
180 ipiv = nr_ivector(1,n);
182 for (j=1;j<=n;j++) ipiv[j] = 0;
184 for (i=1;i<=n;i++) /* This is the main loop over columns */
187 for (j=1;j<=n;j++) /* This is outer loop of pivot search */
193 if (fabs(a[j][k]) >= big)
201 if (ipiv[k] > 1) return -1;
205 /* We now have the pivot element, so we interchange rows, if needed, */
206 /* to put the pivot element on the diagonal. The columns are not */
207 /* physically interchanged, only relabeled: indxc[i], the column of the */
208 /* ith pivot element, is the ith column that is reduced, while indxr[i] */
209 /* is the row in which that pivot element was orignally located. If */
210 /* indxr[i] != indxc[i] there is an implied column interchange. With */
211 /* this form of bookkeeping, the solution b's will end up in the correct */
212 /* order, and the inverse matrix will be scrambed by columns. */
216 /* for (l=1;1<=n;l++) SWAP(a[irow][l],a[icol][l]) */
220 a[irow][l]=a[icol][l];
223 if (bexists) for (l=1;l<=m;l++) SWAP(b[irow][l],b[icol][l])
225 indxr[i] = irow; /* We are now ready to divide the pivot row */
226 indxc[i] = icol; /* by the pivot element, a[irow][icol] */
227 if (a[icol][icol] == 0.0) return -1;
228 pivinv = 1.0/a[icol][icol];
230 for (l=1;l<=n;l++) a[icol][l] *= pivinv;
231 if (bexists) for (l=1;l<=m;l++) b[icol][l] *= pivinv;
232 for (ll=1;ll<=n;ll++) /* Next, we reduce the rows... */
233 if (ll != icol ) /* .. except for the pivot one */
237 for (l=1;l<=n;l++) a[ll][l] -= a[icol][l]*dum;
238 if (bexists) for (l=1;l<=m;l++) b[ll][i] -= b[icol][l]*dum;
242 /* This is the end of the mail loop over columns of the reduction. It
243 only remains to unscrambled the solution in view of the column
244 interchanges. We do this by interchanging pairs of columns in
245 the reverse order that the permutation was built up. */
249 if (indxr[l] != indxc[l])
251 SWAP(a[k][indxr[l]],a[k][indxc[l]])
254 /* and we are done */
256 nr_free_ivector(ipiv,1 /*,n*/ );
257 nr_free_ivector(indxr,1 /*,n*/ );
258 nr_free_ivector(indxc,1 /*,n*/ );
260 return 0; /* indicate success */
263 void nr_copymat(double **orig, int n, double **copy)
264 /* overwrites matrix 'copy' with copy of matrix 'orig' */
268 if ((orig==0)||(copy==0)||(n==0)) return;
272 copy[i][j] = orig[i][j];
275 void nr_multmat(double **m1, int n, double **m2, double **prod)
279 if ((m1==0)||(m2==0)||(prod==0)||(n==0)) return;
285 for(k=1;k<=n;k++) prod[i][j] += m1[i][k]*m2[k][j];
291 void nr_printmat(double **a, int n)
299 printf("% 9.4f ", a[i][j]);
307 void testmat( void ) /* main() for test purposes */
309 double **mat1, **mat2, **mat3;
311 int loop, i, j, n = 20;
312 long maxlong = RAND_MAX;
315 invmaxlong = 1.0/(double)maxlong;
316 mat1 = nr_matrix(1, n, 1, n );
317 mat2 = nr_matrix(1, n, 1, n );
318 mat3 = nr_matrix(1, n, 1, n );
320 /* for(i=1;i<=n;i++) mat1[i][i]= 5.0; */
322 for(loop=0;loop<maxloop;loop++)
327 mat1[i][j] = 2.0 - 4.0*invmaxlong*(double) rand();
329 printf("Original matrix:\n");
330 nr_printmat( mat1, n );
332 nr_copymat( mat1, n, mat2 );
334 i = nr_gaussj( mat2, n, 0, 0 );
336 if (i) printf("Singular matrix.\n");
338 printf("Inverted matrix:\n");
339 nr_printmat( mat2, n );
341 nr_multmat( mat1, n, mat2, mat3 );
343 printf("Original multiplied by inverse:\n");
344 nr_printmat( mat3, n );
346 if (loop < maxloop-1) /* sleep(1) */;
349 nr_free_matrix( mat1, 1, n, 1, n );
350 nr_free_matrix( mat2, 1, n, 1, n );
351 nr_free_matrix( mat3, 1, n, 1, n );