1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
|
/* minv.c
*
* Matrix inversion
*
*
*
* SYNOPSIS:
*
* int n, errcod;
* double A[n*n], X[n*n];
* double B[n];
* int IPS[n];
* int minv();
*
* errcod = minv( A, X, n, B, IPS );
*
*
*
* DESCRIPTION:
*
* Finds the inverse of the n by n matrix A. The result goes
* to X. B and IPS are scratch pad arrays of length n.
* The contents of matrix A are destroyed.
*
* The routine returns nonzero on error; error messages are printed
* by subroutine simq().
*
*/
minv( A, X, n, B, IPS )
double A[], X[];
int n;
double B[];
int IPS[];
{
double *pX;
int i, j, k;
for( i=1; i<n; i++ )
B[i] = 0.0;
B[0] = 1.0;
/* Reduce the matrix and solve for first right hand side vector */
pX = X;
k = simq( A, B, pX, n, 1, IPS );
if( k )
return(-1);
/* Solve for the remaining right hand side vectors */
for( i=1; i<n; i++ )
{
B[i-1] = 0.0;
B[i] = 1.0;
pX += n;
k = simq( A, B, pX, n, -1, IPS );
if( k )
return(-1);
}
/* Transpose the array of solution vectors */
mtransp( n, X, X );
return(0);
}
|
