home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS 1992 September / Simtel20_Sept92.cdr / msdos / fortran / linpklib.arc / STRCO.FOR < prev    next >
Text File  |  1984-01-12  |  5KB  |  161 lines
  1.       SUBROUTINE STRCO(T,LDT,N,RCOND,Z,JOB)
  2.       INTEGER LDT,N,JOB
  3.       REAL T(LDT,1),Z(1)
  4.       REAL RCOND
  5. C
  6. C     STRCO ESTIMATES THE CONDITION OF A REAL 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                = 0         T  IS LOWER TRIANGULAR.
  24. C                = NONZERO   T  IS UPPER TRIANGULAR.
  25. C
  26. C     ON RETURN
  27. C
  28. C        RCOND   REAL
  29. C                AN ESTIMATE OF THE RECIPROCAL CONDITION OF  T .
  30. C                FOR THE SYSTEM  T*X = B , RELATIVE PERTURBATIONS
  31. C                IN  T  AND  B  OF SIZE  EPSILON  MAY CAUSE
  32. C                RELATIVE PERTURBATIONS IN  X  OF SIZE  EPSILON/RCOND .
  33. C                IF  RCOND  IS SO SMALL THAT THE LOGICAL EXPRESSION
  34. C                           1.0 + RCOND .EQ. 1.0
  35. C                IS TRUE, THEN  T  MAY BE SINGULAR TO WORKING
  36. C                PRECISION.  IN PARTICULAR,  RCOND  IS ZERO  IF
  37. C                EXACT SINGULARITY IS DETECTED OR THE ESTIMATE
  38. C                UNDERFLOWS.
  39. C
  40. C        Z       REAL(N)
  41. C                A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT.
  42. C                IF  T  IS CLOSE TO A SINGULAR MATRIX, THEN  Z  IS
  43. C                AN APPROXIMATE NULL VECTOR IN THE SENSE THAT
  44. C                NORM(A*Z) = RCOND*NORM(A)*NORM(Z) .
  45. C
  46. C     LINPACK. THIS VERSION DATED 08/14/78 .
  47. C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
  48. C
  49. C     SUBROUTINES AND FUNCTIONS
  50. C
  51. C     BLAS SAXPY,SSCAL,SASUM
  52. C     FORTRAN ABS,AMAX1,SIGN
  53. C
  54. C     INTERNAL VARIABLES
  55. C
  56.       REAL W,WK,WKM,EK
  57.       REAL TNORM,YNORM,S,SM,SASUM
  58.       INTEGER I1,J,J1,J2,K,KK,L
  59.       LOGICAL LOWER
  60. C
  61.       LOWER = JOB .EQ. 0
  62. C
  63. C     COMPUTE 1-NORM OF T
  64. C
  65.       TNORM = 0.0E0
  66.       DO 10 J = 1, N
  67.          L = J
  68.          IF (LOWER) L = N + 1 - J
  69.          I1 = 1
  70.          IF (LOWER) I1 = J
  71.          TNORM = AMAX1(TNORM,SASUM(L,T(I1,J),1))
  72.    10 CONTINUE
  73. C
  74. C     RCOND = 1/(NORM(T)*(ESTIMATE OF NORM(INVERSE(T)))) .
  75. C     ESTIMATE = NORM(Z)/NORM(Y) WHERE  T*Z = Y  AND  TRANS(T)*Y = E .
  76. C     TRANS(T)  IS THE TRANSPOSE OF T .
  77. C     THE COMPONENTS OF  E  ARE CHOSEN TO CAUSE MAXIMUM LOCAL
  78. C     GROWTH IN THE ELEMENTS OF Y .
  79. C     THE VECTORS ARE FREQUENTLY RESCALED TO AVOID OVERFLOW.
  80. C
  81. C     SOLVE TRANS(T)*Y = E
  82. C
  83.       EK = 1.0E0
  84.       DO 20 J = 1, N
  85.          Z(J) = 0.0E0
  86.    20 CONTINUE
  87.       DO 100 KK = 1, N
  88.          K = KK
  89.          IF (LOWER) K = N + 1 - KK
  90.          IF (Z(K) .NE. 0.0E0) EK = SIGN(EK,-Z(K))
  91.          IF (ABS(EK-Z(K)) .LE. ABS(T(K,K))) GO TO 30
  92.             S = ABS(T(K,K))/ABS(EK-Z(K))
  93.             CALL SSCAL(N,S,Z,1)
  94.             EK = S*EK
  95.    30    CONTINUE
  96.          WK = EK - Z(K)
  97.          WKM = -EK - Z(K)
  98.          S = ABS(WK)
  99.          SM = ABS(WKM)
  100.          IF (T(K,K) .EQ. 0.0E0) GO TO 40
  101.             WK = WK/T(K,K)
  102.             WKM = WKM/T(K,K)
  103.          GO TO 50
  104.    40    CONTINUE
  105.             WK = 1.0E0
  106.             WKM = 1.0E0
  107.    50    CONTINUE
  108.          IF (KK .EQ. N) GO TO 90
  109.             J1 = K + 1
  110.             IF (LOWER) J1 = 1
  111.             J2 = N
  112.             IF (LOWER) J2 = K - 1
  113.             DO 60 J = J1, J2
  114.                SM = SM + ABS(Z(J)+WKM*T(K,J))
  115.                Z(J) = Z(J) + WK*T(K,J)
  116.                S = S + ABS(Z(J))
  117.    60       CONTINUE
  118.             IF (S .GE. SM) GO TO 80
  119.                W = WKM - WK
  120.                WK = WKM
  121.                DO 70 J = J1, J2
  122.                   Z(J) = Z(J) + W*T(K,J)
  123.    70          CONTINUE
  124.    80       CONTINUE
  125.    90    CONTINUE
  126.          Z(K) = WK
  127.   100 CONTINUE
  128.       S = 1.0E0/SASUM(N,Z,1)
  129.       CALL SSCAL(N,S,Z,1)
  130. C
  131.       YNORM = 1.0E0
  132. C
  133. C     SOLVE T*Z = Y
  134. C
  135.       DO 130 KK = 1, N
  136.          K = N + 1 - KK
  137.          IF (LOWER) K = KK
  138.          IF (ABS(Z(K)) .LE. ABS(T(K,K))) GO TO 110
  139.             S = ABS(T(K,K))/ABS(Z(K))
  140.             CALL SSCAL(N,S,Z,1)
  141.             YNORM = S*YNORM
  142.   110    CONTINUE
  143.          IF (T(K,K) .NE. 0.0E0) Z(K) = Z(K)/T(K,K)
  144.          IF (T(K,K) .EQ. 0.0E0) Z(K) = 1.0E0
  145.          I1 = 1
  146.          IF (LOWER) I1 = K + 1
  147.          IF (KK .GE. N) GO TO 120
  148.             W = -Z(K)
  149.             CALL SAXPY(N-KK,W,T(I1,K),1,Z(I1),1)
  150.   120    CONTINUE
  151.   130 CONTINUE
  152. C     MAKE ZNORM = 1.0
  153.       S = 1.0E0/SASUM(N,Z,1)
  154.       CALL SSCAL(N,S,Z,1)
  155.       YNORM = S*YNORM
  156. C
  157.       IF (TNORM .NE. 0.0E0) RCOND = YNORM/TNORM
  158.       IF (TNORM .EQ. 0.0E0) RCOND = 0.0E0
  159.       RETURN
  160.       END
  161.