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