home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / fortran / library / linpack / schdc.for < prev    next >
Encoding:
Text File  |  1985-01-13  |  7.7 KB  |  235 lines
  1.       SUBROUTINE SCHDC(A,LDA,P,WORK,JPVT,JOB,INFO)
  2.       INTEGER LDA,P,JPVT(1),JOB,INFO
  3.       REAL A(LDA,1),WORK(1)
  4. C
  5. C     SCHDC COMPUTES THE CHOLESKY DECOMPOSITION OF A POSITIVE DEFINITE
  6. C     MATRIX.  A PIVOTING OPTION ALLOWS THE USER TO ESTIMATE THE
  7. C     CONDITION OF A POSITIVE DEFINITE MATRIX OR DETERMINE THE RANK
  8. C     OF A POSITIVE SEMIDEFINITE MATRIX.
  9. C
  10. C     ON ENTRY
  11. C
  12. C         A      REAL(LDA,P).
  13. C                A CONTAINS THE MATRIX WHOSE DECOMPOSITION IS TO
  14. C                BE COMPUTED.  ONLT THE UPPER HALF OF A NEED BE STORED.
  15. C                THE LOWER PART OF THE ARRAY A IS NOT REFERENCED.
  16. C
  17. C         LDA    INTEGER.
  18. C                LDA IS THE LEADING DIMENSION OF THE ARRAY A.
  19. C
  20. C         P      INTEGER.
  21. C                P IS THE ORDER OF THE MATRIX.
  22. C
  23. C         WORK   REAL.
  24. C                WORK IS A WORK ARRAY.
  25. C
  26. C         JPVT   INTEGER(P).
  27. C                JPVT CONTAINS INTEGERS THAT CONTROL THE SELECTION
  28. C                OF THE PIVOT ELEMENTS, IF PIVOTING HAS BEEN REQUESTED.
  29. C                EACH DIAGONAL ELEMENT A(K,K)
  30. C                IS PLACED IN ONE OF THREE CLASSES ACCORDING TO THE
  31. C                VALUE OF JPVT(K).
  32. C
  33. C                   IF JPVT(K) .GT. 0, THEN X(K) IS AN INITIAL
  34. C                                      ELEMENT.
  35. C
  36. C                   IF JPVT(K) .EQ. 0, THEN X(K) IS A FREE ELEMENT.
  37. C
  38. C                   IF JPVT(K) .LT. 0, THEN X(K) IS A FINAL ELEMENT.
  39. C
  40. C                BEFORE THE DECOMPOSITION IS COMPUTED, INITIAL ELEMENTS
  41. C                ARE MOVED BY SYMMETRIC ROW AND COLUMN INTERCHANGES TO
  42. C                THE BEGINNING OF THE ARRAY A AND FINAL
  43. C                ELEMENTS TO THE END.  BOTH INITIAL AND FINAL ELEMENTS
  44. C                ARE FROZEN IN PLACE DURING THE COMPUTATION AND ONLY
  45. C                FREE ELEMENTS ARE MOVED.  AT THE K-TH STAGE OF THE
  46. C                REDUCTION, IF A(K,K) IS OCCUPIED BY A FREE ELEMENT
  47. C                IT IS INTERCHANGED WITH THE LARGEST FREE ELEMENT
  48. C                A(L,L) WITH L .GE. K.  JPVT IS NOT REFERENCED IF
  49. C                JOB .EQ. 0.
  50. C
  51. C        JOB     INTEGER.
  52. C                JOB IS AN INTEGER THAT INITIATES COLUMN PIVOTING.
  53. C                IF JOB .EQ. 0, NO PIVOTING IS DONE.
  54. C                IF JOB .NE. 0, PIVOTING IS DONE.
  55. C
  56. C     ON RETURN
  57. C
  58. C         A      A CONTAINS IN ITS UPPER HALF THE CHOLESKY FACTOR
  59. C                OF THE MATRIX A AS IT HAS BEEN PERMUTED BY PIVOTING.
  60. C
  61. C         JPVT   JPVT(J) CONTAINS THE INDEX OF THE DIAGONAL ELEMENT
  62. C                OF A THAT WAS MOVED INTO THE J-TH POSITION,
  63. C                PROVIDED PIVOTING WAS REQUESTED.
  64. C
  65. C         INFO   CONTAINS THE INDEX OF THE LAST POSITIVE DIAGONAL
  66. C                ELEMENT OF THE CHOLESKY FACTOR.
  67. C
  68. C     FOR POSITIVE DEFINITE MATRICES INFO = P IS THE NORMAL RETURN.
  69. C     FOR PIVOTING WITH POSITIVE SEMIDEFINITE MATRICES INFO WILL
  70. C     IN GENERAL BE LESS THAN P.  HOWEVER, INFO MAY BE GREATER THAN
  71. C     THE RANK OF A, SINCE ROUNDING ERROR CAN CAUSE AN OTHERWISE ZERO
  72. C     ELEMENT TO BE POSITIVE. INDEFINITE SYSTEMS WILL ALWAYS CAUSE
  73. C     INFO TO BE LESS THAN P.
  74. C
  75. C     LINPACK. THIS VERSION DATED 03/19/79 .
  76. C     J.J. DONGARRA AND G.W. STEWART, ARGONNE NATIONAL LABORATORY AND
  77. C     UNIVERSITY OF MARYLAND.
  78. C
  79. C
  80. C     BLAS SAXPY,SSWAP
  81. C     FORTRAN SQRT
  82. C
  83. C     INTERNAL VARIABLES
  84. C
  85.       INTEGER PU,PL,PLP1,I,J,JP,JT,K,KB,KM1,KP1,L,MAXL
  86.       REAL TEMP
  87.       REAL MAXDIA
  88.       LOGICAL SWAPK,NEGK
  89. C
  90.       PL = 1
  91.       PU = 0
  92.       INFO = P
  93.       IF (JOB .EQ. 0) GO TO 160
  94. C
  95. C        PIVOTING HAS BEEN REQUESTED. REARRANGE THE
  96. C        THE ELEMENTS ACCORDING TO JPVT.
  97. C
  98.          DO 70 K = 1, P
  99.             SWAPK = JPVT(K) .GT. 0
  100.             NEGK = JPVT(K) .LT. 0
  101.             JPVT(K) = K
  102.             IF (NEGK) JPVT(K) = -JPVT(K)
  103.             IF (.NOT.SWAPK) GO TO 60
  104.                IF (K .EQ. PL) GO TO 50
  105.                   CALL SSWAP(PL-1,A(1,K),1,A(1,PL),1)
  106.                   TEMP = A(K,K)
  107.                   A(K,K) = A(PL,PL)
  108.                   A(PL,PL) = TEMP
  109.                   PLP1 = PL + 1
  110.                   IF (P .LT. PLP1) GO TO 40
  111.                   DO 30 J = PLP1, P
  112.                      IF (J .GE. K) GO TO 10
  113.                         TEMP = A(PL,J)
  114.                         A(PL,J) = A(J,K)
  115.                         A(J,K) = TEMP
  116.                      GO TO 20
  117.    10                CONTINUE
  118.                      IF (J .EQ. K) GO TO 20
  119.                         TEMP = A(K,J)
  120.                         A(K,J) = A(PL,J)
  121.                         A(PL,J) = TEMP
  122.    20                CONTINUE
  123.    30             CONTINUE
  124.    40             CONTINUE
  125.                   JPVT(K) = JPVT(PL)
  126.                   JPVT(PL) = K
  127.    50          CONTINUE
  128.                PL = PL + 1
  129.    60       CONTINUE
  130.    70    CONTINUE
  131.          PU = P
  132.          IF (P .LT. PL) GO TO 150
  133.          DO 140 KB = PL, P
  134.             K = P - KB + PL
  135.             IF (JPVT(K) .GE. 0) GO TO 130
  136.                JPVT(K) = -JPVT(K)
  137.                IF (PU .EQ. K) GO TO 120
  138.                   CALL SSWAP(K-1,A(1,K),1,A(1,PU),1)
  139.                   TEMP = A(K,K)
  140.                   A(K,K) = A(PU,PU)
  141.                   A(PU,PU) = TEMP
  142.                   KP1 = K + 1
  143.                   IF (P .LT. KP1) GO TO 110
  144.                   DO 100 J = KP1, P
  145.                      IF (J .GE. PU) GO TO 80
  146.                         TEMP = A(K,J)
  147.                         A(K,J) = A(J,PU)
  148.                         A(J,PU) = TEMP
  149.                      GO TO 90
  150.    80                CONTINUE
  151.                      IF (J .EQ. PU) GO TO 90
  152.                         TEMP = A(K,J)
  153.                         A(K,J) = A(PU,J)
  154.                         A(PU,J) = TEMP
  155.    90                CONTINUE
  156.   100             CONTINUE
  157.   110             CONTINUE
  158.                   JT = JPVT(K)
  159.                   JPVT(K) = JPVT(PU)
  160.                   JPVT(PU) = JT
  161.   120          CONTINUE
  162.                PU = PU - 1
  163.   130       CONTINUE
  164.   140    CONTINUE
  165.   150    CONTINUE
  166.   160 CONTINUE
  167.       DO 270 K = 1, P
  168. C
  169. C        REDUCTION LOOP.
  170. C
  171.          MAXDIA = A(K,K)
  172.          KP1 = K + 1
  173.          MAXL = K
  174. C
  175. C        DETERMINE THE PIVOT ELEMENT.
  176. C
  177.          IF (K .LT. PL .OR. K .GE. PU) GO TO 190
  178.             DO 180 L = KP1, PU
  179.                IF (A(L,L) .LE. MAXDIA) GO TO 170
  180.                   MAXDIA = A(L,L)
  181.                   MAXL = L
  182.   170          CONTINUE
  183.   180       CONTINUE
  184.   190    CONTINUE
  185. C
  186. C        QUIT IF THE PIVOT ELEMENT IS NOT POSITIVE.
  187. C
  188.          IF (MAXDIA .GT. 0.0E0) GO TO 200
  189.             INFO = K - 1
  190. C     ......EXIT
  191.             GO TO 280
  192.   200    CONTINUE
  193.          IF (K .EQ. MAXL) GO TO 210
  194. C
  195. C           START THE PIVOTING AND UPDATE JPVT.
  196. C
  197.             KM1 = K - 1
  198.             CALL SSWAP(KM1,A(1,K),1,A(1,MAXL),1)
  199.             A(MAXL,MAXL) = A(K,K)
  200.             A(K,K) = MAXDIA
  201.             JP = JPVT(MAXL)
  202.             JPVT(MAXL) = JPVT(K)
  203.             JPVT(K) = JP
  204.   210    CONTINUE
  205. C
  206. C        REDUCTION STEP. PIVOTING IS CONTAINED ACROSS THE ROWS.
  207. C
  208.          WORK(K) = SQRT(A(K,K))
  209.          A(K,K) = WORK(K)
  210.          IF (P .LT. KP1) GO TO 260
  211.          DO 250 J = KP1, P
  212.             IF (K .EQ. MAXL) GO TO 240
  213.                IF (J .GE. MAXL) GO TO 220
  214.                   TEMP = A(K,J)
  215.                   A(K,J) = A(J,MAXL)
  216.                   A(J,MAXL) = TEMP
  217.                GO TO 230
  218.   220          CONTINUE
  219.                IF (J .EQ. MAXL) GO TO 230
  220.                   TEMP = A(K,J)
  221.                   A(K,J) = A(MAXL,J)
  222.                   A(MAXL,J) = TEMP
  223.   230          CONTINUE
  224.   240       CONTINUE
  225.             A(K,J) = A(K,J)/WORK(K)
  226.             WORK(J) = A(K,J)
  227.             TEMP = -A(K,J)
  228.             CALL SAXPY(J-K,TEMP,WORK(KP1),1,A(KP1,J),1)
  229.   250    CONTINUE
  230.   260    CONTINUE
  231.   270 CONTINUE
  232.   280 CONTINUE
  233.       RETURN
  234.       END
  235.