1 //-----------------------------------------------------------------------------
3 // Input: Matrix on stack
5 // Output: Inverse on stack
10 // ((-2,1),(3/2,-1/2))
14 // Uses Gaussian elimination for numerical matrices.
16 //-----------------------------------------------------------------------------
26 else if (p1
->u
.tensor
->ndim
!= 2)
28 else if (p1
->u
.tensor
->dim
[0] != p1
->u
.tensor
->dim
[1])
44 if (check_arg() == 0) {
52 n
= p1
->u
.tensor
->nelem
;
54 a
= p1
->u
.tensor
->elem
;
56 for (i
= 0; i
< n
; i
++)
69 stop("inverse of singular matrix");
84 if (check_arg() == 0) {
97 // inverse using gaussian elimination
104 n
= p1
->u
.tensor
->dim
[0];
108 for (i
= 0; i
< n
; i
++)
109 for (j
= 0; j
< n
; j
++)
115 for (i
= 0; i
< n
* n
; i
++)
116 push(p1
->u
.tensor
->elem
[i
]);
120 p1
= alloc_tensor(n
* n
);
122 p1
->u
.tensor
->ndim
= 2;
123 p1
->u
.tensor
->dim
[0] = n
;
124 p1
->u
.tensor
->dim
[1] = n
;
126 for (i
= 0; i
< n
* n
; i
++)
127 p1
->u
.tensor
->elem
[i
] = stack
[h
+ i
];
134 //-----------------------------------------------------------------------------
136 // Input: n * n unit matrix on stack
138 // n * n operand on stack
140 // Output: n * n inverse matrix on stack
142 // n * n garbage on stack
146 //-----------------------------------------------------------------------------
148 #define A(i, j) stack[a + n * (i) + (j)]
149 #define U(i, j) stack[u + n * (i) + (j)]
160 for (d
= 0; d
< n
; d
++) {
162 // diagonal element zero?
164 if (equal(A(d
, d
), zero
)) {
168 for (i
= d
+ 1; i
< n
; i
++)
169 if (!equal(A(i
, d
), zero
))
173 stop("inverse of singular matrix");
177 for (j
= 0; j
< n
; j
++) {
189 // multiply the pivot row by 1 / pivot
193 for (j
= 0; j
< n
; j
++) {
208 // clear out the column above and below the pivot
210 for (i
= 0; i
< n
; i
++) {
219 // add pivot row to i-th row
221 for (j
= 0; j
< n
; j
++) {