|
molscat
|
|
README,
README.jkl,
README.v12,
dblas.f,
dblas.f.Z,
dcs_save.f,
diag_eispack.f,
ghm_save.f,
ghm_subs.f,
ghm_vib.f,
ident2disting.f,
lapack.f,
lapack.f.Z,
prbr_save.f,
prbr_vib.f,
read_isigu.f,
restrt.v12.f,
sbe.doc,
sbe_save.f,
sig_save.f,
spline.f,
syminv.f,
test1.input,
test1.v12.out,
test1.v14.out,
test2.input,
test2.v12.out,
test2.v14.out,
test3.f,
test3.input,
test3.v12.input,
test3.v12.out,
test3.v14.out,
test5.f,
test5.input,
test5.v12.out,
test5.v14.out,
test6.input,
test6.v12.out,
test6.v14.out,
test8.input,
test8.v12.out,
test8.v14.out,
timers.f,
timers_c.c,
v12.f,
v14.doc.tar,
v14.f,
v14.f.Z,
version_12.doc,
version_14.doc,
version_14.tutorial,
|
|
|
SUBROUTINE SYMINV (A, IA, N, IFAIL)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION A(IA,N)
C
C IN SITU INVERSION OF A REAL SYMMETRIC MATRIX BY L.D.L(T)
C DECOMPOSITION. KOUNT = NUMBER OF NEGATIVE EIGENVALUES.
C ROUTINE FAILS (IFAIL = N+1) IF THE DECOMPOSITION IS UNDEFINED.
C J.H.WILKINSON AND C.REINSCH, LINEAR ALGEBRA, I/1
C
C FOR VERSION 12 - NO LONGER NECESSARY TO RESTORE UPPER TRIANGLE
C
KOUNT = 0
X = A(1,1)
IF (X.LT.0.D0) KOUNT = KOUNT + 1
IF (X.EQ.0.D0) GO TO 320
A(1,1) = 1.D0/X
IF (N.EQ.1) GO TO 300
C FORM L
DO 100 I = 2,N
I1 = I - 1
IF (I1.LT.2) GO TO 60
DO 40 J = 2,I1
J1 = J - 1
X = A(I,J)
DO 20 K = 1,J1
X = X - A(I,K)*A(J,K)
20 CONTINUE
A(I,J) = X
40 CONTINUE
60 X = A(I,I)
DO 80 K = 1,I1
Y = A(I,K)
Z = Y*A(K,K)
A(I,K) = Z
X = X - Y*Z
80 CONTINUE
IF (X.LT.0.D0) KOUNT = KOUNT + 1
IF (X.EQ.0.D0) GO TO 320
A(I,I) = 1.D0/X
100 CONTINUE
C INVERT L
N1 = N - 1
DO 160 I = 1,N1
I1 = I + 1
A(I1,I) = - A(I1,I)
IF (I1.GT.N1) GO TO 160
DO 140 J = I1,N1
J1 = J + 1
X = - A(J1,I)
DO 120 K = I1,J
X = X - A(J1,K)*A(K,I)
120 CONTINUE
A(J1,I) = X
140 CONTINUE
160 CONTINUE
C FORM INVERSE OF A
DO 240 I = 1,N1
I1 = I + 1
X = A(I,I)
DO 180 K = I1,N
Y = A(K,I)
Z = Y*A(K,K)
A(K,I) = Z
X = X + Y*Z
180 CONTINUE
A(I,I) = X
IF (I1.GT.N1) GO TO 240
DO 220 J = I1,N1
J1 = J + 1
X = A(J,I)
DO 200 K = J1,N
X = X + A(K,J)*A(K,I)
200 CONTINUE
A(J,I) = X
220 CONTINUE
240 CONTINUE
C COPY INVERSE INTO UPPER TRIANGLE -----
C DO 280 I = 2,N -
C I1 = I - 1 - NO LONGER NEEDED
C DO 260 J = 1,I1 - IN VERSION 12
C A(J,I) = A(I,J) -
C260 CONTINUE -
C280 CONTINUE -----
300 IFAIL = KOUNT
RETURN
320 IFAIL = N + 1
RETURN
END
|
