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/06/17 18:07:34 curt
52 Revision 1.1.1.1 1999/04/05 21:32:45 curt
53 Start of 0.6.x branch.
55 Revision 1.1 1998/06/27 22:34:57 curt
58 * Revision 1.1 1995/02/27 19:55:44 bjax
62 ----------------------------------------------------------------------------
64 REFERENCES: [1] Press, William H., et. al, Numerical Recipes in
65 C, 2nd edition, Cambridge University Press, 1992
67 ----------------------------------------------------------------------------
71 ----------------------------------------------------------------------------
75 ----------------------------------------------------------------------------
79 ----------------------------------------------------------------------------
83 --------------------------------------------------------------------------*/
97 #include "ls_matrix.h"
100 #define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}
102 static char rcsid[] = "$Id$";
105 int *nr_ivector(long nl, long nh)
109 v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
115 double **nr_matrix(long nrl, long nrh, long ncl, long nch)
116 /* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
118 long i, nrow=nrh-nrl+1, ncol=nch-ncl+1;
121 /* allocate pointers to rows */
122 m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
126 fprintf(stderr, "Memory failure in routine 'nr_matrix'.\n");
133 /* allocate rows and set pointers to them */
134 m[nrl] = (double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
137 fprintf(stderr, "Memory failure in routine 'matrix'\n");
144 for (i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
146 /* return pointer to array of pointers to rows */
151 void nr_free_ivector(int *v, long nl /* , long nh */)
153 free( (char *) (v+nl-NR_END));
157 void nr_free_matrix(double **m, long nrl, long nrh, long ncl, long nch)
158 /* free a double matrix allocated by nr_matrix() */
160 free((char *) (m[nrl]+ncl-NR_END));
161 free((char *) (m+nrl-NR_END));
165 int nr_gaussj(double **a, int n, double **b, int m)
167 /* Linear equation solution by Gauss-Jordan elimination. a[1..n][1..n] is */
168 /* the input matrix. b[1..n][1..m] is input containing the m right-hand */
169 /* side vectors. On output, a is replaced by its matrix invers, and b is */
170 /* replaced by the corresponding set of solution vectors. */
172 /* Note: this routine modified by EBJ to make b optional, if m == 0 */
175 int *indxc, *indxr, *ipiv;
176 int i, icol, irow, j, k, l, ll;
177 double big, dum, pivinv, temp;
179 int bexists = ((m != 0) || (b == 0));
181 indxc = nr_ivector(1,n); /* The integer arrays ipiv, indxr, and */
182 indxr = nr_ivector(1,n); /* indxc are used for pivot bookkeeping */
183 ipiv = nr_ivector(1,n);
185 for (j=1;j<=n;j++) ipiv[j] = 0;
187 for (i=1;i<=n;i++) /* This is the main loop over columns */
190 for (j=1;j<=n;j++) /* This is outer loop of pivot search */
196 if (fabs(a[j][k]) >= big)
204 if (ipiv[k] > 1) return -1;
208 /* We now have the pivot element, so we interchange rows, if needed, */
209 /* to put the pivot element on the diagonal. The columns are not */
210 /* physically interchanged, only relabeled: indxc[i], the column of the */
211 /* ith pivot element, is the ith column that is reduced, while indxr[i] */
212 /* is the row in which that pivot element was orignally located. If */
213 /* indxr[i] != indxc[i] there is an implied column interchange. With */
214 /* this form of bookkeeping, the solution b's will end up in the correct */
215 /* order, and the inverse matrix will be scrambed by columns. */
219 /* for (l=1;1<=n;l++) SWAP(a[irow][l],a[icol][l]) */
223 a[irow][l]=a[icol][l];
226 if (bexists) for (l=1;l<=m;l++) SWAP(b[irow][l],b[icol][l])
228 indxr[i] = irow; /* We are now ready to divide the pivot row */
229 indxc[i] = icol; /* by the pivot element, a[irow][icol] */
230 if (a[icol][icol] == 0.0) return -1;
231 pivinv = 1.0/a[icol][icol];
233 for (l=1;l<=n;l++) a[icol][l] *= pivinv;
234 if (bexists) for (l=1;l<=m;l++) b[icol][l] *= pivinv;
235 for (ll=1;ll<=n;ll++) /* Next, we reduce the rows... */
236 if (ll != icol ) /* .. except for the pivot one */
240 for (l=1;l<=n;l++) a[ll][l] -= a[icol][l]*dum;
241 if (bexists) for (l=1;l<=m;l++) b[ll][i] -= b[icol][l]*dum;
245 /* This is the end of the mail loop over columns of the reduction. It
246 only remains to unscrambled the solution in view of the column
247 interchanges. We do this by interchanging pairs of columns in
248 the reverse order that the permutation was built up. */
252 if (indxr[l] != indxc[l])
254 SWAP(a[k][indxr[l]],a[k][indxc[l]])
257 /* and we are done */
259 nr_free_ivector(ipiv,1 /*,n*/ );
260 nr_free_ivector(indxr,1 /*,n*/ );
261 nr_free_ivector(indxc,1 /*,n*/ );
263 return 0; /* indicate success */
266 void nr_copymat(double **orig, int n, double **copy)
267 /* overwrites matrix 'copy' with copy of matrix 'orig' */
271 if ((orig==0)||(copy==0)||(n==0)) return;
275 copy[i][j] = orig[i][j];
278 void nr_multmat(double **m1, int n, double **m2, double **prod)
282 if ((m1==0)||(m2==0)||(prod==0)||(n==0)) return;
288 for(k=1;k<=n;k++) prod[i][j] += m1[i][k]*m2[k][j];
294 void nr_printmat(double **a, int n)
302 printf("% 9.4f ", a[i][j]);
310 void testmat( void ) /* main() for test purposes */
312 double **mat1, **mat2, **mat3;
314 int loop, i, j, n = 20;
315 long maxlong = RAND_MAX;
318 invmaxlong = 1.0/(double)maxlong;
319 mat1 = nr_matrix(1, n, 1, n );
320 mat2 = nr_matrix(1, n, 1, n );
321 mat3 = nr_matrix(1, n, 1, n );
323 /* for(i=1;i<=n;i++) mat1[i][i]= 5.0; */
325 for(loop=0;loop<maxloop;loop++)
330 mat1[i][j] = 2.0 - 4.0*invmaxlong*(double) rand();
332 printf("Original matrix:\n");
333 nr_printmat( mat1, n );
335 nr_copymat( mat1, n, mat2 );
337 i = nr_gaussj( mat2, n, 0, 0 );
339 if (i) printf("Singular matrix.\n");
341 printf("Inverted matrix:\n");
342 nr_printmat( mat2, n );
344 nr_multmat( mat1, n, mat2, mat3 );
346 printf("Original multiplied by inverse:\n");
347 nr_printmat( mat3, n );
349 if (loop < maxloop-1) /* sleep(1) */;
352 nr_free_matrix( mat1, 1, n, 1, n );
353 nr_free_matrix( mat2, 1, n, 1, n );
354 nr_free_matrix( mat3, 1, n, 1, n );