home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / libcruft / linpack / dgeco.f < prev    next >
Text File  |  1996-09-28  |  6KB  |  194 lines
  1.       SUBROUTINE DGECO(A,LDA,N,IPVT,RCOND,Z)
  2.       INTEGER LDA,N,IPVT(1)
  3.       DOUBLE PRECISION A(LDA,1),Z(1)
  4.       DOUBLE PRECISION RCOND
  5. C
  6. C     DGECO FACTORS A DOUBLE PRECISION MATRIX BY GAUSSIAN ELIMINATION
  7. C     AND ESTIMATES THE CONDITION OF THE MATRIX.
  8. C
  9. C     IF  RCOND  IS NOT NEEDED, DGEFA IS SLIGHTLY FASTER.
  10. C     TO SOLVE  A*X = B , FOLLOW DGECO BY DGESL.
  11. C     TO COMPUTE  INVERSE(A)*C , FOLLOW DGECO BY DGESL.
  12. C     TO COMPUTE  DETERMINANT(A) , FOLLOW DGECO BY DGEDI.
  13. C     TO COMPUTE  INVERSE(A) , FOLLOW DGECO BY DGEDI.
  14. C
  15. C     ON ENTRY
  16. C
  17. C        A       DOUBLE PRECISION(LDA, N)
  18. C                THE MATRIX TO BE FACTORED.
  19. C
  20. C        LDA     INTEGER
  21. C                THE LEADING DIMENSION OF THE ARRAY  A .
  22. C
  23. C        N       INTEGER
  24. C                THE ORDER OF THE MATRIX  A .
  25. C
  26. C     ON RETURN
  27. C
  28. C        A       AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS
  29. C                WHICH WERE USED TO OBTAIN IT.
  30. C                THE FACTORIZATION CAN BE WRITTEN  A = L*U  WHERE
  31. C                L  IS A PRODUCT OF PERMUTATION AND UNIT LOWER
  32. C                TRIANGULAR MATRICES AND  U  IS UPPER TRIANGULAR.
  33. C
  34. C        IPVT    INTEGER(N)
  35. C                AN INTEGER VECTOR OF PIVOT INDICES.
  36. C
  37. C        RCOND   DOUBLE PRECISION
  38. C                AN ESTIMATE OF THE RECIPROCAL CONDITION OF  A .
  39. C                FOR THE SYSTEM  A*X = B , RELATIVE PERTURBATIONS
  40. C                IN  A  AND  B  OF SIZE  EPSILON  MAY CAUSE
  41. C                RELATIVE PERTURBATIONS IN  X  OF SIZE  EPSILON/RCOND .
  42. C                IF  RCOND  IS SO SMALL THAT THE LOGICAL EXPRESSION
  43. C                           1.0 + RCOND .EQ. 1.0
  44. C                IS TRUE, THEN  A  MAY BE SINGULAR TO WORKING
  45. C                PRECISION.  IN PARTICULAR,  RCOND  IS ZERO  IF
  46. C                EXACT SINGULARITY IS DETECTED OR THE ESTIMATE
  47. C                UNDERFLOWS.
  48. C
  49. C        Z       DOUBLE PRECISION(N)
  50. C                A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT.
  51. C                IF  A  IS CLOSE TO A SINGULAR MATRIX, THEN  Z  IS
  52. C                AN APPROXIMATE NULL VECTOR IN THE SENSE THAT
  53. C                NORM(A*Z) = RCOND*NORM(A)*NORM(Z) .
  54. C
  55. C     LINPACK. THIS VERSION DATED 08/14/78 .
  56. C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
  57. C
  58. C     SUBROUTINES AND FUNCTIONS
  59. C
  60. C     LINPACK DGEFA
  61. C     BLAS DAXPY,DDOT,DSCAL,DASUM
  62. C     FORTRAN DABS,DMAX1,DSIGN
  63. C
  64. C     INTERNAL VARIABLES
  65. C
  66.       DOUBLE PRECISION DDOT,EK,T,WK,WKM
  67.       DOUBLE PRECISION ANORM,S,DASUM,SM,YNORM
  68.       INTEGER INFO,J,K,KB,KP1,L
  69. C
  70. C
  71. C     COMPUTE 1-NORM OF A
  72. C
  73.       ANORM = 0.0D0
  74.       DO 10 J = 1, N
  75.          ANORM = DMAX1(ANORM,DASUM(N,A(1,J),1))
  76.    10 CONTINUE
  77. C
  78. C     FACTOR
  79. C
  80.       CALL DGEFA(A,LDA,N,IPVT,INFO)
  81. C
  82. C     RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) .
  83. C     ESTIMATE = NORM(Z)/NORM(Y) WHERE  A*Z = Y  AND  TRANS(A)*Y = E .
  84. C     TRANS(A)  IS THE TRANSPOSE OF A .  THE COMPONENTS OF  E  ARE
  85. C     CHOSEN TO CAUSE MAXIMUM LOCAL GROWTH IN THE ELEMENTS OF W  WHERE
  86. C     TRANS(U)*W = E .  THE VECTORS ARE FREQUENTLY RESCALED TO AVOID
  87. C     OVERFLOW.
  88. C
  89. C     SOLVE TRANS(U)*W = E
  90. C
  91.       EK = 1.0D0
  92.       DO 20 J = 1, N
  93.          Z(J) = 0.0D0
  94.    20 CONTINUE
  95.       DO 100 K = 1, N
  96.          IF (Z(K) .NE. 0.0D0) EK = DSIGN(EK,-Z(K))
  97.          IF (DABS(EK-Z(K)) .LE. DABS(A(K,K))) GO TO 30
  98.             S = DABS(A(K,K))/DABS(EK-Z(K))
  99.             CALL DSCAL(N,S,Z,1)
  100.             EK = S*EK
  101.    30    CONTINUE
  102.          WK = EK - Z(K)
  103.          WKM = -EK - Z(K)
  104.          S = DABS(WK)
  105.          SM = DABS(WKM)
  106.          IF (A(K,K) .EQ. 0.0D0) GO TO 40
  107.             WK = WK/A(K,K)
  108.             WKM = WKM/A(K,K)
  109.          GO TO 50
  110.    40    CONTINUE
  111.             WK = 1.0D0
  112.             WKM = 1.0D0
  113.    50    CONTINUE
  114.          KP1 = K + 1
  115.          IF (KP1 .GT. N) GO TO 90
  116.             DO 60 J = KP1, N
  117.                SM = SM + DABS(Z(J)+WKM*A(K,J))
  118.                Z(J) = Z(J) + WK*A(K,J)
  119.                S = S + DABS(Z(J))
  120.    60       CONTINUE
  121.             IF (S .GE. SM) GO TO 80
  122.                T = WKM - WK
  123.                WK = WKM
  124.                DO 70 J = KP1, N
  125.                   Z(J) = Z(J) + T*A(K,J)
  126.    70          CONTINUE
  127.    80       CONTINUE
  128.    90    CONTINUE
  129.          Z(K) = WK
  130.   100 CONTINUE
  131.       S = 1.0D0/DASUM(N,Z,1)
  132.       CALL DSCAL(N,S,Z,1)
  133. C
  134. C     SOLVE TRANS(L)*Y = W
  135. C
  136.       DO 120 KB = 1, N
  137.          K = N + 1 - KB
  138.          IF (K .LT. N) Z(K) = Z(K) + DDOT(N-K,A(K+1,K),1,Z(K+1),1)
  139.          IF (DABS(Z(K)) .LE. 1.0D0) GO TO 110
  140.             S = 1.0D0/DABS(Z(K))
  141.             CALL DSCAL(N,S,Z,1)
  142.   110    CONTINUE
  143.          L = IPVT(K)
  144.          T = Z(L)
  145.          Z(L) = Z(K)
  146.          Z(K) = T
  147.   120 CONTINUE
  148.       S = 1.0D0/DASUM(N,Z,1)
  149.       CALL DSCAL(N,S,Z,1)
  150. C
  151.       YNORM = 1.0D0
  152. C
  153. C     SOLVE L*V = Y
  154. C
  155.       DO 140 K = 1, N
  156.          L = IPVT(K)
  157.          T = Z(L)
  158.          Z(L) = Z(K)
  159.          Z(K) = T
  160.          IF (K .LT. N) CALL DAXPY(N-K,T,A(K+1,K),1,Z(K+1),1)
  161.          IF (DABS(Z(K)) .LE. 1.0D0) GO TO 130
  162.             S = 1.0D0/DABS(Z(K))
  163.             CALL DSCAL(N,S,Z,1)
  164.             YNORM = S*YNORM
  165.   130    CONTINUE
  166.   140 CONTINUE
  167.       S = 1.0D0/DASUM(N,Z,1)
  168.       CALL DSCAL(N,S,Z,1)
  169.       YNORM = S*YNORM
  170. C
  171. C     SOLVE  U*Z = V
  172. C
  173.       DO 160 KB = 1, N
  174.          K = N + 1 - KB
  175.          IF (DABS(Z(K)) .LE. DABS(A(K,K))) GO TO 150
  176.             S = DABS(A(K,K))/DABS(Z(K))
  177.             CALL DSCAL(N,S,Z,1)
  178.             YNORM = S*YNORM
  179.   150    CONTINUE
  180.          IF (A(K,K) .NE. 0.0D0) Z(K) = Z(K)/A(K,K)
  181.          IF (A(K,K) .EQ. 0.0D0) Z(K) = 1.0D0
  182.          T = -Z(K)
  183.          CALL DAXPY(K-1,T,A(1,K),1,Z(1),1)
  184.   160 CONTINUE
  185. C     MAKE ZNORM = 1.0
  186.       S = 1.0D0/DASUM(N,Z,1)
  187.       CALL DSCAL(N,S,Z,1)
  188.       YNORM = S*YNORM
  189. C
  190.       IF (ANORM .NE. 0.0D0) RCOND = YNORM/ANORM
  191.       IF (ANORM .EQ. 0.0D0) RCOND = 0.0D0
  192.       RETURN
  193.       END
  194.