1 SUBROUTINE DPSTF2( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
2 *
3 * -- LAPACK PROTOTYPE routine (version 3.2.2) --
4 * Craig Lucas, University of Manchester / NAG Ltd.
5 * October, 2008
6 *
7 * .. Scalar Arguments ..
8 DOUBLE PRECISION TOL
9 INTEGER INFO, LDA, N, RANK
10 CHARACTER UPLO
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION A( LDA, * ), WORK( 2*N )
14 INTEGER PIV( N )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DPSTF2 computes the Cholesky factorization with complete
21 * pivoting of a real symmetric positive semidefinite matrix A.
22 *
23 * The factorization has the form
24 * P**T * A * P = U**T * U , if UPLO = 'U',
25 * P**T * A * P = L * L**T, if UPLO = 'L',
26 * where U is an upper triangular matrix and L is lower triangular, and
27 * P is stored as vector PIV.
28 *
29 * This algorithm does not attempt to check that A is positive
30 * semidefinite. This version of the algorithm calls level 2 BLAS.
31 *
32 * Arguments
33 * =========
34 *
35 * UPLO (input) CHARACTER*1
36 * Specifies whether the upper or lower triangular part of the
37 * symmetric matrix A is stored.
38 * = 'U': Upper triangular
39 * = 'L': Lower triangular
40 *
41 * N (input) INTEGER
42 * The order of the matrix A. N >= 0.
43 *
44 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
45 * On entry, the symmetric matrix A. If UPLO = 'U', the leading
46 * n by n upper triangular part of A contains the upper
47 * triangular part of the matrix A, and the strictly lower
48 * triangular part of A is not referenced. If UPLO = 'L', the
49 * leading n by n lower triangular part of A contains the lower
50 * triangular part of the matrix A, and the strictly upper
51 * triangular part of A is not referenced.
52 *
53 * On exit, if INFO = 0, the factor U or L from the Cholesky
54 * factorization as above.
55 *
56 * PIV (output) INTEGER array, dimension (N)
57 * PIV is such that the nonzero entries are P( PIV(K), K ) = 1.
58 *
59 * RANK (output) INTEGER
60 * The rank of A given by the number of steps the algorithm
61 * completed.
62 *
63 * TOL (input) DOUBLE PRECISION
64 * User defined tolerance. If TOL < 0, then N*U*MAX( A( K,K ) )
65 * will be used. The algorithm terminates at the (K-1)st step
66 * if the pivot <= TOL.
67 *
68 * LDA (input) INTEGER
69 * The leading dimension of the array A. LDA >= max(1,N).
70 *
71 * WORK (workspace) DOUBLE PRECISION array, dimension (2*N)
72 * Work space.
73 *
74 * INFO (output) INTEGER
75 * < 0: If INFO = -K, the K-th argument had an illegal value,
76 * = 0: algorithm completed successfully, and
77 * > 0: the matrix A is either rank deficient with computed rank
78 * as returned in RANK, or is indefinite. See Section 7 of
79 * LAPACK Working Note #161 for further information.
80 *
81 * =====================================================================
82 *
83 * .. Parameters ..
84 DOUBLE PRECISION ONE, ZERO
85 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
86 * ..
87 * .. Local Scalars ..
88 DOUBLE PRECISION AJJ, DSTOP, DTEMP
89 INTEGER I, ITEMP, J, PVT
90 LOGICAL UPPER
91 * ..
92 * .. External Functions ..
93 DOUBLE PRECISION DLAMCH
94 LOGICAL LSAME, DISNAN
95 EXTERNAL DLAMCH, LSAME, DISNAN
96 * ..
97 * .. External Subroutines ..
98 EXTERNAL DGEMV, DSCAL, DSWAP, XERBLA
99 * ..
100 * .. Intrinsic Functions ..
101 INTRINSIC MAX, SQRT, MAXLOC
102 * ..
103 * .. Executable Statements ..
104 *
105 * Test the input parameters
106 *
107 INFO = 0
108 UPPER = LSAME( UPLO, 'U' )
109 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
110 INFO = -1
111 ELSE IF( N.LT.0 ) THEN
112 INFO = -2
113 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
114 INFO = -4
115 END IF
116 IF( INFO.NE.0 ) THEN
117 CALL XERBLA( 'DPSTF2', -INFO )
118 RETURN
119 END IF
120 *
121 * Quick return if possible
122 *
123 IF( N.EQ.0 )
124 $ RETURN
125 *
126 * Initialize PIV
127 *
128 DO 100 I = 1, N
129 PIV( I ) = I
130 100 CONTINUE
131 *
132 * Compute stopping value
133 *
134 PVT = 1
135 AJJ = A( PVT, PVT )
136 DO I = 2, N
137 IF( A( I, I ).GT.AJJ ) THEN
138 PVT = I
139 AJJ = A( PVT, PVT )
140 END IF
141 END DO
142 IF( AJJ.EQ.ZERO.OR.DISNAN( AJJ ) ) THEN
143 RANK = 0
144 INFO = 1
145 GO TO 170
146 END IF
147 *
148 * Compute stopping value if not supplied
149 *
150 IF( TOL.LT.ZERO ) THEN
151 DSTOP = N * DLAMCH( 'Epsilon' ) * AJJ
152 ELSE
153 DSTOP = TOL
154 END IF
155 *
156 * Set first half of WORK to zero, holds dot products
157 *
158 DO 110 I = 1, N
159 WORK( I ) = 0
160 110 CONTINUE
161 *
162 IF( UPPER ) THEN
163 *
164 * Compute the Cholesky factorization P**T * A * P = U**T * U
165 *
166 DO 130 J = 1, N
167 *
168 * Find pivot, test for exit, else swap rows and columns
169 * Update dot products, compute possible pivots which are
170 * stored in the second half of WORK
171 *
172 DO 120 I = J, N
173 *
174 IF( J.GT.1 ) THEN
175 WORK( I ) = WORK( I ) + A( J-1, I )**2
176 END IF
177 WORK( N+I ) = A( I, I ) - WORK( I )
178 *
179 120 CONTINUE
180 *
181 IF( J.GT.1 ) THEN
182 ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
183 PVT = ITEMP + J - 1
184 AJJ = WORK( N+PVT )
185 IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN
186 A( J, J ) = AJJ
187 GO TO 160
188 END IF
189 END IF
190 *
191 IF( J.NE.PVT ) THEN
192 *
193 * Pivot OK, so can now swap pivot rows and columns
194 *
195 A( PVT, PVT ) = A( J, J )
196 CALL DSWAP( J-1, A( 1, J ), 1, A( 1, PVT ), 1 )
197 IF( PVT.LT.N )
198 $ CALL DSWAP( N-PVT, A( J, PVT+1 ), LDA,
199 $ A( PVT, PVT+1 ), LDA )
200 CALL DSWAP( PVT-J-1, A( J, J+1 ), LDA, A( J+1, PVT ), 1 )
201 *
202 * Swap dot products and PIV
203 *
204 DTEMP = WORK( J )
205 WORK( J ) = WORK( PVT )
206 WORK( PVT ) = DTEMP
207 ITEMP = PIV( PVT )
208 PIV( PVT ) = PIV( J )
209 PIV( J ) = ITEMP
210 END IF
211 *
212 AJJ = SQRT( AJJ )
213 A( J, J ) = AJJ
214 *
215 * Compute elements J+1:N of row J
216 *
217 IF( J.LT.N ) THEN
218 CALL DGEMV( 'Trans', J-1, N-J, -ONE, A( 1, J+1 ), LDA,
219 $ A( 1, J ), 1, ONE, A( J, J+1 ), LDA )
220 CALL DSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA )
221 END IF
222 *
223 130 CONTINUE
224 *
225 ELSE
226 *
227 * Compute the Cholesky factorization P**T * A * P = L * L**T
228 *
229 DO 150 J = 1, N
230 *
231 * Find pivot, test for exit, else swap rows and columns
232 * Update dot products, compute possible pivots which are
233 * stored in the second half of WORK
234 *
235 DO 140 I = J, N
236 *
237 IF( J.GT.1 ) THEN
238 WORK( I ) = WORK( I ) + A( I, J-1 )**2
239 END IF
240 WORK( N+I ) = A( I, I ) - WORK( I )
241 *
242 140 CONTINUE
243 *
244 IF( J.GT.1 ) THEN
245 ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
246 PVT = ITEMP + J - 1
247 AJJ = WORK( N+PVT )
248 IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN
249 A( J, J ) = AJJ
250 GO TO 160
251 END IF
252 END IF
253 *
254 IF( J.NE.PVT ) THEN
255 *
256 * Pivot OK, so can now swap pivot rows and columns
257 *
258 A( PVT, PVT ) = A( J, J )
259 CALL DSWAP( J-1, A( J, 1 ), LDA, A( PVT, 1 ), LDA )
260 IF( PVT.LT.N )
261 $ CALL DSWAP( N-PVT, A( PVT+1, J ), 1, A( PVT+1, PVT ),
262 $ 1 )
263 CALL DSWAP( PVT-J-1, A( J+1, J ), 1, A( PVT, J+1 ), LDA )
264 *
265 * Swap dot products and PIV
266 *
267 DTEMP = WORK( J )
268 WORK( J ) = WORK( PVT )
269 WORK( PVT ) = DTEMP
270 ITEMP = PIV( PVT )
271 PIV( PVT ) = PIV( J )
272 PIV( J ) = ITEMP
273 END IF
274 *
275 AJJ = SQRT( AJJ )
276 A( J, J ) = AJJ
277 *
278 * Compute elements J+1:N of column J
279 *
280 IF( J.LT.N ) THEN
281 CALL DGEMV( 'No Trans', N-J, J-1, -ONE, A( J+1, 1 ), LDA,
282 $ A( J, 1 ), LDA, ONE, A( J+1, J ), 1 )
283 CALL DSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 )
284 END IF
285 *
286 150 CONTINUE
287 *
288 END IF
289 *
290 * Ran to completion, A has full rank
291 *
292 RANK = N
293 *
294 GO TO 170
295 160 CONTINUE
296 *
297 * Rank is number of steps completed. Set INFO = 1 to signal
298 * that the factorization cannot be used to solve a system.
299 *
300 RANK = J - 1
301 INFO = 1
302 *
303 170 CONTINUE
304 RETURN
305 *
306 * End of DPSTF2
307 *
308 END
2 *
3 * -- LAPACK PROTOTYPE routine (version 3.2.2) --
4 * Craig Lucas, University of Manchester / NAG Ltd.
5 * October, 2008
6 *
7 * .. Scalar Arguments ..
8 DOUBLE PRECISION TOL
9 INTEGER INFO, LDA, N, RANK
10 CHARACTER UPLO
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION A( LDA, * ), WORK( 2*N )
14 INTEGER PIV( N )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DPSTF2 computes the Cholesky factorization with complete
21 * pivoting of a real symmetric positive semidefinite matrix A.
22 *
23 * The factorization has the form
24 * P**T * A * P = U**T * U , if UPLO = 'U',
25 * P**T * A * P = L * L**T, if UPLO = 'L',
26 * where U is an upper triangular matrix and L is lower triangular, and
27 * P is stored as vector PIV.
28 *
29 * This algorithm does not attempt to check that A is positive
30 * semidefinite. This version of the algorithm calls level 2 BLAS.
31 *
32 * Arguments
33 * =========
34 *
35 * UPLO (input) CHARACTER*1
36 * Specifies whether the upper or lower triangular part of the
37 * symmetric matrix A is stored.
38 * = 'U': Upper triangular
39 * = 'L': Lower triangular
40 *
41 * N (input) INTEGER
42 * The order of the matrix A. N >= 0.
43 *
44 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
45 * On entry, the symmetric matrix A. If UPLO = 'U', the leading
46 * n by n upper triangular part of A contains the upper
47 * triangular part of the matrix A, and the strictly lower
48 * triangular part of A is not referenced. If UPLO = 'L', the
49 * leading n by n lower triangular part of A contains the lower
50 * triangular part of the matrix A, and the strictly upper
51 * triangular part of A is not referenced.
52 *
53 * On exit, if INFO = 0, the factor U or L from the Cholesky
54 * factorization as above.
55 *
56 * PIV (output) INTEGER array, dimension (N)
57 * PIV is such that the nonzero entries are P( PIV(K), K ) = 1.
58 *
59 * RANK (output) INTEGER
60 * The rank of A given by the number of steps the algorithm
61 * completed.
62 *
63 * TOL (input) DOUBLE PRECISION
64 * User defined tolerance. If TOL < 0, then N*U*MAX( A( K,K ) )
65 * will be used. The algorithm terminates at the (K-1)st step
66 * if the pivot <= TOL.
67 *
68 * LDA (input) INTEGER
69 * The leading dimension of the array A. LDA >= max(1,N).
70 *
71 * WORK (workspace) DOUBLE PRECISION array, dimension (2*N)
72 * Work space.
73 *
74 * INFO (output) INTEGER
75 * < 0: If INFO = -K, the K-th argument had an illegal value,
76 * = 0: algorithm completed successfully, and
77 * > 0: the matrix A is either rank deficient with computed rank
78 * as returned in RANK, or is indefinite. See Section 7 of
79 * LAPACK Working Note #161 for further information.
80 *
81 * =====================================================================
82 *
83 * .. Parameters ..
84 DOUBLE PRECISION ONE, ZERO
85 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
86 * ..
87 * .. Local Scalars ..
88 DOUBLE PRECISION AJJ, DSTOP, DTEMP
89 INTEGER I, ITEMP, J, PVT
90 LOGICAL UPPER
91 * ..
92 * .. External Functions ..
93 DOUBLE PRECISION DLAMCH
94 LOGICAL LSAME, DISNAN
95 EXTERNAL DLAMCH, LSAME, DISNAN
96 * ..
97 * .. External Subroutines ..
98 EXTERNAL DGEMV, DSCAL, DSWAP, XERBLA
99 * ..
100 * .. Intrinsic Functions ..
101 INTRINSIC MAX, SQRT, MAXLOC
102 * ..
103 * .. Executable Statements ..
104 *
105 * Test the input parameters
106 *
107 INFO = 0
108 UPPER = LSAME( UPLO, 'U' )
109 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
110 INFO = -1
111 ELSE IF( N.LT.0 ) THEN
112 INFO = -2
113 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
114 INFO = -4
115 END IF
116 IF( INFO.NE.0 ) THEN
117 CALL XERBLA( 'DPSTF2', -INFO )
118 RETURN
119 END IF
120 *
121 * Quick return if possible
122 *
123 IF( N.EQ.0 )
124 $ RETURN
125 *
126 * Initialize PIV
127 *
128 DO 100 I = 1, N
129 PIV( I ) = I
130 100 CONTINUE
131 *
132 * Compute stopping value
133 *
134 PVT = 1
135 AJJ = A( PVT, PVT )
136 DO I = 2, N
137 IF( A( I, I ).GT.AJJ ) THEN
138 PVT = I
139 AJJ = A( PVT, PVT )
140 END IF
141 END DO
142 IF( AJJ.EQ.ZERO.OR.DISNAN( AJJ ) ) THEN
143 RANK = 0
144 INFO = 1
145 GO TO 170
146 END IF
147 *
148 * Compute stopping value if not supplied
149 *
150 IF( TOL.LT.ZERO ) THEN
151 DSTOP = N * DLAMCH( 'Epsilon' ) * AJJ
152 ELSE
153 DSTOP = TOL
154 END IF
155 *
156 * Set first half of WORK to zero, holds dot products
157 *
158 DO 110 I = 1, N
159 WORK( I ) = 0
160 110 CONTINUE
161 *
162 IF( UPPER ) THEN
163 *
164 * Compute the Cholesky factorization P**T * A * P = U**T * U
165 *
166 DO 130 J = 1, N
167 *
168 * Find pivot, test for exit, else swap rows and columns
169 * Update dot products, compute possible pivots which are
170 * stored in the second half of WORK
171 *
172 DO 120 I = J, N
173 *
174 IF( J.GT.1 ) THEN
175 WORK( I ) = WORK( I ) + A( J-1, I )**2
176 END IF
177 WORK( N+I ) = A( I, I ) - WORK( I )
178 *
179 120 CONTINUE
180 *
181 IF( J.GT.1 ) THEN
182 ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
183 PVT = ITEMP + J - 1
184 AJJ = WORK( N+PVT )
185 IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN
186 A( J, J ) = AJJ
187 GO TO 160
188 END IF
189 END IF
190 *
191 IF( J.NE.PVT ) THEN
192 *
193 * Pivot OK, so can now swap pivot rows and columns
194 *
195 A( PVT, PVT ) = A( J, J )
196 CALL DSWAP( J-1, A( 1, J ), 1, A( 1, PVT ), 1 )
197 IF( PVT.LT.N )
198 $ CALL DSWAP( N-PVT, A( J, PVT+1 ), LDA,
199 $ A( PVT, PVT+1 ), LDA )
200 CALL DSWAP( PVT-J-1, A( J, J+1 ), LDA, A( J+1, PVT ), 1 )
201 *
202 * Swap dot products and PIV
203 *
204 DTEMP = WORK( J )
205 WORK( J ) = WORK( PVT )
206 WORK( PVT ) = DTEMP
207 ITEMP = PIV( PVT )
208 PIV( PVT ) = PIV( J )
209 PIV( J ) = ITEMP
210 END IF
211 *
212 AJJ = SQRT( AJJ )
213 A( J, J ) = AJJ
214 *
215 * Compute elements J+1:N of row J
216 *
217 IF( J.LT.N ) THEN
218 CALL DGEMV( 'Trans', J-1, N-J, -ONE, A( 1, J+1 ), LDA,
219 $ A( 1, J ), 1, ONE, A( J, J+1 ), LDA )
220 CALL DSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA )
221 END IF
222 *
223 130 CONTINUE
224 *
225 ELSE
226 *
227 * Compute the Cholesky factorization P**T * A * P = L * L**T
228 *
229 DO 150 J = 1, N
230 *
231 * Find pivot, test for exit, else swap rows and columns
232 * Update dot products, compute possible pivots which are
233 * stored in the second half of WORK
234 *
235 DO 140 I = J, N
236 *
237 IF( J.GT.1 ) THEN
238 WORK( I ) = WORK( I ) + A( I, J-1 )**2
239 END IF
240 WORK( N+I ) = A( I, I ) - WORK( I )
241 *
242 140 CONTINUE
243 *
244 IF( J.GT.1 ) THEN
245 ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
246 PVT = ITEMP + J - 1
247 AJJ = WORK( N+PVT )
248 IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN
249 A( J, J ) = AJJ
250 GO TO 160
251 END IF
252 END IF
253 *
254 IF( J.NE.PVT ) THEN
255 *
256 * Pivot OK, so can now swap pivot rows and columns
257 *
258 A( PVT, PVT ) = A( J, J )
259 CALL DSWAP( J-1, A( J, 1 ), LDA, A( PVT, 1 ), LDA )
260 IF( PVT.LT.N )
261 $ CALL DSWAP( N-PVT, A( PVT+1, J ), 1, A( PVT+1, PVT ),
262 $ 1 )
263 CALL DSWAP( PVT-J-1, A( J+1, J ), 1, A( PVT, J+1 ), LDA )
264 *
265 * Swap dot products and PIV
266 *
267 DTEMP = WORK( J )
268 WORK( J ) = WORK( PVT )
269 WORK( PVT ) = DTEMP
270 ITEMP = PIV( PVT )
271 PIV( PVT ) = PIV( J )
272 PIV( J ) = ITEMP
273 END IF
274 *
275 AJJ = SQRT( AJJ )
276 A( J, J ) = AJJ
277 *
278 * Compute elements J+1:N of column J
279 *
280 IF( J.LT.N ) THEN
281 CALL DGEMV( 'No Trans', N-J, J-1, -ONE, A( J+1, 1 ), LDA,
282 $ A( J, 1 ), LDA, ONE, A( J+1, J ), 1 )
283 CALL DSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 )
284 END IF
285 *
286 150 CONTINUE
287 *
288 END IF
289 *
290 * Ran to completion, A has full rank
291 *
292 RANK = N
293 *
294 GO TO 170
295 160 CONTINUE
296 *
297 * Rank is number of steps completed. Set INFO = 1 to signal
298 * that the factorization cannot be used to solve a system.
299 *
300 RANK = J - 1
301 INFO = 1
302 *
303 170 CONTINUE
304 RETURN
305 *
306 * End of DPSTF2
307 *
308 END