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 / dgefa.f < prev    next >
Text File  |  1996-09-28  |  3KB  |  104 lines
  1.       SUBROUTINE DGEFA(A,LDA,N,IPVT,INFO)
  2.       INTEGER LDA,N,IPVT(1),INFO
  3.       DOUBLE PRECISION A(LDA,1)
  4. C
  5. C     DGEFA FACTORS A DOUBLE PRECISION MATRIX BY GAUSSIAN ELIMINATION.
  6. C
  7. C     DGEFA IS USUALLY CALLED BY DGECO, BUT IT CAN BE CALLED
  8. C     DIRECTLY WITH A SAVING IN TIME IF  RCOND  IS NOT NEEDED.
  9. C     (TIME FOR DGECO) = (1 + 9/N)*(TIME FOR DGEFA) .
  10. C
  11. C     ON ENTRY
  12. C
  13. C        A       DOUBLE PRECISION(LDA, N)
  14. C                THE MATRIX TO BE FACTORED.
  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     ON RETURN
  23. C
  24. C        A       AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS
  25. C                WHICH WERE USED TO OBTAIN IT.
  26. C                THE FACTORIZATION CAN BE WRITTEN  A = L*U  WHERE
  27. C                L  IS A PRODUCT OF PERMUTATION AND UNIT LOWER
  28. C                TRIANGULAR MATRICES AND  U  IS UPPER TRIANGULAR.
  29. C
  30. C        IPVT    INTEGER(N)
  31. C                AN INTEGER VECTOR OF PIVOT INDICES.
  32. C
  33. C        INFO    INTEGER
  34. C                = 0  NORMAL VALUE.
  35. C                = K  IF  U(K,K) .EQ. 0.0 .  THIS IS NOT AN ERROR
  36. C                     CONDITION FOR THIS SUBROUTINE, BUT IT DOES
  37. C                     INDICATE THAT DGESL OR DGEDI WILL DIVIDE BY ZERO
  38. C                     IF CALLED.  USE  RCOND  IN DGECO FOR A RELIABLE
  39. C                     INDICATION OF SINGULARITY.
  40. C
  41. C     LINPACK. THIS VERSION DATED 08/14/78 .
  42. C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
  43. C
  44. C     SUBROUTINES AND FUNCTIONS
  45. C
  46. C     BLAS DAXPY,DSCAL,IDAMAX
  47. C
  48. C     INTERNAL VARIABLES
  49. C
  50.       DOUBLE PRECISION T
  51.       INTEGER IDAMAX,J,K,KP1,L,NM1
  52. C
  53. C
  54. C     GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
  55. C
  56.       INFO = 0
  57.       NM1 = N - 1
  58.       IF (NM1 .LT. 1) GO TO 70
  59.       DO 60 K = 1, NM1
  60.          KP1 = K + 1
  61. C
  62. C        FIND L = PIVOT INDEX
  63. C
  64.          L = IDAMAX(N-K+1,A(K,K),1) + K - 1
  65.          IPVT(K) = L
  66. C
  67. C        ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
  68. C
  69.          IF (A(L,K) .EQ. 0.0D0) GO TO 40
  70. C
  71. C           INTERCHANGE IF NECESSARY
  72. C
  73.             IF (L .EQ. K) GO TO 10
  74.                T = A(L,K)
  75.                A(L,K) = A(K,K)
  76.                A(K,K) = T
  77.    10       CONTINUE
  78. C
  79. C           COMPUTE MULTIPLIERS
  80. C
  81.             T = -1.0D0/A(K,K)
  82.             CALL DSCAL(N-K,T,A(K+1,K),1)
  83. C
  84. C           ROW ELIMINATION WITH COLUMN INDEXING
  85. C
  86.             DO 30 J = KP1, N
  87.                T = A(L,J)
  88.                IF (L .EQ. K) GO TO 20
  89.                   A(L,J) = A(K,J)
  90.                   A(K,J) = T
  91.    20          CONTINUE
  92.                CALL DAXPY(N-K,T,A(K+1,K),1,A(K+1,J),1)
  93.    30       CONTINUE
  94.          GO TO 50
  95.    40    CONTINUE
  96.             INFO = K
  97.    50    CONTINUE
  98.    60 CONTINUE
  99.    70 CONTINUE
  100.       IPVT(N) = N
  101.       IF (A(N,N) .EQ. 0.0D0) INFO = N
  102.       RETURN
  103.       END
  104.