2
SUBROUTINE SCHDC (A, LDA, P, WORK, JPVT, JOB, INFO)
3
C***BEGIN PROLOGUE SCHDC
4
C***PURPOSE Compute the Cholesky decomposition of a positive definite
5
C matrix. A pivoting option allows the user to estimate the
6
C condition number of a positive definite matrix or determine
7
C the rank of a positive semidefinite matrix.
8
C***LIBRARY SLATEC (LINPACK)
10
C***TYPE SINGLE PRECISION (SCHDC-S, DCHDC-D, CCHDC-C)
11
C***KEYWORDS CHOLESKY DECOMPOSITION, LINEAR ALGEBRA, LINPACK, MATRIX,
13
C***AUTHOR Dongarra, J., (ANL)
14
C Stewart, G. W., (U. of Maryland)
17
C SCHDC computes the Cholesky decomposition of a positive definite
18
C matrix. A pivoting option allows the user to estimate the
19
C condition of a positive definite matrix or determine the rank
20
C of a positive semidefinite matrix.
25
C A contains the matrix whose decomposition is to
26
C be computed. Only the upper half of A need be stored.
27
C The lower part of the array A is not referenced.
30
C LDA is the leading dimension of the array A.
33
C P is the order of the matrix.
36
C WORK is a work array.
39
C JPVT contains integers that control the selection
40
C of the pivot elements, if pivoting has been requested.
41
C Each diagonal element A(K,K)
42
C is placed in one of three classes according to the
45
C If JPVT(K) .GT. 0, then X(K) is an initial
48
C If JPVT(K) .EQ. 0, then X(K) is a free element.
50
C If JPVT(K) .LT. 0, then X(K) is a final element.
52
C Before the decomposition is computed, initial elements
53
C are moved by symmetric row and column interchanges to
54
C the beginning of the array A and final
55
C elements to the end. Both initial and final elements
56
C are frozen in place during the computation and only
57
C free elements are moved. At the K-th stage of the
58
C reduction, if A(K,K) is occupied by a free element
59
C it is interchanged with the largest free element
60
C A(L,L) with L .GE. K. JPVT is not referenced if
64
C JOB is an integer that initiates column pivoting.
65
C If JOB .EQ. 0, no pivoting is done.
66
C If JOB .NE. 0, pivoting is done.
70
C A A contains in its upper half the Cholesky factor
71
C of the matrix A as it has been permuted by pivoting.
73
C JPVT JPVT(J) contains the index of the diagonal element
74
C of a that was moved into the J-th position,
75
C provided pivoting was requested.
77
C INFO contains the index of the last positive diagonal
78
C element of the Cholesky factor.
80
C For positive definite matrices INFO = P is the normal return.
81
C For pivoting with positive semidefinite matrices INFO will
82
C in general be less than P. However, INFO may be greater than
83
C the rank of A, since rounding error can cause an otherwise zero
84
C element to be positive. Indefinite systems will always cause
85
C INFO to be less than P.
87
C***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
88
C Stewart, LINPACK Users' Guide, SIAM, 1979.
89
C***ROUTINES CALLED SAXPY, SSWAP
90
C***REVISION HISTORY (YYMMDD)
92
C 890313 REVISION DATE from Version 3.2
93
C 891214 Prologue converted to Version 4.0 format. (BAB)
94
C 900326 Removed duplicate information from DESCRIPTION section.
96
C 920501 Reformatted the REFERENCES section. (WRB)
97
C***END PROLOGUE SCHDC
98
INTEGER LDA,P,JPVT(*),JOB,INFO
101
INTEGER PU,PL,PLP1,J,JP,JT,K,KB,KM1,KP1,L,MAXL
105
C***FIRST EXECUTABLE STATEMENT SCHDC
109
IF (JOB .EQ. 0) GO TO 160
111
C PIVOTING HAS BEEN REQUESTED. REARRANGE THE
112
C THE ELEMENTS ACCORDING TO JPVT.
115
SWAPK = JPVT(K) .GT. 0
116
NEGK = JPVT(K) .LT. 0
118
IF (NEGK) JPVT(K) = -JPVT(K)
119
IF (.NOT.SWAPK) GO TO 60
120
IF (K .EQ. PL) GO TO 50
121
CALL SSWAP(PL-1,A(1,K),1,A(1,PL),1)
126
IF (P .LT. PLP1) GO TO 40
128
IF (J .GE. K) GO TO 10
134
IF (J .EQ. K) GO TO 20
148
IF (P .LT. PL) GO TO 150
151
IF (JPVT(K) .GE. 0) GO TO 130
153
IF (PU .EQ. K) GO TO 120
154
CALL SSWAP(K-1,A(1,K),1,A(1,PU),1)
159
IF (P .LT. KP1) GO TO 110
161
IF (J .GE. PU) GO TO 80
167
IF (J .EQ. PU) GO TO 90
191
C DETERMINE THE PIVOT ELEMENT.
193
IF (K .LT. PL .OR. K .GE. PU) GO TO 190
195
IF (A(L,L) .LE. MAXDIA) GO TO 170
202
C QUIT IF THE PIVOT ELEMENT IS NOT POSITIVE.
204
IF (MAXDIA .GT. 0.0E0) GO TO 200
208
IF (K .EQ. MAXL) GO TO 210
210
C START THE PIVOTING AND UPDATE JPVT.
213
CALL SSWAP(KM1,A(1,K),1,A(1,MAXL),1)
214
A(MAXL,MAXL) = A(K,K)
221
C REDUCTION STEP. PIVOTING IS CONTAINED ACROSS THE ROWS.
223
WORK(K) = SQRT(A(K,K))
225
IF (P .LT. KP1) GO TO 260
227
IF (K .EQ. MAXL) GO TO 240
228
IF (J .GE. MAXL) GO TO 220
234
IF (J .EQ. MAXL) GO TO 230
240
A(K,J) = A(K,J)/WORK(K)
243
CALL SAXPY(J-K,TEMP,WORK(KP1),1,A(KP1,J),1)