File: invert.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 (130 lines) | stat: -rw-r--r-- 2,668 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
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
# File rvsao/Util/invert.x
# December 2, 1991
# By Doug Mink, Harvard-Smithsonian Center for Astrophysics
# After John Tonry (9/2/80)

# Copyright(c) 1991 Smithsonian Astrophysical Observatory
# You may do anything you like with this file except remove this copyright.
# The Smithsonian Astrophysical Observatory makes no representations about
# the suitability of this software for any purpose.  It is provided "as is"
# without express or implied warranty.

procedure invert (n,a,rst,det)

int	n		# Dimension of the square matrix
double	a[n,n]		# Matrix to be inverted.
			# Upon successful inversion, a contains the inverse.
int	rst[2,n]	# Scratch vector (the row status vector)
double	det		# Determinant of the matrix
			# Set to 0 for a singular matrix, a then is garbage.


double	save, pivot, onrow, cprev, cnow, decr
int	mrank,isign,i,j,k,l,nrow,ncol

begin
	mrank = 0
	isign = 1
	det = 0.d0
	do j = 1, n {
	    do i = 1, 2 {
		rst[i,j] = 0
		}
	    }

#	Loop over columns, reducing each
	do i = 1, n {

#	Find the pivot element
	    pivot = 0
	    nrow = 0
	    ncol = 0
	    do j = 1, n {
		if (rst[1,j] != 0) next
		do k = 1, n {
		    if (rst[1,k] != 0) next
		    if (pivot >= dabs (a[j,k])) next
		    pivot = dabs (a[j,k])
		    nrow = j
		    ncol = k
		    }
		}
	    pivot = a[nrow,ncol]
	    if (pivot == 0) {
		det = 0
		return
		}
	    rst[1,ncol] = nrow
	    rst[2,ncol] = i

#	Swap pivot element onto the diagonal
	    do k = 1, n {
		save = a[nrow,k]
		a[nrow,k] = a[ncol,k]
		a[ncol,k] = save
		}

#	Reduce pivot column
	    do j = 1, n {
		a[j,ncol] = -a[j,ncol] / pivot
		}
	    a[ncol,ncol] = 1 / pivot

#	Reduce other columns
	    do k = 1, n {
		if (k == ncol) next

#	Find maximum of column to check for singularity
		cprev = 0
		do j = 1, n {
		    cprev = dmax1(cprev,dabs(a[j,k]))
		    }

#	Reduce the column
		onrow = a[ncol,k]
		a[ncol,k] = 0
		do j = 1, n {
		    a[j,k] = a[j,k] + onrow * a[j,ncol]
		    }

#	Find the new maximum of the column
		cnow = 0
		do j = 1, n {
		    cnow = dmax1 (cnow,dabs(a[j,k]))
		    }

#	Quit if too many figures accuracy were lost (singular)
		if (cnow == 0)  {
		    det = 0
		    return
		    }
		decr = cprev / cnow
		if (decr > 1.e8) {
		    det = 0
		    return
		    }

		}

	    det = det + dlog10 (dabs (pivot))
	    if (pivot < 0) isign = -isign
	    mrank = mrank + 1
	    }

#     now untangle the mess
	do j = 1, n {
	    do k = 1, n {
		if (rst[2,k] != (n + 1 - j)) next
		ncol = rst[1,k]
		if(ncol == k) break
		do l = 1, n {
		    save = a[l,ncol]
		    a[l,ncol] = a[l,k]
		    a[l,k] = save
		    }
		break
		}
	    }
	if (abs(det) < 35) det = isign * (10. ** det)
	return
end