home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS 1992 September / Simtel20_Sept92.cdr / msdos / fortran / linpklib.arc / STRDI.FOR < prev    next >
Text File  |  1984-01-12  |  4KB  |  139 lines
  1.       SUBROUTINE STRDI(T,LDT,N,DET,JOB,INFO)
  2.       INTEGER LDT,N,JOB,INFO
  3.       REAL T(LDT,1),DET(2)
  4. C
  5. C     STRDI COMPUTES THE DETERMINANT AND INVERSE OF A REAL
  6. C     TRIANGULAR MATRIX.
  7. C
  8. C     ON ENTRY
  9. C
  10. C        T       REAL(LDT,N)
  11. C                T CONTAINS THE TRIANGULAR MATRIX. THE ZERO
  12. C                ELEMENTS OF THE MATRIX ARE NOT REFERENCED, AND
  13. C                THE CORRESPONDING ELEMENTS OF THE ARRAY CAN BE
  14. C                USED TO STORE OTHER INFORMATION.
  15. C
  16. C        LDT     INTEGER
  17. C                LDT IS THE LEADING DIMENSION OF THE ARRAY T.
  18. C
  19. C        N       INTEGER
  20. C                N IS THE ORDER OF THE SYSTEM.
  21. C
  22. C        JOB     INTEGER
  23. C                = 010       NO DET, INVERSE OF LOWER TRIANGULAR.
  24. C                = 011       NO DET, INVERSE OF UPPER TRIANGULAR.
  25. C                = 100       DET, NO INVERSE.
  26. C                = 110       DET, INVERSE OF LOWER TRIANGULAR.
  27. C                = 111       DET, INVERSE OF UPPER TRIANGULAR.
  28. C
  29. C     ON RETURN
  30. C
  31. C        T       INVERSE OF ORIGINAL MATRIX IF REQUESTED.
  32. C                OTHERWISE UNCHANGED.
  33. C
  34. C        DET     REAL(2)
  35. C                DETERMINANT OF ORIGINAL MATRIX IF REQUESTED.
  36. C                OTHERWISE NOT REFERENCED.
  37. C                DETERMINANT = DET(1) * 10.0**DET(2)
  38. C                WITH  1.0 .LE. ABS(DET(1)) .LT. 10.0
  39. C                OR  DET(1) .EQ. 0.0 .
  40. C
  41. C        INFO    INTEGER
  42. C                INFO CONTAINS ZERO IF THE SYSTEM IS NONSINGULAR
  43. C                AND THE INVERSE IS REQUESTED.
  44. C                OTHERWISE INFO CONTAINS THE INDEX OF
  45. C                A ZERO DIAGONAL ELEMENT OF T.
  46. C
  47. C
  48. C     LINPACK. THIS VERSION DATED 08/14/78 .
  49. C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
  50. C
  51. C     SUBROUTINES AND FUNCTIONS
  52. C
  53. C     BLAS SAXPY,SSCAL
  54. C     FORTRAN ABS,MOD
  55. C
  56. C     INTERNAL VARIABLES
  57. C
  58.       REAL TEMP
  59.       REAL TEN
  60.       INTEGER I,J,K,KB,KM1,KP1
  61. C
  62. C     BEGIN BLOCK PERMITTING ...EXITS TO 180
  63. C
  64. C        COMPUTE DETERMINANT
  65. C
  66.          IF (JOB/100 .EQ. 0) GO TO 70
  67.             DET(1) = 1.0E0
  68.             DET(2) = 0.0E0
  69.             TEN = 10.0E0
  70.             DO 50 I = 1, N
  71.                DET(1) = T(I,I)*DET(1)
  72. C           ...EXIT
  73.                IF (DET(1) .EQ. 0.0E0) GO TO 60
  74.    10          IF (ABS(DET(1)) .GE. 1.0E0) GO TO 20
  75.                   DET(1) = TEN*DET(1)
  76.                   DET(2) = DET(2) - 1.0E0
  77.                GO TO 10
  78.    20          CONTINUE
  79.    30          IF (ABS(DET(1)) .LT. TEN) GO TO 40
  80.                   DET(1) = DET(1)/TEN
  81.                   DET(2) = DET(2) + 1.0E0
  82.                GO TO 30
  83.    40          CONTINUE
  84.    50       CONTINUE
  85.    60       CONTINUE
  86.    70    CONTINUE
  87. C
  88. C        COMPUTE INVERSE OF UPPER TRIANGULAR
  89. C
  90.          IF (MOD(JOB/10,10) .EQ. 0) GO TO 170
  91.             IF (MOD(JOB,10) .EQ. 0) GO TO 120
  92. C              BEGIN BLOCK PERMITTING ...EXITS TO 110
  93.                   DO 100 K = 1, N
  94.                      INFO = K
  95. C              ......EXIT
  96.                      IF (T(K,K) .EQ. 0.0E0) GO TO 110
  97.                      T(K,K) = 1.0E0/T(K,K)
  98.                      TEMP = -T(K,K)
  99.                      CALL SSCAL(K-1,TEMP,T(1,K),1)
  100.                      KP1 = K + 1
  101.                      IF (N .LT. KP1) GO TO 90
  102.                      DO 80 J = KP1, N
  103.                         TEMP = T(K,J)
  104.                         T(K,J) = 0.0E0
  105.                         CALL SAXPY(K,TEMP,T(1,K),1,T(1,J),1)
  106.    80                CONTINUE
  107.    90                CONTINUE
  108.   100             CONTINUE
  109.                   INFO = 0
  110.   110          CONTINUE
  111.             GO TO 160
  112.   120       CONTINUE
  113. C
  114. C              COMPUTE INVERSE OF LOWER TRIANGULAR
  115. C
  116.                DO 150 KB = 1, N
  117.                   K = N + 1 - KB
  118.                   INFO = K
  119. C     ............EXIT
  120.                   IF (T(K,K) .EQ. 0.0E0) GO TO 180
  121.                   T(K,K) = 1.0E0/T(K,K)
  122.                   TEMP = -T(K,K)
  123.                   IF (K .NE. N) CALL SSCAL(N-K,TEMP,T(K+1,K),1)
  124.                   KM1 = K - 1
  125.                   IF (KM1 .LT. 1) GO TO 140
  126.                   DO 130 J = 1, KM1
  127.                      TEMP = T(K,J)
  128.                      T(K,J) = 0.0E0
  129.                      CALL SAXPY(N-K+1,TEMP,T(K,K),1,T(K,J),1)
  130.   130             CONTINUE
  131.   140             CONTINUE
  132.   150          CONTINUE
  133.                INFO = 0
  134.   160       CONTINUE
  135.   170    CONTINUE
  136.   180 CONTINUE
  137.       RETURN
  138.       END
  139.