File: matinv.x

package info (click to toggle)
iraf-rvsao 2.8.3-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, buster, sid
  • size: 16,456 kB
  • sloc: ansic: 963; lisp: 651; fortran: 397; makefile: 27
file content (97 lines) | stat: -rw-r--r-- 1,934 bytes parent folder | download
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
# File rvsao/Util/matinv.x
# October 29, 1991
# By Doug Mink, Harvard-Smithsonian Center for Astrophysics
# After Bevington, pages 302-303.

#  Invert a symmetric matrix up to 10x10 and calculate its determinant

procedure matinv (array,norder,det)

double	array[16,16]	# Input matrix which is replaced by its inverse
int	norder		# Degree of matrix (order of determinant)
double	det		# Determinant of input matrix

double	amax,save
int	ik[10],jk[10]
int	i,j,k,l

begin
	det = 1.d0
	do k = 1,norder {

# find largest element array[i,j] in rest of matrix
	    amax = 0. 
21	    do i = k, norder {
		do j = k, norder {
		    if (dabs(amax) <= dabs(array[i,j])) {
			amax = array[i,j] 
			ik[k] = i 
			jk[k] = j 
			}
		    }
		}

# interchange rows and columns to put amax in array[k,k]
	    if (amax == 0) {
		det = 0.
		return
		}
	    i = ik[k] 
	    if (i < k) go to 21
	    if (i > k) {
		do j = 1, norder {
		    save = array[k,j] 
		    array[k,j] = array[i,j]
		    array[i,j] = -save
		    }
		}
	    j = jk[k] 
	    if (j < k) go to 21
	    if (j > k) {
		do i = 1, norder {
		    save = array[i,k] 
		    array[i,k] = array[i,j]
		    array[i,j] = -save
		    }
		}

# accumulate elements of inverse matrix 
	    do i = 1,norder {
		if (i != k) array[i,k] = -array[i,k]/amax
		}
	    do i = 1, norder {
		do j = 1, norder {
		    if (i != k && j != k) 
			array[i,j] = array[i,j]+array[i,k]*array[k,j]
		    }
		}
	    do j = 1,norder {
		if (j != k) array[k,j] = array[k,j]/amax
		}
	    array[k,k] = 1./amax
	    det = det*amax
	    }

#  Restore ordering of matrix
	do l = 1, norder {
	    k = norder-l+1
	    j = ik[k] 
	    if (j > k) {
		do i = 1, norder {
		    save = array[i,k] 
		    array[i,k] = -array[i,j]
		    array[i,j] = save 
		    }
		}
	    i = jk[k] 
	    if (i > k) {
		do j = 1, norder {
		    save = array[k,j] 
		    array[k,j] = -array[i,j]
		    array[i,j] = save 
		    }
		}
	    }

	return
	end