home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS 1992 June / SIMTEL_0692.cdr / msdos / fortran / linpklib.arc / SPPDI.FOR < prev    next >
Encoding:
Text File  |  1984-01-07  |  3.4 KB  |  129 lines
  1.       SUBROUTINE SPPDI(AP,N,DET,JOB)
  2.       INTEGER N,JOB
  3.       REAL AP(1)
  4.       REAL DET(2)
  5. C
  6. C     SPPDI COMPUTES THE DETERMINANT AND INVERSE
  7. C     OF A REAL SYMMETRIC POSITIVE DEFINITE MATRIX
  8. C     USING THE FACTORS COMPUTED BY SPPCO OR SPPFA .
  9. C
  10. C     ON ENTRY
  11. C
  12. C        AP      REAL (N*(N+1)/2)
  13. C                THE OUTPUT FROM SPPCO OR SPPFA.
  14. C
  15. C        N       INTEGER
  16. C                THE ORDER OF THE MATRIX  A .
  17. C
  18. C        JOB     INTEGER
  19. C                = 11   BOTH DETERMINANT AND INVERSE.
  20. C                = 01   INVERSE ONLY.
  21. C                = 10   DETERMINANT ONLY.
  22. C
  23. C     ON RETURN
  24. C
  25. C        AP      THE UPPER TRIANGULAR HALF OF THE INVERSE .
  26. C                THE STRICT LOWER TRIANGLE IS UNALTERED.
  27. C
  28. C        DET     REAL(2)
  29. C                DETERMINANT OF ORIGINAL MATRIX IF REQUESTED.
  30. C                OTHERWISE NOT REFERENCED.
  31. C                DETERMINANT = DET(1) * 10.0**DET(2)
  32. C                WITH  1.0 .LE. DET(1) .LT. 10.0
  33. C                OR  DET(1) .EQ. 0.0 .
  34. C
  35. C     ERROR CONDITION
  36. C
  37. C        A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS
  38. C        A ZERO ON THE DIAGONAL AND THE INVERSE IS REQUESTED.
  39. C        IT WILL NOT OCCUR IF THE SUBROUTINES ARE CALLED CORRECTLY
  40. C        AND IF DPOCO OR DPOFA HAS SET INFO .EQ. 0 .
  41. C
  42. C     LINPACK.  THIS VERSION DATED 08/14/78 .
  43. C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
  44. C
  45. C     SUBROUTINES AND FUNCTIONS
  46. C
  47. C     BLAS SAXPY,SSCAL
  48. C     FORTRAN MOD
  49. C
  50. C     INTERNAL VARIABLES
  51. C
  52.       REAL T
  53.       REAL S
  54.       INTEGER I,II,J,JJ,JM1,J1,K,KJ,KK,KP1,K1
  55. C
  56. C     COMPUTE DETERMINANT
  57. C
  58.       IF (JOB/10 .EQ. 0) GO TO 70
  59.          DET(1) = 1.0E0
  60.          DET(2) = 0.0E0
  61.          S = 10.0E0
  62.          II = 0
  63.          DO 50 I = 1, N
  64.             II = II + I
  65.             DET(1) = AP(II)**2*DET(1)
  66. C        ...EXIT
  67.             IF (DET(1) .EQ. 0.0E0) GO TO 60
  68.    10       IF (DET(1) .GE. 1.0E0) GO TO 20
  69.                DET(1) = S*DET(1)
  70.                DET(2) = DET(2) - 1.0E0
  71.             GO TO 10
  72.    20       CONTINUE
  73.    30       IF (DET(1) .LT. S) GO TO 40
  74.                DET(1) = DET(1)/S
  75.                DET(2) = DET(2) + 1.0E0
  76.             GO TO 30
  77.    40       CONTINUE
  78.    50    CONTINUE
  79.    60    CONTINUE
  80.    70 CONTINUE
  81. C
  82. C     COMPUTE INVERSE(R)
  83. C
  84.       IF (MOD(JOB,10) .EQ. 0) GO TO 140
  85.          KK = 0
  86.          DO 100 K = 1, N
  87.             K1 = KK + 1
  88.             KK = KK + K
  89.             AP(KK) = 1.0E0/AP(KK)
  90.             T = -AP(KK)
  91.             CALL SSCAL(K-1,T,AP(K1),1)
  92.             KP1 = K + 1
  93.             J1 = KK + 1
  94.             KJ = KK + K
  95.             IF (N .LT. KP1) GO TO 90
  96.             DO 80 J = KP1, N
  97.                T = AP(KJ)
  98.                AP(KJ) = 0.0E0
  99.                CALL SAXPY(K,T,AP(K1),1,AP(J1),1)
  100.                J1 = J1 + J
  101.                KJ = KJ + J
  102.    80       CONTINUE
  103.    90       CONTINUE
  104.   100    CONTINUE
  105. C
  106. C        FORM  INVERSE(R) * TRANS(INVERSE(R))
  107. C
  108.          JJ = 0
  109.          DO 130 J = 1, N
  110.             J1 = JJ + 1
  111.             JJ = JJ + J
  112.             JM1 = J - 1
  113.             K1 = 1
  114.             KJ = J1
  115.             IF (JM1 .LT. 1) GO TO 120
  116.             DO 110 K = 1, JM1
  117.                T = AP(KJ)
  118.                CALL SAXPY(K,T,AP(J1),1,AP(K1),1)
  119.                K1 = K1 + K
  120.                KJ = KJ + 1
  121.   110       CONTINUE
  122.   120       CONTINUE
  123.             T = AP(JJ)
  124.             CALL SSCAL(J,T,AP(J1),1)
  125.   130    CONTINUE
  126.   140 CONTINUE
  127.       RETURN
  128.       END
  129.