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