GRASS Programmer's Manual
6.4.2(2012)
Main Page
Related Pages
Namespaces
Data Structures
Files
File List
Globals
All
Data Structures
Namespaces
Files
Functions
Variables
Typedefs
Enumerations
Enumerator
Macros
Pages
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
}
lib
vector
transform
inverse.c
Generated on Sun Sep 9 2012 18:55:33 for GRASS Programmer's Manual by
1.8.1.2