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