Go to the documentation of this file.00001
00002 #include <math.h>
00003 #include <grass/libtrans.h>
00004
00005 #define EPSILON 1.0e-16
00006
00007
00008 #define N DIM_matrix
00009
00010
00011
00012
00013
00014
00015
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
00030 for (i = 0; i < N; i++)
00031 ipivot[i] = 0;
00032
00033 for (i = 0; i < N; i++) {
00034 t = 0.0;
00035
00036 for (j = 0; j < N; j++) {
00037 if (ipivot[j] == 1)
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:
00055 return (-1);
00056 break;
00057 }
00058 }
00059
00060 ipivot[ic] += 1;
00061 if (ipivot[ic] > 1) {
00062 return (-1);
00063 }
00064
00065
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
00078 if (fabs(pivot[i]) < EPSILON) {
00079 return (-1);
00080 }
00081
00082
00083 m[ic][ic] = 1.0;
00084
00085 for (j = 0; j < N; j++)
00086 m[ic][j] /= pivot[i];
00087
00088
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
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
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 }