home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS 1992 June / SIMTEL_0692.cdr / msdos / fortran / linpklib.arc / SSIDI.FOR < prev    next >
Encoding:
Text File  |  1984-01-08  |  6.5 KB  |  213 lines
  1.       SUBROUTINE SSIDI(A,LDA,N,KPVT,DET,INERT,WORK,JOB)
  2.       INTEGER LDA,N,JOB
  3.       REAL A(LDA,1),WORK(1)
  4.       REAL DET(2)
  5.       INTEGER KPVT(1),INERT(3)
  6. C
  7. C     SSIDI COMPUTES THE DETERMINANT, INERTIA AND INVERSE
  8. C     OF A REAL SYMMETRIC MATRIX USING THE FACTORS FROM SSIFA.
  9. C
  10. C     ON ENTRY
  11. C
  12. C        A       REAL(LDA,N)
  13. C                THE OUTPUT FROM SSIFA.
  14. C
  15. C        LDA     INTEGER
  16. C                THE LEADING DIMENSION OF THE ARRAY A.
  17. C
  18. C        N       INTEGER
  19. C                THE ORDER OF THE MATRIX A.
  20. C
  21. C        KPVT    INTEGER(N)
  22. C                THE PIVOT VECTOR FROM SSIFA.
  23. C
  24. C        WORK    REAL(N)
  25. C                WORK VECTOR.  CONTENTS DESTROYED.
  26. C
  27. C        JOB     INTEGER
  28. C                JOB HAS THE DECIMAL EXPANSION  ABC  WHERE
  29. C                   IF  C .NE. 0, THE INVERSE IS COMPUTED,
  30. C                   IF  B .NE. 0, THE DETERMINANT IS COMPUTED,
  31. C                   IF  A .NE. 0, THE INERTIA IS COMPUTED.
  32. C
  33. C                FOR EXAMPLE, JOB = 111  GIVES ALL THREE.
  34. C
  35. C     ON RETURN
  36. C
  37. C        VARIABLES NOT REQUESTED BY JOB ARE NOT USED.
  38. C
  39. C        A      CONTAINS THE UPPER TRIANGLE OF THE INVERSE OF
  40. C               THE ORIGINAL MATRIX.  THE STRICT LOWER TRIANGLE
  41. C               IS NEVER REFERENCED.
  42. C
  43. C        DET    REAL(2)
  44. C               DETERMINANT OF ORIGINAL MATRIX.
  45. C               DETERMINANT = DET(1) * 10.0**DET(2)
  46. C               WITH 1.0 .LE. ABS(DET(1)) .LT. 10.0
  47. C               OR DET(1) = 0.0.
  48. C
  49. C        INERT  INTEGER(3)
  50. C               THE INERTIA OF THE ORIGINAL MATRIX.
  51. C               INERT(1)  =  NUMBER OF POSITIVE EIGENVALUES.
  52. C               INERT(2)  =  NUMBER OF NEGATIVE EIGENVALUES.
  53. C               INERT(3)  =  NUMBER OF ZERO EIGENVALUES.
  54. C
  55. C     ERROR CONDITION
  56. C
  57. C        A DIVISION BY ZERO MAY OCCUR IF THE INVERSE IS REQUESTED
  58. C        AND  SSICO  HAS SET RCOND .EQ. 0.0
  59. C        OR  SSIFA  HAS SET  INFO .NE. 0 .
  60. C
  61. C     LINPACK. THIS VERSION DATED 08/14/78 .
  62. C     JAMES BUNCH, UNIV. CALIF. SAN DIEGO, ARGONNE NAT. LAB
  63. C
  64. C     SUBROUTINES AND FUNCTIONS
  65. C
  66. C     BLAS SAXPY,SCOPY,SDOT,SSWAP
  67. C     FORTRAN ABS,IABS,MOD
  68. C
  69. C     INTERNAL VARIABLES.
  70. C
  71.       REAL AKKP1,SDOT,TEMP
  72.       REAL TEN,D,T,AK,AKP1
  73.       INTEGER J,JB,K,KM1,KS,KSTEP
  74.       LOGICAL NOINV,NODET,NOERT
  75. C
  76.       NOINV = MOD(JOB,10) .EQ. 0
  77.       NODET = MOD(JOB,100)/10 .EQ. 0
  78.       NOERT = MOD(JOB,1000)/100 .EQ. 0
  79. C
  80.       IF (NODET .AND. NOERT) GO TO 140
  81.          IF (NOERT) GO TO 10
  82.             INERT(1) = 0
  83.             INERT(2) = 0
  84.             INERT(3) = 0
  85.    10    CONTINUE
  86.          IF (NODET) GO TO 20
  87.             DET(1) = 1.0E0
  88.             DET(2) = 0.0E0
  89.             TEN = 10.0E0
  90.    20    CONTINUE
  91.          T = 0.0E0
  92.          DO 130 K = 1, N
  93.             D = A(K,K)
  94. C
  95. C           CHECK IF 1 BY 1
  96. C
  97.             IF (KPVT(K) .GT. 0) GO TO 50
  98. C
  99. C              2 BY 2 BLOCK
  100. C              USE DET (D  S)  =  (D/T * C - T) * T  ,  T = ABS(S)
  101. C                      (S  C)
  102. C              TO AVOID UNDERFLOW/OVERFLOW TROUBLES.
  103. C              TAKE TWO PASSES THROUGH SCALING.  USE  T  FOR FLAG.
  104. C
  105.                IF (T .NE. 0.0E0) GO TO 30
  106.                   T = ABS(A(K,K+1))
  107.                   D = (D/T)*A(K+1,K+1) - T
  108.                GO TO 40
  109.    30          CONTINUE
  110.                   D = T
  111.                   T = 0.0E0
  112.    40          CONTINUE
  113.    50       CONTINUE
  114. C
  115.             IF (NOERT) GO TO 60
  116.                IF (D .GT. 0.0E0) INERT(1) = INERT(1) + 1
  117.                IF (D .LT. 0.0E0) INERT(2) = INERT(2) + 1
  118.                IF (D .EQ. 0.0E0) INERT(3) = INERT(3) + 1
  119.    60       CONTINUE
  120. C
  121.             IF (NODET) GO TO 120
  122.                DET(1) = D*DET(1)
  123.                IF (DET(1) .EQ. 0.0E0) GO TO 110
  124.    70             IF (ABS(DET(1)) .GE. 1.0E0) GO TO 80
  125.                      DET(1) = TEN*DET(1)
  126.                      DET(2) = DET(2) - 1.0E0
  127.                   GO TO 70
  128.    80             CONTINUE
  129.    90             IF (ABS(DET(1)) .LT. TEN) GO TO 100
  130.                      DET(1) = DET(1)/TEN
  131.                      DET(2) = DET(2) + 1.0E0
  132.                   GO TO 90
  133.   100             CONTINUE
  134.   110          CONTINUE
  135.   120       CONTINUE
  136.   130    CONTINUE
  137.   140 CONTINUE
  138. C
  139. C     COMPUTE INVERSE(A)
  140. C
  141.       IF (NOINV) GO TO 270
  142.          K = 1
  143.   150    IF (K .GT. N) GO TO 260
  144.             KM1 = K - 1
  145.             IF (KPVT(K) .LT. 0) GO TO 180
  146. C
  147. C              1 BY 1
  148. C
  149.                A(K,K) = 1.0E0/A(K,K)
  150.                IF (KM1 .LT. 1) GO TO 170
  151.                   CALL SCOPY(KM1,A(1,K),1,WORK,1)
  152.                   DO 160 J = 1, KM1
  153.                      A(J,K) = SDOT(J,A(1,J),1,WORK,1)
  154.                      CALL SAXPY(J-1,WORK(J),A(1,J),1,A(1,K),1)
  155.   160             CONTINUE
  156.                   A(K,K) = A(K,K) + SDOT(KM1,WORK,1,A(1,K),1)
  157.   170          CONTINUE
  158.                KSTEP = 1
  159.             GO TO 220
  160.   180       CONTINUE
  161. C
  162. C              2 BY 2
  163. C
  164.                T = ABS(A(K,K+1))
  165.                AK = A(K,K)/T
  166.                AKP1 = A(K+1,K+1)/T
  167.                AKKP1 = A(K,K+1)/T
  168.                D = T*(AK*AKP1 - 1.0E0)
  169.                A(K,K) = AKP1/D
  170.                A(K+1,K+1) = AK/D
  171.                A(K,K+1) = -AKKP1/D
  172.                IF (KM1 .LT. 1) GO TO 210
  173.                   CALL SCOPY(KM1,A(1,K+1),1,WORK,1)
  174.                   DO 190 J = 1, KM1
  175.                      A(J,K+1) = SDOT(J,A(1,J),1,WORK,1)
  176.                      CALL SAXPY(J-1,WORK(J),A(1,J),1,A(1,K+1),1)
  177.   190             CONTINUE
  178.                   A(K+1,K+1) = A(K+1,K+1) + SDOT(KM1,WORK,1,A(1,K+1),1)
  179.                   A(K,K+1) = A(K,K+1) + SDOT(KM1,A(1,K),1,A(1,K+1),1)
  180.                   CALL SCOPY(KM1,A(1,K),1,WORK,1)
  181.                   DO 200 J = 1, KM1
  182.                      A(J,K) = SDOT(J,A(1,J),1,WORK,1)
  183.                      CALL SAXPY(J-1,WORK(J),A(1,J),1,A(1,K),1)
  184.   200             CONTINUE
  185.                   A(K,K) = A(K,K) + SDOT(KM1,WORK,1,A(1,K),1)
  186.   210          CONTINUE
  187.                KSTEP = 2
  188.   220       CONTINUE
  189. C
  190. C           SWAP
  191. C
  192.             KS = IABS(KPVT(K))
  193.             IF (KS .EQ. K) GO TO 250
  194.                CALL SSWAP(KS,A(1,KS),1,A(1,K),1)
  195.                DO 230 JB = KS, K
  196.                   J = K + KS - JB
  197.                   TEMP = A(J,K)
  198.                   A(J,K) = A(KS,J)
  199.                   A(KS,J) = TEMP
  200.   230          CONTINUE
  201.                IF (KSTEP .EQ. 1) GO TO 240
  202.                   TEMP = A(KS,K+1)
  203.                   A(KS,K+1) = A(K,K+1)
  204.                   A(K,K+1) = TEMP
  205.   240          CONTINUE
  206.   250       CONTINUE
  207.             K = K + KSTEP
  208.          GO TO 150
  209.   260    CONTINUE
  210.   270 CONTINUE
  211.       RETURN
  212.       END
  213.