1 SUBROUTINE DPSTRF( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
2 *
3 * -- LAPACK 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 * DPSTRF 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 3 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 * LDA (input) INTEGER
57 * The leading dimension of the array A. LDA >= max(1,N).
58 *
59 * PIV (output) INTEGER array, dimension (N)
60 * PIV is such that the nonzero entries are P( PIV(K), K ) = 1.
61 *
62 * RANK (output) INTEGER
63 * The rank of A given by the number of steps the algorithm
64 * completed.
65 *
66 * TOL (input) DOUBLE PRECISION
67 * User defined tolerance. If TOL < 0, then N*U*MAX( A(K,K) )
68 * will be used. The algorithm terminates at the (K-1)st step
69 * if the pivot <= TOL.
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, JB, K, NB, PVT
90 LOGICAL UPPER
91 * ..
92 * .. External Functions ..
93 DOUBLE PRECISION DLAMCH
94 INTEGER ILAENV
95 LOGICAL LSAME, DISNAN
96 EXTERNAL DLAMCH, ILAENV, LSAME, DISNAN
97 * ..
98 * .. External Subroutines ..
99 EXTERNAL DGEMV, DPSTF2, DSCAL, DSWAP, DSYRK, XERBLA
100 * ..
101 * .. Intrinsic Functions ..
102 INTRINSIC MAX, MIN, SQRT, MAXLOC
103 * ..
104 * .. Executable Statements ..
105 *
106 * Test the input parameters.
107 *
108 INFO = 0
109 UPPER = LSAME( UPLO, 'U' )
110 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
111 INFO = -1
112 ELSE IF( N.LT.0 ) THEN
113 INFO = -2
114 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
115 INFO = -4
116 END IF
117 IF( INFO.NE.0 ) THEN
118 CALL XERBLA( 'DPSTRF', -INFO )
119 RETURN
120 END IF
121 *
122 * Quick return if possible
123 *
124 IF( N.EQ.0 )
125 $ RETURN
126 *
127 * Get block size
128 *
129 NB = ILAENV( 1, 'DPOTRF', UPLO, N, -1, -1, -1 )
130 IF( NB.LE.1 .OR. NB.GE.N ) THEN
131 *
132 * Use unblocked code
133 *
134 CALL DPSTF2( UPLO, N, A( 1, 1 ), LDA, PIV, RANK, TOL, WORK,
135 $ INFO )
136 GO TO 200
137 *
138 ELSE
139 *
140 * Initialize PIV
141 *
142 DO 100 I = 1, N
143 PIV( I ) = I
144 100 CONTINUE
145 *
146 * Compute stopping value
147 *
148 PVT = 1
149 AJJ = A( PVT, PVT )
150 DO I = 2, N
151 IF( A( I, I ).GT.AJJ ) THEN
152 PVT = I
153 AJJ = A( PVT, PVT )
154 END IF
155 END DO
156 IF( AJJ.EQ.ZERO.OR.DISNAN( AJJ ) ) THEN
157 RANK = 0
158 INFO = 1
159 GO TO 200
160 END IF
161 *
162 * Compute stopping value if not supplied
163 *
164 IF( TOL.LT.ZERO ) THEN
165 DSTOP = N * DLAMCH( 'Epsilon' ) * AJJ
166 ELSE
167 DSTOP = TOL
168 END IF
169 *
170 *
171 IF( UPPER ) THEN
172 *
173 * Compute the Cholesky factorization P**T * A * P = U**T * U
174 *
175 DO 140 K = 1, N, NB
176 *
177 * Account for last block not being NB wide
178 *
179 JB = MIN( NB, N-K+1 )
180 *
181 * Set relevant part of first half of WORK to zero,
182 * holds dot products
183 *
184 DO 110 I = K, N
185 WORK( I ) = 0
186 110 CONTINUE
187 *
188 DO 130 J = K, K + JB - 1
189 *
190 * Find pivot, test for exit, else swap rows and columns
191 * Update dot products, compute possible pivots which are
192 * stored in the second half of WORK
193 *
194 DO 120 I = J, N
195 *
196 IF( J.GT.K ) THEN
197 WORK( I ) = WORK( I ) + A( J-1, I )**2
198 END IF
199 WORK( N+I ) = A( I, I ) - WORK( I )
200 *
201 120 CONTINUE
202 *
203 IF( J.GT.1 ) THEN
204 ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
205 PVT = ITEMP + J - 1
206 AJJ = WORK( N+PVT )
207 IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN
208 A( J, J ) = AJJ
209 GO TO 190
210 END IF
211 END IF
212 *
213 IF( J.NE.PVT ) THEN
214 *
215 * Pivot OK, so can now swap pivot rows and columns
216 *
217 A( PVT, PVT ) = A( J, J )
218 CALL DSWAP( J-1, A( 1, J ), 1, A( 1, PVT ), 1 )
219 IF( PVT.LT.N )
220 $ CALL DSWAP( N-PVT, A( J, PVT+1 ), LDA,
221 $ A( PVT, PVT+1 ), LDA )
222 CALL DSWAP( PVT-J-1, A( J, J+1 ), LDA,
223 $ A( J+1, PVT ), 1 )
224 *
225 * Swap dot products and PIV
226 *
227 DTEMP = WORK( J )
228 WORK( J ) = WORK( PVT )
229 WORK( PVT ) = DTEMP
230 ITEMP = PIV( PVT )
231 PIV( PVT ) = PIV( J )
232 PIV( J ) = ITEMP
233 END IF
234 *
235 AJJ = SQRT( AJJ )
236 A( J, J ) = AJJ
237 *
238 * Compute elements J+1:N of row J.
239 *
240 IF( J.LT.N ) THEN
241 CALL DGEMV( 'Trans', J-K, N-J, -ONE, A( K, J+1 ),
242 $ LDA, A( K, J ), 1, ONE, A( J, J+1 ),
243 $ LDA )
244 CALL DSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA )
245 END IF
246 *
247 130 CONTINUE
248 *
249 * Update trailing matrix, J already incremented
250 *
251 IF( K+JB.LE.N ) THEN
252 CALL DSYRK( 'Upper', 'Trans', N-J+1, JB, -ONE,
253 $ A( K, J ), LDA, ONE, A( J, J ), LDA )
254 END IF
255 *
256 140 CONTINUE
257 *
258 ELSE
259 *
260 * Compute the Cholesky factorization P**T * A * P = L * L**T
261 *
262 DO 180 K = 1, N, NB
263 *
264 * Account for last block not being NB wide
265 *
266 JB = MIN( NB, N-K+1 )
267 *
268 * Set relevant part of first half of WORK to zero,
269 * holds dot products
270 *
271 DO 150 I = K, N
272 WORK( I ) = 0
273 150 CONTINUE
274 *
275 DO 170 J = K, K + JB - 1
276 *
277 * Find pivot, test for exit, else swap rows and columns
278 * Update dot products, compute possible pivots which are
279 * stored in the second half of WORK
280 *
281 DO 160 I = J, N
282 *
283 IF( J.GT.K ) THEN
284 WORK( I ) = WORK( I ) + A( I, J-1 )**2
285 END IF
286 WORK( N+I ) = A( I, I ) - WORK( I )
287 *
288 160 CONTINUE
289 *
290 IF( J.GT.1 ) THEN
291 ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
292 PVT = ITEMP + J - 1
293 AJJ = WORK( N+PVT )
294 IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN
295 A( J, J ) = AJJ
296 GO TO 190
297 END IF
298 END IF
299 *
300 IF( J.NE.PVT ) THEN
301 *
302 * Pivot OK, so can now swap pivot rows and columns
303 *
304 A( PVT, PVT ) = A( J, J )
305 CALL DSWAP( J-1, A( J, 1 ), LDA, A( PVT, 1 ), LDA )
306 IF( PVT.LT.N )
307 $ CALL DSWAP( N-PVT, A( PVT+1, J ), 1,
308 $ A( PVT+1, PVT ), 1 )
309 CALL DSWAP( PVT-J-1, A( J+1, J ), 1, A( PVT, J+1 ),
310 $ LDA )
311 *
312 * Swap dot products and PIV
313 *
314 DTEMP = WORK( J )
315 WORK( J ) = WORK( PVT )
316 WORK( PVT ) = DTEMP
317 ITEMP = PIV( PVT )
318 PIV( PVT ) = PIV( J )
319 PIV( J ) = ITEMP
320 END IF
321 *
322 AJJ = SQRT( AJJ )
323 A( J, J ) = AJJ
324 *
325 * Compute elements J+1:N of column J.
326 *
327 IF( J.LT.N ) THEN
328 CALL DGEMV( 'No Trans', N-J, J-K, -ONE,
329 $ A( J+1, K ), LDA, A( J, K ), LDA, ONE,
330 $ A( J+1, J ), 1 )
331 CALL DSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 )
332 END IF
333 *
334 170 CONTINUE
335 *
336 * Update trailing matrix, J already incremented
337 *
338 IF( K+JB.LE.N ) THEN
339 CALL DSYRK( 'Lower', 'No Trans', N-J+1, JB, -ONE,
340 $ A( J, K ), LDA, ONE, A( J, J ), LDA )
341 END IF
342 *
343 180 CONTINUE
344 *
345 END IF
346 END IF
347 *
348 * Ran to completion, A has full rank
349 *
350 RANK = N
351 *
352 GO TO 200
353 190 CONTINUE
354 *
355 * Rank is the number of steps completed. Set INFO = 1 to signal
356 * that the factorization cannot be used to solve a system.
357 *
358 RANK = J - 1
359 INFO = 1
360 *
361 200 CONTINUE
362 RETURN
363 *
364 * End of DPSTRF
365 *
366 END
2 *
3 * -- LAPACK 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 * DPSTRF 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 3 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 * LDA (input) INTEGER
57 * The leading dimension of the array A. LDA >= max(1,N).
58 *
59 * PIV (output) INTEGER array, dimension (N)
60 * PIV is such that the nonzero entries are P( PIV(K), K ) = 1.
61 *
62 * RANK (output) INTEGER
63 * The rank of A given by the number of steps the algorithm
64 * completed.
65 *
66 * TOL (input) DOUBLE PRECISION
67 * User defined tolerance. If TOL < 0, then N*U*MAX( A(K,K) )
68 * will be used. The algorithm terminates at the (K-1)st step
69 * if the pivot <= TOL.
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, JB, K, NB, PVT
90 LOGICAL UPPER
91 * ..
92 * .. External Functions ..
93 DOUBLE PRECISION DLAMCH
94 INTEGER ILAENV
95 LOGICAL LSAME, DISNAN
96 EXTERNAL DLAMCH, ILAENV, LSAME, DISNAN
97 * ..
98 * .. External Subroutines ..
99 EXTERNAL DGEMV, DPSTF2, DSCAL, DSWAP, DSYRK, XERBLA
100 * ..
101 * .. Intrinsic Functions ..
102 INTRINSIC MAX, MIN, SQRT, MAXLOC
103 * ..
104 * .. Executable Statements ..
105 *
106 * Test the input parameters.
107 *
108 INFO = 0
109 UPPER = LSAME( UPLO, 'U' )
110 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
111 INFO = -1
112 ELSE IF( N.LT.0 ) THEN
113 INFO = -2
114 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
115 INFO = -4
116 END IF
117 IF( INFO.NE.0 ) THEN
118 CALL XERBLA( 'DPSTRF', -INFO )
119 RETURN
120 END IF
121 *
122 * Quick return if possible
123 *
124 IF( N.EQ.0 )
125 $ RETURN
126 *
127 * Get block size
128 *
129 NB = ILAENV( 1, 'DPOTRF', UPLO, N, -1, -1, -1 )
130 IF( NB.LE.1 .OR. NB.GE.N ) THEN
131 *
132 * Use unblocked code
133 *
134 CALL DPSTF2( UPLO, N, A( 1, 1 ), LDA, PIV, RANK, TOL, WORK,
135 $ INFO )
136 GO TO 200
137 *
138 ELSE
139 *
140 * Initialize PIV
141 *
142 DO 100 I = 1, N
143 PIV( I ) = I
144 100 CONTINUE
145 *
146 * Compute stopping value
147 *
148 PVT = 1
149 AJJ = A( PVT, PVT )
150 DO I = 2, N
151 IF( A( I, I ).GT.AJJ ) THEN
152 PVT = I
153 AJJ = A( PVT, PVT )
154 END IF
155 END DO
156 IF( AJJ.EQ.ZERO.OR.DISNAN( AJJ ) ) THEN
157 RANK = 0
158 INFO = 1
159 GO TO 200
160 END IF
161 *
162 * Compute stopping value if not supplied
163 *
164 IF( TOL.LT.ZERO ) THEN
165 DSTOP = N * DLAMCH( 'Epsilon' ) * AJJ
166 ELSE
167 DSTOP = TOL
168 END IF
169 *
170 *
171 IF( UPPER ) THEN
172 *
173 * Compute the Cholesky factorization P**T * A * P = U**T * U
174 *
175 DO 140 K = 1, N, NB
176 *
177 * Account for last block not being NB wide
178 *
179 JB = MIN( NB, N-K+1 )
180 *
181 * Set relevant part of first half of WORK to zero,
182 * holds dot products
183 *
184 DO 110 I = K, N
185 WORK( I ) = 0
186 110 CONTINUE
187 *
188 DO 130 J = K, K + JB - 1
189 *
190 * Find pivot, test for exit, else swap rows and columns
191 * Update dot products, compute possible pivots which are
192 * stored in the second half of WORK
193 *
194 DO 120 I = J, N
195 *
196 IF( J.GT.K ) THEN
197 WORK( I ) = WORK( I ) + A( J-1, I )**2
198 END IF
199 WORK( N+I ) = A( I, I ) - WORK( I )
200 *
201 120 CONTINUE
202 *
203 IF( J.GT.1 ) THEN
204 ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
205 PVT = ITEMP + J - 1
206 AJJ = WORK( N+PVT )
207 IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN
208 A( J, J ) = AJJ
209 GO TO 190
210 END IF
211 END IF
212 *
213 IF( J.NE.PVT ) THEN
214 *
215 * Pivot OK, so can now swap pivot rows and columns
216 *
217 A( PVT, PVT ) = A( J, J )
218 CALL DSWAP( J-1, A( 1, J ), 1, A( 1, PVT ), 1 )
219 IF( PVT.LT.N )
220 $ CALL DSWAP( N-PVT, A( J, PVT+1 ), LDA,
221 $ A( PVT, PVT+1 ), LDA )
222 CALL DSWAP( PVT-J-1, A( J, J+1 ), LDA,
223 $ A( J+1, PVT ), 1 )
224 *
225 * Swap dot products and PIV
226 *
227 DTEMP = WORK( J )
228 WORK( J ) = WORK( PVT )
229 WORK( PVT ) = DTEMP
230 ITEMP = PIV( PVT )
231 PIV( PVT ) = PIV( J )
232 PIV( J ) = ITEMP
233 END IF
234 *
235 AJJ = SQRT( AJJ )
236 A( J, J ) = AJJ
237 *
238 * Compute elements J+1:N of row J.
239 *
240 IF( J.LT.N ) THEN
241 CALL DGEMV( 'Trans', J-K, N-J, -ONE, A( K, J+1 ),
242 $ LDA, A( K, J ), 1, ONE, A( J, J+1 ),
243 $ LDA )
244 CALL DSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA )
245 END IF
246 *
247 130 CONTINUE
248 *
249 * Update trailing matrix, J already incremented
250 *
251 IF( K+JB.LE.N ) THEN
252 CALL DSYRK( 'Upper', 'Trans', N-J+1, JB, -ONE,
253 $ A( K, J ), LDA, ONE, A( J, J ), LDA )
254 END IF
255 *
256 140 CONTINUE
257 *
258 ELSE
259 *
260 * Compute the Cholesky factorization P**T * A * P = L * L**T
261 *
262 DO 180 K = 1, N, NB
263 *
264 * Account for last block not being NB wide
265 *
266 JB = MIN( NB, N-K+1 )
267 *
268 * Set relevant part of first half of WORK to zero,
269 * holds dot products
270 *
271 DO 150 I = K, N
272 WORK( I ) = 0
273 150 CONTINUE
274 *
275 DO 170 J = K, K + JB - 1
276 *
277 * Find pivot, test for exit, else swap rows and columns
278 * Update dot products, compute possible pivots which are
279 * stored in the second half of WORK
280 *
281 DO 160 I = J, N
282 *
283 IF( J.GT.K ) THEN
284 WORK( I ) = WORK( I ) + A( I, J-1 )**2
285 END IF
286 WORK( N+I ) = A( I, I ) - WORK( I )
287 *
288 160 CONTINUE
289 *
290 IF( J.GT.1 ) THEN
291 ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
292 PVT = ITEMP + J - 1
293 AJJ = WORK( N+PVT )
294 IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN
295 A( J, J ) = AJJ
296 GO TO 190
297 END IF
298 END IF
299 *
300 IF( J.NE.PVT ) THEN
301 *
302 * Pivot OK, so can now swap pivot rows and columns
303 *
304 A( PVT, PVT ) = A( J, J )
305 CALL DSWAP( J-1, A( J, 1 ), LDA, A( PVT, 1 ), LDA )
306 IF( PVT.LT.N )
307 $ CALL DSWAP( N-PVT, A( PVT+1, J ), 1,
308 $ A( PVT+1, PVT ), 1 )
309 CALL DSWAP( PVT-J-1, A( J+1, J ), 1, A( PVT, J+1 ),
310 $ LDA )
311 *
312 * Swap dot products and PIV
313 *
314 DTEMP = WORK( J )
315 WORK( J ) = WORK( PVT )
316 WORK( PVT ) = DTEMP
317 ITEMP = PIV( PVT )
318 PIV( PVT ) = PIV( J )
319 PIV( J ) = ITEMP
320 END IF
321 *
322 AJJ = SQRT( AJJ )
323 A( J, J ) = AJJ
324 *
325 * Compute elements J+1:N of column J.
326 *
327 IF( J.LT.N ) THEN
328 CALL DGEMV( 'No Trans', N-J, J-K, -ONE,
329 $ A( J+1, K ), LDA, A( J, K ), LDA, ONE,
330 $ A( J+1, J ), 1 )
331 CALL DSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 )
332 END IF
333 *
334 170 CONTINUE
335 *
336 * Update trailing matrix, J already incremented
337 *
338 IF( K+JB.LE.N ) THEN
339 CALL DSYRK( 'Lower', 'No Trans', N-J+1, JB, -ONE,
340 $ A( J, K ), LDA, ONE, A( J, J ), LDA )
341 END IF
342 *
343 180 CONTINUE
344 *
345 END IF
346 END IF
347 *
348 * Ran to completion, A has full rank
349 *
350 RANK = N
351 *
352 GO TO 200
353 190 CONTINUE
354 *
355 * Rank is the number of steps completed. Set INFO = 1 to signal
356 * that the factorization cannot be used to solve a system.
357 *
358 RANK = J - 1
359 INFO = 1
360 *
361 200 CONTINUE
362 RETURN
363 *
364 * End of DPSTRF
365 *
366 END