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
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
# MW_INVERTD -- Invert a square matrix, double precision version. The matrix
# need not be symmetric. The input and output matrices cannot be the same.
procedure mw_invertd (o_ltm, n_ltm, ndim)
double o_ltm[ndim,ndim] #I input matrix
double n_ltm[ndim,ndim] #O output (inverted) matrix
int ndim #I dimensionality of system
pointer sp, ix, ltm
int nelem, i, j
begin
call smark (sp)
nelem = ndim * ndim
call salloc (ix, ndim, TY_INT)
call salloc (ltm, nelem, TY_DOUBLE)
# Make scratch copy (to be modified) of input matrix.
call amovd (o_ltm, Memd[ltm], nelem)
# Set up identity matrix.
do i = 1, ndim {
do j = 1, ndim
n_ltm[i,j] = 0.0
n_ltm[i,i] = 1.0
}
# Perform the LU decomposition.
call mw_ludecompose (Memd[ltm], Memi[ix], ndim)
# Compute the inverse matrix by backsubstitution.
do j = 1, ndim
call mw_lubacksub (Memd[ltm], Memi[ix], n_ltm[1,j], ndim)
call sfree (sp)
end
|
