GRASS Programmer's Manual  6.5.svn(2014)-r66266
inverse.c
Go to the documentation of this file.
1 /* @(#)inverse.c 2.1 6/26/87 */
2 #include <math.h>
3 #include <grass/libtrans.h>
4
5 #define EPSILON 1.0e-16
6
7 /* DIM_matrix is defined in "libtrans.h" */
8 #define N DIM_matrix
9
10 /*
11  * inverse: invert a square matrix (puts pivot elements on main diagonal).
12  * returns arg2 as the inverse of arg1.
13  *
14  * This routine is based on a routine found in Andrei Rogers, "Matrix
15  * Methods in Urban and Regional Analysis", (1971), pp. 143-153.
16  */
17 int inverse(double m[N][N])
18 {
19  int i, j, k, l, ir = 0, ic = 0;
20  int ipivot[N], itemp[N][2];
21  double pivot[N], t;
22  double fabs();
23
24
25  if (isnull(m))
26  return (-1);
27
28
29  /* initialization */
30  for (i = 0; i < N; i++)
31  ipivot[i] = 0;
32
33  for (i = 0; i < N; i++) {
34  t = 0.0; /* search for pivot element */
35
36  for (j = 0; j < N; j++) {
37  if (ipivot[j] == 1) /* found pivot */
38  continue;
39
40  for (k = 0; k < N; k++)
41  switch (ipivot[k] - 1) {
42  case 0:
43  break;
44  case -1:
45  if (fabs(t) < fabs(m[j][k])) {
46  ir = j;
47  ic = k;
48  t = m[j][k];
49  }
50  break;
51  case 1:
52  return (-1);
53  break;
54  default: /* shouldn't get here */
55  return (-1);
56  break;
57  }
58  }
59
60  ipivot[ic] += 1;
61  if (ipivot[ic] > 1) { /* check for dependency */
62  return (-1);
63  }
64
65  /* interchange rows to put pivot element on diagonal */
66  if (ir != ic)
67  for (l = 0; l < N; l++) {
68  t = m[ir][l];
69  m[ir][l] = m[ic][l];
70  m[ic][l] = t;
71  }
72
73  itemp[i][0] = ir;
74  itemp[i][1] = ic;
75  pivot[i] = m[ic][ic];
76
77  /* check for zero pivot */
78  if (fabs(pivot[i]) < EPSILON) {
79  return (-1);
80  }
81
82  /* divide pivot row by pivot element */
83  m[ic][ic] = 1.0;
84
85  for (j = 0; j < N; j++)
86  m[ic][j] /= pivot[i];
87
88  /* reduce nonpivot rows */
89  for (k = 0; k < N; k++)
90  if (k != ic) {
91  t = m[k][ic];
92  m[k][ic] = 0.0;
93
94  for (l = 0; l < N; l++)
95  m[k][l] -= (m[ic][l] * t);
96  }
97  }
98
99  /* interchange columns */
100  for (i = 0; i < N; i++) {
101  l = N - i - 1;
102  if (itemp[l][0] == itemp[l][1])
103  continue;
104
105  ir = itemp[l][0];
106  ic = itemp[l][1];
107
108  for (k = 0; k < N; k++) {
109  t = m[k][ir];
110  m[k][ir] = m[k][ic];
111  m[k][ic] = t;
112  }
113  }
114
115  return 1;
116 }
117
118
119
120
121 #define ZERO 1.0e-8
122
123 /*
124  * isnull: returns 1 if matrix is null, else 0.
125  */
126
127 int isnull(double a[N][N])
128 {
129  register int i, j;
130  double fabs();
131
132
133  for (i = 0; i < N; i++)
134  for (j = 0; j < N; j++)
135  if ((fabs(a[i][j]) - ZERO) > ZERO)
136  return 0;
137
138  return 1;
139 }
int l