inverse.c

Go to the documentation of this file.
00001 /*  @(#)inverse.c       2.1  6/26/87  */
00002 #include <math.h>
00003 #include <grass/libtrans.h>
00004 
00005 #define EPSILON 1.0e-16
00006 
00007 /* DIM_matrix is defined in "libtrans.h" */
00008 #define N       DIM_matrix
00009 
00010 /*
00011  * inverse: invert a square matrix (puts pivot elements on main diagonal).
00012  *          returns arg2 as the inverse of arg1.
00013  *
00014  *  This routine is based on a routine found in Andrei Rogers, "Matrix
00015  *  Methods in Urban and Regional Analysis", (1971), pp. 143-153.
00016  */
00017 int inverse(double m[N][N])
00018 {
00019     int i, j, k, l, ir = 0, ic = 0;
00020     int ipivot[N], itemp[N][2];
00021     double pivot[N], t;
00022     double fabs();
00023 
00024 
00025     if (isnull(m))
00026         return (-1);
00027 
00028 
00029     /* initialization */
00030     for (i = 0; i < N; i++)
00031         ipivot[i] = 0;
00032 
00033     for (i = 0; i < N; i++) {
00034         t = 0.0;                /* search for pivot element */
00035 
00036         for (j = 0; j < N; j++) {
00037             if (ipivot[j] == 1) /* found pivot */
00038                 continue;
00039 
00040             for (k = 0; k < N; k++)
00041                 switch (ipivot[k] - 1) {
00042                 case 0:
00043                     break;
00044                 case -1:
00045                     if (fabs(t) < fabs(m[j][k])) {
00046                         ir = j;
00047                         ic = k;
00048                         t = m[j][k];
00049                     }
00050                     break;
00051                 case 1:
00052                     return (-1);
00053                     break;
00054                 default:        /* shouldn't get here */
00055                     return (-1);
00056                     break;
00057                 }
00058         }
00059 
00060         ipivot[ic] += 1;
00061         if (ipivot[ic] > 1) {   /* check for dependency */
00062             return (-1);
00063         }
00064 
00065         /* interchange rows to put pivot element on diagonal */
00066         if (ir != ic)
00067             for (l = 0; l < N; l++) {
00068                 t = m[ir][l];
00069                 m[ir][l] = m[ic][l];
00070                 m[ic][l] = t;
00071             }
00072 
00073         itemp[i][0] = ir;
00074         itemp[i][1] = ic;
00075         pivot[i] = m[ic][ic];
00076 
00077         /* check for zero pivot */
00078         if (fabs(pivot[i]) < EPSILON) {
00079             return (-1);
00080         }
00081 
00082         /* divide pivot row by pivot element */
00083         m[ic][ic] = 1.0;
00084 
00085         for (j = 0; j < N; j++)
00086             m[ic][j] /= pivot[i];
00087 
00088         /* reduce nonpivot rows */
00089         for (k = 0; k < N; k++)
00090             if (k != ic) {
00091                 t = m[k][ic];
00092                 m[k][ic] = 0.0;
00093 
00094                 for (l = 0; l < N; l++)
00095                     m[k][l] -= (m[ic][l] * t);
00096             }
00097     }
00098 
00099     /* interchange columns */
00100     for (i = 0; i < N; i++) {
00101         l = N - i - 1;
00102         if (itemp[l][0] == itemp[l][1])
00103             continue;
00104 
00105         ir = itemp[l][0];
00106         ic = itemp[l][1];
00107 
00108         for (k = 0; k < N; k++) {
00109             t = m[k][ir];
00110             m[k][ir] = m[k][ic];
00111             m[k][ic] = t;
00112         }
00113     }
00114 
00115     return 1;
00116 }
00117 
00118 
00119 
00120 
00121 #define ZERO 1.0e-8
00122 
00123 /*
00124  * isnull: returns 1 if matrix is null, else 0.
00125  */
00126 
00127 int isnull(double a[N][N])
00128 {
00129     register int i, j;
00130     double fabs();
00131 
00132 
00133     for (i = 0; i < N; i++)
00134         for (j = 0; j < N; j++)
00135             if ((fabs(a[i][j]) - ZERO) > ZERO)
00136                 return 0;
00137 
00138     return 1;
00139 }

Generated on Sat Oct 24 03:25:20 2009 for GRASS Programmer's Manual by  doxygen 1.6.1