forked from leto/parrot-lapack
/
inverse.winxed
99 lines (71 loc) · 1.95 KB
/
inverse.winxed
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
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
//DGETRI computes the inverse of a matrix using the LU factorization computed by DGETRF.
//This method inverts U and then computes inv(A) by solving the system
//inv(A)*L = inv(U) for inv(A).
$load "Double.pbc";
namespace dgetri_func{
const int PRINT_DEBUG_STUFF = 0;
//a is DOUBLE PRECISION array, dimension (LDA,N)
function dgetri_exec(var a)
{
var pla = loadlib("linalg_group");
var lapack = loadlib('liblapack.so');
if (lapack == null || !lapack)
die("Cannot find liblapack");
say("Given Matrix:");
say(a);
/*
A is DOUBLE PRECISION array, dimension (LDA,N)
On entry, the factors L and U from the factorization
A = P*L*U as computed by DGETRF.
On exit, if INFO = 0, the inverse of the original matrix A.
*/
int n,m,ipiv_size;
m=a.rows;
n=a.cols;
ipiv_size=min(m,n);
int ipiv[n];
int i;
for(i=0;i<n;i++)
ipiv[i]=0;
var dgetri = dlfunc(lapack, "dgetri_", "vppppppp");
if(dgetri == null || !dgetri)
die("Not DGETRI");
int lda,lwork,info;
lda=max(1,m);
lda=max(lda,n);
lwork=max(1,n);
// say(n);
// say(m);
// say(lda);
// say("lwork=",lwork);
int work[lwork];
for(i=0;i<n;++i)
work[i]=0;
//calling the function
//This method inverts U and then computes inv(A) by solving the system inv(A)*L = inv(U) for inv(A).
dgetri(n,a,lda,ipiv,work,lwork,info);
/*
INFO is INTEGER
= 0: successful exit
< 0: if INFO = -i, the i-th argument had an illegal value
> 0: if INFO = i, U(i,i) is exactly zero; the matrix is
singular and its inverse could not be computed.
*/
return info;
}
function max(var a,var b)
{
return a>b?a:b;
}
function min(var a,var b)
{
return a<b?a:b;
}
function debug(var matrix, string msg, var args [slurpy])
{
if (PRINT_DEBUG_STUFF) {
say(sprintf(msg, args));
say(matrix);
}
}
}