CCL Home Page
Up Directory CCL syminv
      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
Modified: Wed Mar 8 17:00:00 1995 GMT
Page accessed 5314 times since Sat Apr 17 21:25:26 1999 GMT