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