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/04/05 21:32:45  curt
  50 Initial revision
  51
  52 Revision 1.1  1998/06/27 22:34:57  curt
  53 Initial revision.
  54
  55  * Revision 1.1  1995/02/27  19:55:44  bjax
  56  * Initial revision
  57  *
  58
  59 ----------------------------------------------------------------------------
  60
  61         REFERENCES:     [1] Press, William H., et. al, Numerical Recipes in
  62                             C, 2nd edition, Cambridge University Press, 1992
  63
  64 ----------------------------------------------------------------------------
  65
  66         CALLED BY:
  67
  68 ----------------------------------------------------------------------------
  69
  70         CALLS TO:
  71
  72 ----------------------------------------------------------------------------
  73
  74         INPUTS:
  75
  76 ----------------------------------------------------------------------------
  77
  78         OUTPUTS:
  79
  80 --------------------------------------------------------------------------*/
  81
  82 #ifdef HAVE_CONFIG_H
  83 #  include <config.h>
  84 #endif
  85
  86 #include <stdlib.h>
  87 #include <stdio.h>
  88 #include <math.h>
  89
  90 #ifdef HAVE_UNISTD_H
  91 #  include <unistd.h>
  92 #endif
  93
  94 #include "ls_matrix.h"
  95
  96
  97 #define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}
  98
  99 static char rcsid[] = "$Id$";
 100
 101
 102 int *nr_ivector(long nl, long nh)
 103 {
 104     int *v;
 105
 106     v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
 107     return v-nl+NR_END;
 108 }
 109
 110
 111
 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] */
 114 {
 115     long i, nrow=nrh-nrl+1, ncol=nch-ncl+1;
 116     double **m;
 117
 118     /* allocate pointers to rows */
 119     m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
 120
 121     if (!m)
 122         {
 123             fprintf(stderr, "Memory failure in routine 'nr_matrix'.\n");
 124             exit(1);
 125         }
 126
 127     m += NR_END;
 128     m -= nrl;
 129
 130     /* allocate rows and set pointers to them */
 131     m[nrl] = (double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
 132     if (!m[nrl])
 133         {
 134             fprintf(stderr, "Memory failure in routine 'matrix'\n");
 135             exit(1);
 136         }
 137
 138     m[nrl] += NR_END;
 139     m[nrl] -= ncl;
 140
 141     for (i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
 142
 143     /* return pointer to array of pointers to rows */
 144     return m;
 145 }
 146
 147
 148 void nr_free_ivector(int *v, long nl /* , long nh */)
 149 {
 150     free( (char *) (v+nl-NR_END));
 151 }
 152
 153
 154 void nr_free_matrix(double **m, long nrl, long nrh, long ncl, long nch)
 155 /* free a double matrix allocated by nr_matrix() */
 156 {
 157     free((char *) (m[nrl]+ncl-NR_END));
 158     free((char *) (m+nrl-NR_END));
 159 }
 160
 161
 162 int nr_gaussj(double **a, int n, double **b, int m)
 163
 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.                 */
 168
 169 /* Note: this routine modified by EBJ to make b optional, if m == 0 */
 170
 171 {
 172     int         *indxc, *indxr, *ipiv;
 173     int         i, icol, irow, j, k, l, ll;
 174     double       big, dum, pivinv, temp;
 175
 176     int         bexists = ((m != 0) || (b == 0));
 177
 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);
 181
 182     for (j=1;j<=n;j++) ipiv[j] = 0;
 183
 184     for (i=1;i<=n;i++)          /* This is the main loop over columns   */
 185         {
 186             big = 0.0;
 187             for (j=1;j<=n;j++)  /* This is outer loop of pivot search   */
 188                 if (ipiv[j] != 1)
 189                     for (k=1;k<=n;k++)
 190                         {
 191                             if (ipiv[k] == 0)
 192                                 {
 193                                     if (fabs(a[j][k]) >= big)
 194                                         {
 195                                             big = fabs(a[j][k]);
 196                                             irow = j;
 197                                             icol = k;
 198                                         }
 199                                 }
 200                             else
 201                                 if (ipiv[k] > 1) return -1;
 202                         }
 203             ++(ipiv[icol]);
 204
 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.            */
 213
 214             if (irow != icol)
 215                 {
 216 /*                  for (l=1;1<=n;l++) SWAP(a[irow][l],a[icol][l]) */
 217                         for (l=1;l<=n;l++)
 218                           {
 219                                 temp=a[irow][l];
 220                                 a[irow][l]=a[icol][l];
 221                                 a[icol][l]=temp;
 222                           }
 223                     if (bexists) for (l=1;l<=m;l++) SWAP(b[irow][l],b[icol][l])
 224                 }
 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];
 229             a[icol][icol] = 1.0;
 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  */
 234                     {
 235                         dum = a[ll][icol];
 236                         a[ll][icol] = 0.0;
 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;
 239                     }
 240         }
 241
 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. */
 246
 247     for (l=n;l>=1;l--)
 248         {
 249             if (indxr[l] != indxc[l])
 250                 for (k=1;k<=n;k++)
 251                     SWAP(a[k][indxr[l]],a[k][indxc[l]])
 252         }
 253
 254 /* and we are done */
 255
 256     nr_free_ivector(ipiv,1 /*,n*/ );
 257     nr_free_ivector(indxr,1 /*,n*/ );
 258     nr_free_ivector(indxc,1 /*,n*/ );
 259
 260     return 0;   /* indicate success */
 261 }
 262
 263 void nr_copymat(double **orig, int n, double **copy)
 264 /* overwrites matrix 'copy' with copy of matrix 'orig' */
 265 {
 266         long i, j;
 267
 268         if ((orig==0)||(copy==0)||(n==0)) return;
 269
 270         for (i=1;i<=n;i++)
 271                 for (j=1;j<=n;j++)
 272                         copy[i][j] = orig[i][j];
 273 }
 274
 275 void nr_multmat(double **m1, int n, double **m2, double **prod)
 276 {
 277         long i, j, k;
 278
 279         if ((m1==0)||(m2==0)||(prod==0)||(n==0)) return;
 280
 281         for (i=1;i<=n;i++)
 282                 for (j=1;j<=n;j++)
 283                         {
 284                                 prod[i][j] = 0.0;
 285                                 for(k=1;k<=n;k++) prod[i][j] += m1[i][k]*m2[k][j];
 286                         }
 287 }
 288
 289
 290
 291 void nr_printmat(double **a, int n)
 292 {
 293     int i,j;
 294
 295     printf("\n");
 296     for(i=1;i<=n;i++)
 297       {
 298           for(j=1;j<=n;j++)
 299               printf("% 9.4f ", a[i][j]);
 300           printf("\n");
 301       }
 302     printf("\n");
 303
 304 }
 305
 306
 307 void testmat( void ) /* main() for test purposes */
 308 {
 309     double **mat1, **mat2, **mat3;
 310     double invmaxlong;
 311     int loop, i, j, n = 20;
 312     long maxlong = RAND_MAX;
 313     int maxloop = 2;
 314
 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 );
 319
 320 /*    for(i=1;i<=n;i++) mat1[i][i]= 5.0; */
 321
 322         for(loop=0;loop<maxloop;loop++)
 323         {
 324             if (loop != 0)
 325                 for(i=1;i<=n;i++)
 326                     for(j=1;j<=n;j++)
 327                         mat1[i][j] = 2.0 - 4.0*invmaxlong*(double) rand();
 328
 329                 printf("Original matrix:\n");
 330             nr_printmat( mat1, n );
 331
 332             nr_copymat( mat1, n, mat2 );
 333
 334             i = nr_gaussj( mat2, n, 0, 0 );
 335
 336             if (i) printf("Singular matrix.\n");
 337
 338                 printf("Inverted matrix:\n");
 339             nr_printmat( mat2, n );
 340
 341             nr_multmat( mat1, n, mat2, mat3 );
 342
 343             printf("Original multiplied by inverse:\n");
 344             nr_printmat( mat3, n );
 345
 346             if (loop < maxloop-1) /* sleep(1) */;
 347         }
 348
 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 );
 352 }