src/FDM/LaRCsim/ls_matrix.c

   1 /***************************************************************************
   2
   3         TITLE:          ls_matrix.c
   4
   5 ----------------------------------------------------------------------------
   6
   7         FUNCTION:       general real matrix routines; includes
   8                                 gaussj() for matrix inversion using
   9                                 Gauss-Jordan method with full pivoting.
  10
  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.
  22
  23 ----------------------------------------------------------------------------
  24
  25         MODULE STATUS:  developmental
  26
  27 ----------------------------------------------------------------------------
  28
  29         GENEALOGY:      Created 950222 E. B. Jackson
  30
  31 ----------------------------------------------------------------------------
  32
  33         DESIGNED BY:    from Numerical Recipes in C, by Press, et. al.
  34
  35         CODED BY:       Bruce Jackson
  36
  37         MAINTAINED BY:
  38
  39 ----------------------------------------------------------------------------
  40
  41         MODIFICATION HISTORY:
  42
  43         DATE    PURPOSE                                         BY
  44
  45         CURRENT RCS HEADER:
  46
  47 $Header$
  48 $Log$
  49 Revision 1.1  1999/06/17 18:07:34  curt
  50 Initial revision
  51
  52 Revision 1.1.1.1  1999/04/05 21:32:45  curt
  53 Start of 0.6.x branch.
  54
  55 Revision 1.1  1998/06/27 22:34:57  curt
  56 Initial revision.
  57
  58  * Revision 1.1  1995/02/27  19:55:44  bjax
  59  * Initial revision
  60  *
  61
  62 ----------------------------------------------------------------------------
  63
  64         REFERENCES:     [1] Press, William H., et. al, Numerical Recipes in
  65                             C, 2nd edition, Cambridge University Press, 1992
  66
  67 ----------------------------------------------------------------------------
  68
  69         CALLED BY:
  70
  71 ----------------------------------------------------------------------------
  72
  73         CALLS TO:
  74
  75 ----------------------------------------------------------------------------
  76
  77         INPUTS:
  78
  79 ----------------------------------------------------------------------------
  80
  81         OUTPUTS:
  82
  83 --------------------------------------------------------------------------*/
  84
  85 #ifdef HAVE_CONFIG_H
  86 #  include <config.h>
  87 #endif
  88
  89 #include <stdlib.h>
  90 #include <stdio.h>
  91 #include <math.h>
  92
  93 #ifdef HAVE_UNISTD_H
  94 #  include <unistd.h>
  95 #endif
  96
  97 #include "ls_matrix.h"
  98
  99
 100 #define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}
 101
 102 static char rcsid[] = "$Id$";
 103
 104
 105 int *nr_ivector(long nl, long nh)
 106 {
 107     int *v;
 108
 109     v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
 110     return v-nl+NR_END;
 111 }
 112
 113
 114
 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] */
 117 {
 118     long i, nrow=nrh-nrl+1, ncol=nch-ncl+1;
 119     double **m;
 120
 121     /* allocate pointers to rows */
 122     m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
 123
 124     if (!m)
 125         {
 126             fprintf(stderr, "Memory failure in routine 'nr_matrix'.\n");
 127             exit(1);
 128         }
 129
 130     m += NR_END;
 131     m -= nrl;
 132
 133     /* allocate rows and set pointers to them */
 134     m[nrl] = (double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
 135     if (!m[nrl])
 136         {
 137             fprintf(stderr, "Memory failure in routine 'matrix'\n");
 138             exit(1);
 139         }
 140
 141     m[nrl] += NR_END;
 142     m[nrl] -= ncl;
 143
 144     for (i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
 145
 146     /* return pointer to array of pointers to rows */
 147     return m;
 148 }
 149
 150
 151 void nr_free_ivector(int *v, long nl /* , long nh */)
 152 {
 153     free( (char *) (v+nl-NR_END));
 154 }
 155
 156
 157 void nr_free_matrix(double **m, long nrl, long nrh, long ncl, long nch)
 158 /* free a double matrix allocated by nr_matrix() */
 159 {
 160     free((char *) (m[nrl]+ncl-NR_END));
 161     free((char *) (m+nrl-NR_END));
 162 }
 163
 164
 165 int nr_gaussj(double **a, int n, double **b, int m)
 166
 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.                 */
 171
 172 /* Note: this routine modified by EBJ to make b optional, if m == 0 */
 173
 174 {
 175     int         *indxc, *indxr, *ipiv;
 176     int         i, icol, irow, j, k, l, ll;
 177     double       big, dum, pivinv, temp;
 178
 179     int         bexists = ((m != 0) || (b == 0));
 180
 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);
 184
 185     for (j=1;j<=n;j++) ipiv[j] = 0;
 186
 187     for (i=1;i<=n;i++)          /* This is the main loop over columns   */
 188         {
 189             big = 0.0;
 190             for (j=1;j<=n;j++)  /* This is outer loop of pivot search   */
 191                 if (ipiv[j] != 1)
 192                     for (k=1;k<=n;k++)
 193                         {
 194                             if (ipiv[k] == 0)
 195                                 {
 196                                     if (fabs(a[j][k]) >= big)
 197                                         {
 198                                             big = fabs(a[j][k]);
 199                                             irow = j;
 200                                             icol = k;
 201                                         }
 202                                 }
 203                             else
 204                                 if (ipiv[k] > 1) return -1;
 205                         }
 206             ++(ipiv[icol]);
 207
 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.            */
 216
 217             if (irow != icol)
 218                 {
 219 /*                  for (l=1;1<=n;l++) SWAP(a[irow][l],a[icol][l]) */
 220                         for (l=1;l<=n;l++)
 221                           {
 222                                 temp=a[irow][l];
 223                                 a[irow][l]=a[icol][l];
 224                                 a[icol][l]=temp;
 225                           }
 226                     if (bexists) for (l=1;l<=m;l++) SWAP(b[irow][l],b[icol][l])
 227                 }
 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];
 232             a[icol][icol] = 1.0;
 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  */
 237                     {
 238                         dum = a[ll][icol];
 239                         a[ll][icol] = 0.0;
 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;
 242                     }
 243         }
 244
 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. */
 249
 250     for (l=n;l>=1;l--)
 251         {
 252             if (indxr[l] != indxc[l])
 253                 for (k=1;k<=n;k++)
 254                     SWAP(a[k][indxr[l]],a[k][indxc[l]])
 255         }
 256
 257 /* and we are done */
 258
 259     nr_free_ivector(ipiv,1 /*,n*/ );
 260     nr_free_ivector(indxr,1 /*,n*/ );
 261     nr_free_ivector(indxc,1 /*,n*/ );
 262
 263     return 0;   /* indicate success */
 264 }
 265
 266 void nr_copymat(double **orig, int n, double **copy)
 267 /* overwrites matrix 'copy' with copy of matrix 'orig' */
 268 {
 269         long i, j;
 270
 271         if ((orig==0)||(copy==0)||(n==0)) return;
 272
 273         for (i=1;i<=n;i++)
 274                 for (j=1;j<=n;j++)
 275                         copy[i][j] = orig[i][j];
 276 }
 277
 278 void nr_multmat(double **m1, int n, double **m2, double **prod)
 279 {
 280         long i, j, k;
 281
 282         if ((m1==0)||(m2==0)||(prod==0)||(n==0)) return;
 283
 284         for (i=1;i<=n;i++)
 285                 for (j=1;j<=n;j++)
 286                         {
 287                                 prod[i][j] = 0.0;
 288                                 for(k=1;k<=n;k++) prod[i][j] += m1[i][k]*m2[k][j];
 289                         }
 290 }
 291
 292
 293
 294 void nr_printmat(double **a, int n)
 295 {
 296     int i,j;
 297
 298     printf("\n");
 299     for(i=1;i<=n;i++)
 300       {
 301           for(j=1;j<=n;j++)
 302               printf("% 9.4f ", a[i][j]);
 303           printf("\n");
 304       }
 305     printf("\n");
 306
 307 }
 308
 309
 310 void testmat( void ) /* main() for test purposes */
 311 {
 312     double **mat1, **mat2, **mat3;
 313     double invmaxlong;
 314     int loop, i, j, n = 20;
 315     long maxlong = RAND_MAX;
 316     int maxloop = 2;
 317
 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 );
 322
 323 /*    for(i=1;i<=n;i++) mat1[i][i]= 5.0; */
 324
 325         for(loop=0;loop<maxloop;loop++)
 326         {
 327             if (loop != 0)
 328                 for(i=1;i<=n;i++)
 329                     for(j=1;j<=n;j++)
 330                         mat1[i][j] = 2.0 - 4.0*invmaxlong*(double) rand();
 331
 332                 printf("Original matrix:\n");
 333             nr_printmat( mat1, n );
 334
 335             nr_copymat( mat1, n, mat2 );
 336
 337             i = nr_gaussj( mat2, n, 0, 0 );
 338
 339             if (i) printf("Singular matrix.\n");
 340
 341                 printf("Inverted matrix:\n");
 342             nr_printmat( mat2, n );
 343
 344             nr_multmat( mat1, n, mat2, mat3 );
 345
 346             printf("Original multiplied by inverse:\n");
 347             nr_printmat( mat3, n );
 348
 349             if (loop < maxloop-1) /* sleep(1) */;
 350         }
 351
 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 );
 355 }