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
62
63
64
65
66
67
68
69
|
/***************************************************
*
* routine for inverting a matrix, extracted from
* Ziheng Yang's source files.
*
****************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <math.h>
#include <limits.h>
#include <float.h>
#define PYCOGENT_VERSION "1.4.1"
/*****
*
* taking inv of matrix
*
*
******/
int matinv (double x[], int n, int m, double space[])
{
/* x[n*m] ... m>=n
*/
register int i,j,k;
int *irow=(int*) space;
double ee=1.0e-30, t,t1,xmax;
double det=1.0;
for (i = 0; i < n; i++) {
xmax = 0.;
for (j= i; j < n; j++)
if (xmax < fabs(x[j*m+i])) { xmax = fabs(x[j*m+i]); irow[i]=j; }
det *= xmax;
if (xmax < ee) {
/* printf("\nDet = %.4e close to zero at %3d!\t\n", xmax,i+1); */
return(-1);
}
if (irow[i] != i) {
for (j = 0; j < m; j++) {
t = x[i*m+j];
x[i*m+j] = x[irow[i]*m+j];
x[irow[i]*m+j] = t;
}
}
t = 1./x[i*m+i];
for (j = 0; j < n; j++) {
if (j == i) continue;
t1 = t*x[j*m+i];
for(k = 0; k < m; k++) x[j*m+k] -= t1*x[i*m+k];
x[j*m+i] = -t1;
}
for(j = 0; j < m; j++) x[i*m+j] *= t;
x[i*m+i] = t;
} /* i */
for (i=n-1; i>=0; i--) {
if (irow[i] == i) continue;
for(j = 0; j < n; j++) {
t = x[j*m+i];
x[j*m+i] = x[j*m + irow[i]];
x[j*m + irow[i]] = t;
}
}
return (0);
}
|
