GRASS Programmer's Manual
6.4.1(2011)
|
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 }