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
|
C Copyright 1981-2007 ECMWF
C
C Licensed under the GNU Lesser General Public License which
C incorporates the terms and conditions of version 3 of the GNU
C General Public License.
C See LICENSE and gpl-3.0.txt for details.
C
SUBROUTINE SYMINV(A,NDIM,N,COND,V,DEPS)
DIMENSION A(NDIM,NDIM),V(NDIM)
C INVERT IN PLACE THE LOWER TRIANGLE OF A (I.E. A(I,J) I.GE.J)
C A IS A SYMMETRIC POSITIVE DEFINITE MATRIX
C THE UPPER TRIANGLE OF A (IE. A(I,J) I.LT.J) IS NOT USED OR ALTERED
C
C THIS VERSION IS OPTIMIZED FOR THE CDC FTN COMPILER.
C IT REQUIRES THE FUNCTION 'NUMARG' FROM 'ECLIB'
C
IF(N.LT.1)RETURN
EPS=0.
c IF(NUMARG().EQ.6) EPS=DEPS
NDIA=NDIM+1
IMX=ISAMAX(N,A(1,1),NDIA)
ZMX=A(IMX,IMX)*FLOAT(N*N)
C SECTION 1 & 2 COMBINED IN THIS VERSION
C**** 1. ROOT FREE CHOLESKY DECOMPOSITION A =L D L(TRANSPOSE)
J=1
IF(A(1,1).LE.EPS)GOTO91
V(1)=1./A(1,1)
IF(N.EQ.1)GOTO20
X=A(2,1)*V(1)
A(2,2)=A(2,2)-X*A(2,1)
A(2,1)=X
IF(A(2,2).LE.EPS)GOTO91
V(2)=1./A(2,2)
IF(N.EQ.2)GOTO20
DO 14 I=3,N
DO 12 J=3,I
S=A(I,J-1)
DO 11 K=3,J
11 S=S-A(I,K-2)*A(J-1,K-2)
A(I,J-1)=S
12 CONTINUE
S=A(I,I)
DO 13 J=2,I
X=A(I,J-1)*V(J-1)
S=S-X*A(I,J-1)
13 A(I,J-1)=X
A(I,I)=S
C CHECK FOR POSITIVE-DEFINITENESS AND INVERT DIAGONAL MATRIX D.
IF(A(I,I).LE.EPS)GOTO91
V(I)=1./A(I,I)
14 CONTINUE
C
C**** 2. COPY INVERSE OF D WHICH HAS ALREADY BEEN CALCULATED.
20 DO 21 J=1,N
21 A(J,J)=V(J)
IF(N.EQ.1) GO TO 50
C
C**** 3. INVERSION OF L
30 A(2,1)=-A(2,1)
NM1=N-1
IF(N.EQ.2)GOTO40
DO 33 I=2,NM1
DO 32 J=2,I
S=A(I+1,J-1)
DO 31 K=J,I
31 S=S+A(I+1,K)*A(K,J-1)
32 A(I+1,J-1)=-S
33 A(I+1,I)=-A(I+1,I)
C
C**** 4. INV A = INV L(TRANSPOSE) * INV D * INV L
40 DO 44 J=2,N
S=A(J-1,J-1)
DO 41 I=J,N
X=A(I,I)*A(I,J-1)
S=S+A(I,J-1)*X
41 A(I,J-1)=X
A(J-1,J-1)=S
IF(J.EQ.N) GO TO 50
DO 43 I=J,NM1
S=A(I,J-1)
DO 42 K=I,NM1
42 S=S+A(K+1,I)*A(K+1,J-1)
A(I,J-1)=S
43 CONTINUE
44 CONTINUE
C
50 IMX=ISAMAX(N,A(1,1),NDIA)
COND=1./ABS(A(IMX,IMX)*ZMX)
RETURN
C
91 COND=-FLOAT(J)
RETURN
END
INTEGER FUNCTION ISAMAX(N,A,M)
C
C FIND THE LARGEST ABSOLUTE ELEMENT OF A , SPACED M WORDS APART
C
DIMENSION A(*)
C
LARGE=1
IF(N.LE.1) GO TO 9
INDEX=1+M
DO 1 I=2,N
IF(ABS(A(INDEX)).GE.ABS(A(LARGE))) LARGE=I
1 INDEX=INDEX+M
9 ISAMAX=LARGE
RETURN
END
|
