1
2      SUBROUTINE SCHDC(A,LDA,P,WORK,JPVT,JOB,INFO)
3C***BEGIN PROLOGUE  SCHDC
4C***DATE WRITTEN   790319   (YYMMDD)
5C***REVISION DATE  820801   (YYMMDD)
6C***REVISION HISTORY  (YYMMDD)
7C   000330  Modified array declarations.  (JEC)
8C***CATEGORY NO.  D2B1B
9C***KEYWORDS  CHOLESKY DECOMPOSITION,LINEAR ALGEBRA,LINPACK,MATRIX,
10C             POSITIVE DEFINITE
11C***AUTHOR  DONGARRA, J., (ANL)
12C           STEWART, G. W., (U. OF MARYLAND)
13C***PURPOSE  Computes the Cholesky decomposition of A POSITIVE DEFINITE
14C            matrix.  A pivoting option allows the user to estimate the
15C            condition of a POSITIVE DEFINITE matrix or determine the
16C            rank of a POSITIVE SEMIDEFINITE matrix.
17C***DESCRIPTION
18C
19C     SCHDC computes the Cholesky decomposition of a positive definite
20C     matrix.  A pivoting option allows the user to estimate the
21C     condition of a positive definite matrix or determine the rank
22C     of a positive semidefinite matrix.
23C
24C     On Entry
25C
26C         A      REAL(LDA,P).
27C                A contains the matrix whose decomposition is to
28C                be computed.  Only the upper half of A need be stored.
29C                The lower part of the array A is not referenced.
30C
31C         LDA    INTEGER.
32C                LDA is the leading dimension of the array A.
33C
34C         P      INTEGER.
35C                P is the order of the matrix.
36C
37C         WORK   REAL.
38C                WORK is a work array.
39C
40C         JPVT   INTEGER(P).
41C                JPVT contains integers that control the selection
42C                of the pivot elements, if pivoting has been requested.
43C                Each diagonal element A(K,K)
44C                is placed in one of three classes according to the
45C                value of JPVT(K).
46C
47C                   If JPVT(K) .GT. 0, then X(K) is an initial
48C                                      element.
49C
50C                   If JPVT(K) .EQ. 0, then X(K) is a free element.
51C
52C                   If JPVT(K) .LT. 0, then X(K) is a final element.
53C
54C                Before the decomposition is computed, initial elements
55C                are moved by symmetric row and column interchanges to
56C                the beginning of the array A and final
57C                elements to the end.  Both initial and final elements
58C                are frozen in place during the computation and only
59C                free elements are moved.  At the K-th stage of the
60C                reduction, if A(K,K) is occupied by a free element
61C                it is interchanged with the largest free element
62C                A(L,L) with L .GE. K.  JPVT is not referenced if
63C                JOB .EQ. 0.
64C
65C        JOB     INTEGER.
66C                JOB is an integer that initiates column pivoting.
67C                If JOB .EQ. 0, no pivoting is done.
68C                If JOB .NE. 0, pivoting is done.
69C
70C     On Return
71C
72C         A      A contains in its upper half the Cholesky factor
73C                of the matrix A as it has been permuted by pivoting.
74C
75C         JPVT   JPVT(J) contains the index of the diagonal element
76C                of a that was moved into the J-th position,
77C                provided pivoting was requested.
78C
79C         INFO   contains the index of the last positive diagonal
80C                element of the Cholesky factor.
81C
82C     For positive definite matrices INFO = P is the normal return.
83C     For pivoting with positive semidefinite matrices INFO will
84C     in general be less than P.  However, INFO may be greater than
85C     the rank of A, since rounding error can cause an otherwise zero
86C     element to be positive.  Indefinite systems will always cause
87C     INFO to be less than P.
88C
89C     LINPACK.  This version dated 03/19/79 .
90C     J. J. Dongarra and G. W. Stewart, Argonne National Laboratory and
91C     University of Maryland.
92C
93C
94C     BLAS SAXPY,SSWAP
95C     Fortran SQRT
96C***REFERENCES  DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W.,
97C                 *LINPACK USERS  GUIDE*, SIAM, 1979.
98C***ROUTINES CALLED  SAXPY,SSWAP
99C***END PROLOGUE  SCHDC
100      INTEGER LDA,P,JPVT(*),JOB,INFO
101      REAL A(LDA,*),WORK(*)
102C
103      INTEGER PU,PL,PLP1,I,J,JP,JT,K,KB,KM1,KP1,L,MAXL
104      REAL TEMP
105      REAL MAXDIA
106      LOGICAL SWAPK,NEGK
107C***FIRST EXECUTABLE STATEMENT  SCHDC
108      PL = 1
109      PU = 0
110      INFO = P
111      IF (JOB .EQ. 0) GO TO 160
112C
113C        PIVOTING HAS BEEN REQUESTED. REARRANGE THE
114C        THE ELEMENTS ACCORDING TO JPVT.
115C
116         DO 70 K = 1, P
117            SWAPK = JPVT(K) .GT. 0
118            NEGK = JPVT(K) .LT. 0
119            JPVT(K) = K
120            IF (NEGK) JPVT(K) = -JPVT(K)
121            IF (.NOT.SWAPK) GO TO 60
122               IF (K .EQ. PL) GO TO 50
123                  CALL SSWAP(PL-1,A(1,K),1,A(1,PL),1)
124                  TEMP = A(K,K)
125                  A(K,K) = A(PL,PL)
126                  A(PL,PL) = TEMP
127                  PLP1 = PL + 1
128                  IF (P .LT. PLP1) GO TO 40
129                  DO 30 J = PLP1, P
130                     IF (J .GE. K) GO TO 10
131                        TEMP = A(PL,J)
132                        A(PL,J) = A(J,K)
133                        A(J,K) = TEMP
134                     GO TO 20
135   10                CONTINUE
136                     IF (J .EQ. K) GO TO 20
137                        TEMP = A(K,J)
138                        A(K,J) = A(PL,J)
139                        A(PL,J) = TEMP
140   20                CONTINUE
141   30             CONTINUE
142   40             CONTINUE
143                  JPVT(K) = JPVT(PL)
144                  JPVT(PL) = K
145   50          CONTINUE
146               PL = PL + 1
147   60       CONTINUE
148   70    CONTINUE
149         PU = P
150         IF (P .LT. PL) GO TO 150
151         DO 140 KB = PL, P
152            K = P - KB + PL
153            IF (JPVT(K) .GE. 0) GO TO 130
154               JPVT(K) = -JPVT(K)
155               IF (PU .EQ. K) GO TO 120
156                  CALL SSWAP(K-1,A(1,K),1,A(1,PU),1)
157                  TEMP = A(K,K)
158                  A(K,K) = A(PU,PU)
159                  A(PU,PU) = TEMP
160                  KP1 = K + 1
161                  IF (P .LT. KP1) GO TO 110
162                  DO 100 J = KP1, P
163                     IF (J .GE. PU) GO TO 80
164                        TEMP = A(K,J)
165                        A(K,J) = A(J,PU)
166                        A(J,PU) = TEMP
167                     GO TO 90
168   80                CONTINUE
169                     IF (J .EQ. PU) GO TO 90
170                        TEMP = A(K,J)
171                        A(K,J) = A(PU,J)
172                        A(PU,J) = TEMP
173   90                CONTINUE
174  100             CONTINUE
175  110             CONTINUE
176                  JT = JPVT(K)
177                  JPVT(K) = JPVT(PU)
178                  JPVT(PU) = JT
179  120          CONTINUE
180               PU = PU - 1
181  130       CONTINUE
182  140    CONTINUE
183  150    CONTINUE
184  160 CONTINUE
185      DO 270 K = 1, P
186C
187C        REDUCTION LOOP.
188C
189         MAXDIA = A(K,K)
190         KP1 = K + 1
191         MAXL = K
192C
193C        DETERMINE THE PIVOT ELEMENT.
194C
195         IF (K .LT. PL .OR. K .GE. PU) GO TO 190
196            DO 180 L = KP1, PU
197               IF (A(L,L) .LE. MAXDIA) GO TO 170
198                  MAXDIA = A(L,L)
199                  MAXL = L
200  170          CONTINUE
201  180       CONTINUE
202  190    CONTINUE
203C
204C        QUIT IF THE PIVOT ELEMENT IS NOT POSITIVE.
205C
206         IF (MAXDIA .GT. 0.0E0) GO TO 200
207            INFO = K - 1
208C     ......EXIT
209            GO TO 280
210  200    CONTINUE
211         IF (K .EQ. MAXL) GO TO 210
212C
213C           START THE PIVOTING AND UPDATE JPVT.
214C
215            KM1 = K - 1
216            CALL SSWAP(KM1,A(1,K),1,A(1,MAXL),1)
217            A(MAXL,MAXL) = A(K,K)
218            A(K,K) = MAXDIA
219            JP = JPVT(MAXL)
220            JPVT(MAXL) = JPVT(K)
221            JPVT(K) = JP
222  210    CONTINUE
223C
224C        REDUCTION STEP. PIVOTING IS CONTAINED ACROSS THE ROWS.
225C
226         WORK(K) = SQRT(A(K,K))
227         A(K,K) = WORK(K)
228         IF (P .LT. KP1) GO TO 260
229         DO 250 J = KP1, P
230            IF (K .EQ. MAXL) GO TO 240
231               IF (J .GE. MAXL) GO TO 220
232                  TEMP = A(K,J)
233                  A(K,J) = A(J,MAXL)
234                  A(J,MAXL) = TEMP
235               GO TO 230
236  220          CONTINUE
237               IF (J .EQ. MAXL) GO TO 230
238                  TEMP = A(K,J)
239                  A(K,J) = A(MAXL,J)
240                  A(MAXL,J) = TEMP
241  230          CONTINUE
242  240       CONTINUE
243            A(K,J) = A(K,J)/WORK(K)
244            WORK(J) = A(K,J)
245            TEMP = -A(K,J)
246            CALL SAXPY(J-K,TEMP,WORK(KP1),1,A(KP1,J),1)
247  250    CONTINUE
248  260    CONTINUE
249  270 CONTINUE
250  280 CONTINUE
251      RETURN
252      END
253