1 SUBROUTINE ZHETF2( UPLO, N, A, LDA, IPIV, INFO )
2 *
3 * -- LAPACK routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
7 *
8 * .. Scalar Arguments ..
9 CHARACTER UPLO
10 INTEGER INFO, LDA, N
11 * ..
12 * .. Array Arguments ..
13 INTEGER IPIV( * )
14 COMPLEX*16 A( LDA, * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZHETF2 computes the factorization of a complex Hermitian matrix A
21 * using the Bunch-Kaufman diagonal pivoting method:
22 *
23 * A = U*D*U**H or A = L*D*L**H
24 *
25 * where U (or L) is a product of permutation and unit upper (lower)
26 * triangular matrices, U**H is the conjugate transpose of U, and D is
27 * Hermitian and block diagonal with 1-by-1 and 2-by-2 diagonal blocks.
28 *
29 * This is the unblocked version of the algorithm, calling Level 2 BLAS.
30 *
31 * Arguments
32 * =========
33 *
34 * UPLO (input) CHARACTER*1
35 * Specifies whether the upper or lower triangular part of the
36 * Hermitian matrix A is stored:
37 * = 'U': Upper triangular
38 * = 'L': Lower triangular
39 *
40 * N (input) INTEGER
41 * The order of the matrix A. N >= 0.
42 *
43 * A (input/output) COMPLEX*16 array, dimension (LDA,N)
44 * On entry, the Hermitian matrix A. If UPLO = 'U', the leading
45 * n-by-n upper triangular part of A contains the upper
46 * triangular part of the matrix A, and the strictly lower
47 * triangular part of A is not referenced. If UPLO = 'L', the
48 * leading n-by-n lower triangular part of A contains the lower
49 * triangular part of the matrix A, and the strictly upper
50 * triangular part of A is not referenced.
51 *
52 * On exit, the block diagonal matrix D and the multipliers used
53 * to obtain the factor U or L (see below for further details).
54 *
55 * LDA (input) INTEGER
56 * The leading dimension of the array A. LDA >= max(1,N).
57 *
58 * IPIV (output) INTEGER array, dimension (N)
59 * Details of the interchanges and the block structure of D.
60 * If IPIV(k) > 0, then rows and columns k and IPIV(k) were
61 * interchanged and D(k,k) is a 1-by-1 diagonal block.
62 * If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
63 * columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
64 * is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) =
65 * IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
66 * interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
67 *
68 * INFO (output) INTEGER
69 * = 0: successful exit
70 * < 0: if INFO = -k, the k-th argument had an illegal value
71 * > 0: if INFO = k, D(k,k) is exactly zero. The factorization
72 * has been completed, but the block diagonal matrix D is
73 * exactly singular, and division by zero will occur if it
74 * is used to solve a system of equations.
75 *
76 * Further Details
77 * ===============
78 *
79 * 09-29-06 - patch from
80 * Bobby Cheng, MathWorks
81 *
82 * Replace l.210 and l.393
83 * IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
84 * by
85 * IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
86 *
87 * 01-01-96 - Based on modifications by
88 * J. Lewis, Boeing Computer Services Company
89 * A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
90 *
91 * If UPLO = 'U', then A = U*D*U**H, where
92 * U = P(n)*U(n)* ... *P(k)U(k)* ...,
93 * i.e., U is a product of terms P(k)*U(k), where k decreases from n to
94 * 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
95 * and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
96 * defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
97 * that if the diagonal block D(k) is of order s (s = 1 or 2), then
98 *
99 * ( I v 0 ) k-s
100 * U(k) = ( 0 I 0 ) s
101 * ( 0 0 I ) n-k
102 * k-s s n-k
103 *
104 * If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
105 * If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
106 * and A(k,k), and v overwrites A(1:k-2,k-1:k).
107 *
108 * If UPLO = 'L', then A = L*D*L**H, where
109 * L = P(1)*L(1)* ... *P(k)*L(k)* ...,
110 * i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
111 * n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
112 * and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
113 * defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
114 * that if the diagonal block D(k) is of order s (s = 1 or 2), then
115 *
116 * ( I 0 0 ) k-1
117 * L(k) = ( 0 I 0 ) s
118 * ( 0 v I ) n-k-s+1
119 * k-1 s n-k-s+1
120 *
121 * If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
122 * If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
123 * and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
124 *
125 * =====================================================================
126 *
127 * .. Parameters ..
128 DOUBLE PRECISION ZERO, ONE
129 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
130 DOUBLE PRECISION EIGHT, SEVTEN
131 PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
132 * ..
133 * .. Local Scalars ..
134 LOGICAL UPPER
135 INTEGER I, IMAX, J, JMAX, K, KK, KP, KSTEP
136 DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D, D11, D22, R1, ROWMAX,
137 $ TT
138 COMPLEX*16 D12, D21, T, WK, WKM1, WKP1, ZDUM
139 * ..
140 * .. External Functions ..
141 LOGICAL LSAME, DISNAN
142 INTEGER IZAMAX
143 DOUBLE PRECISION DLAPY2
144 EXTERNAL LSAME, IZAMAX, DLAPY2, DISNAN
145 * ..
146 * .. External Subroutines ..
147 EXTERNAL XERBLA, ZDSCAL, ZHER, ZSWAP
148 * ..
149 * .. Intrinsic Functions ..
150 INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, SQRT
151 * ..
152 * .. Statement Functions ..
153 DOUBLE PRECISION CABS1
154 * ..
155 * .. Statement Function definitions ..
156 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
157 * ..
158 * .. Executable Statements ..
159 *
160 * Test the input parameters.
161 *
162 INFO = 0
163 UPPER = LSAME( UPLO, 'U' )
164 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
165 INFO = -1
166 ELSE IF( N.LT.0 ) THEN
167 INFO = -2
168 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
169 INFO = -4
170 END IF
171 IF( INFO.NE.0 ) THEN
172 CALL XERBLA( 'ZHETF2', -INFO )
173 RETURN
174 END IF
175 *
176 * Initialize ALPHA for use in choosing pivot block size.
177 *
178 ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
179 *
180 IF( UPPER ) THEN
181 *
182 * Factorize A as U*D*U**H using the upper triangle of A
183 *
184 * K is the main loop index, decreasing from N to 1 in steps of
185 * 1 or 2
186 *
187 K = N
188 10 CONTINUE
189 *
190 * If K < 1, exit from loop
191 *
192 IF( K.LT.1 )
193 $ GO TO 90
194 KSTEP = 1
195 *
196 * Determine rows and columns to be interchanged and whether
197 * a 1-by-1 or 2-by-2 pivot block will be used
198 *
199 ABSAKK = ABS( DBLE( A( K, K ) ) )
200 *
201 * IMAX is the row-index of the largest off-diagonal element in
202 * column K, and COLMAX is its absolute value
203 *
204 IF( K.GT.1 ) THEN
205 IMAX = IZAMAX( K-1, A( 1, K ), 1 )
206 COLMAX = CABS1( A( IMAX, K ) )
207 ELSE
208 COLMAX = ZERO
209 END IF
210 *
211 IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
212 *
213 * Column K is zero or contains a NaN: set INFO and continue
214 *
215 IF( INFO.EQ.0 )
216 $ INFO = K
217 KP = K
218 A( K, K ) = DBLE( A( K, K ) )
219 ELSE
220 IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
221 *
222 * no interchange, use 1-by-1 pivot block
223 *
224 KP = K
225 ELSE
226 *
227 * JMAX is the column-index of the largest off-diagonal
228 * element in row IMAX, and ROWMAX is its absolute value
229 *
230 JMAX = IMAX + IZAMAX( K-IMAX, A( IMAX, IMAX+1 ), LDA )
231 ROWMAX = CABS1( A( IMAX, JMAX ) )
232 IF( IMAX.GT.1 ) THEN
233 JMAX = IZAMAX( IMAX-1, A( 1, IMAX ), 1 )
234 ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) )
235 END IF
236 *
237 IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
238 *
239 * no interchange, use 1-by-1 pivot block
240 *
241 KP = K
242 ELSE IF( ABS( DBLE( A( IMAX, IMAX ) ) ).GE.ALPHA*ROWMAX )
243 $ THEN
244 *
245 * interchange rows and columns K and IMAX, use 1-by-1
246 * pivot block
247 *
248 KP = IMAX
249 ELSE
250 *
251 * interchange rows and columns K-1 and IMAX, use 2-by-2
252 * pivot block
253 *
254 KP = IMAX
255 KSTEP = 2
256 END IF
257 END IF
258 *
259 KK = K - KSTEP + 1
260 IF( KP.NE.KK ) THEN
261 *
262 * Interchange rows and columns KK and KP in the leading
263 * submatrix A(1:k,1:k)
264 *
265 CALL ZSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
266 DO 20 J = KP + 1, KK - 1
267 T = DCONJG( A( J, KK ) )
268 A( J, KK ) = DCONJG( A( KP, J ) )
269 A( KP, J ) = T
270 20 CONTINUE
271 A( KP, KK ) = DCONJG( A( KP, KK ) )
272 R1 = DBLE( A( KK, KK ) )
273 A( KK, KK ) = DBLE( A( KP, KP ) )
274 A( KP, KP ) = R1
275 IF( KSTEP.EQ.2 ) THEN
276 A( K, K ) = DBLE( A( K, K ) )
277 T = A( K-1, K )
278 A( K-1, K ) = A( KP, K )
279 A( KP, K ) = T
280 END IF
281 ELSE
282 A( K, K ) = DBLE( A( K, K ) )
283 IF( KSTEP.EQ.2 )
284 $ A( K-1, K-1 ) = DBLE( A( K-1, K-1 ) )
285 END IF
286 *
287 * Update the leading submatrix
288 *
289 IF( KSTEP.EQ.1 ) THEN
290 *
291 * 1-by-1 pivot block D(k): column k now holds
292 *
293 * W(k) = U(k)*D(k)
294 *
295 * where U(k) is the k-th column of U
296 *
297 * Perform a rank-1 update of A(1:k-1,1:k-1) as
298 *
299 * A := A - U(k)*D(k)*U(k)**H = A - W(k)*1/D(k)*W(k)**H
300 *
301 R1 = ONE / DBLE( A( K, K ) )
302 CALL ZHER( UPLO, K-1, -R1, A( 1, K ), 1, A, LDA )
303 *
304 * Store U(k) in column k
305 *
306 CALL ZDSCAL( K-1, R1, A( 1, K ), 1 )
307 ELSE
308 *
309 * 2-by-2 pivot block D(k): columns k and k-1 now hold
310 *
311 * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
312 *
313 * where U(k) and U(k-1) are the k-th and (k-1)-th columns
314 * of U
315 *
316 * Perform a rank-2 update of A(1:k-2,1:k-2) as
317 *
318 * A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**H
319 * = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )**H
320 *
321 IF( K.GT.2 ) THEN
322 *
323 D = DLAPY2( DBLE( A( K-1, K ) ),
324 $ DIMAG( A( K-1, K ) ) )
325 D22 = DBLE( A( K-1, K-1 ) ) / D
326 D11 = DBLE( A( K, K ) ) / D
327 TT = ONE / ( D11*D22-ONE )
328 D12 = A( K-1, K ) / D
329 D = TT / D
330 *
331 DO 40 J = K - 2, 1, -1
332 WKM1 = D*( D11*A( J, K-1 )-DCONJG( D12 )*
333 $ A( J, K ) )
334 WK = D*( D22*A( J, K )-D12*A( J, K-1 ) )
335 DO 30 I = J, 1, -1
336 A( I, J ) = A( I, J ) - A( I, K )*DCONJG( WK ) -
337 $ A( I, K-1 )*DCONJG( WKM1 )
338 30 CONTINUE
339 A( J, K ) = WK
340 A( J, K-1 ) = WKM1
341 A( J, J ) = DCMPLX( DBLE( A( J, J ) ), 0.0D+0 )
342 40 CONTINUE
343 *
344 END IF
345 *
346 END IF
347 END IF
348 *
349 * Store details of the interchanges in IPIV
350 *
351 IF( KSTEP.EQ.1 ) THEN
352 IPIV( K ) = KP
353 ELSE
354 IPIV( K ) = -KP
355 IPIV( K-1 ) = -KP
356 END IF
357 *
358 * Decrease K and return to the start of the main loop
359 *
360 K = K - KSTEP
361 GO TO 10
362 *
363 ELSE
364 *
365 * Factorize A as L*D*L**H using the lower triangle of A
366 *
367 * K is the main loop index, increasing from 1 to N in steps of
368 * 1 or 2
369 *
370 K = 1
371 50 CONTINUE
372 *
373 * If K > N, exit from loop
374 *
375 IF( K.GT.N )
376 $ GO TO 90
377 KSTEP = 1
378 *
379 * Determine rows and columns to be interchanged and whether
380 * a 1-by-1 or 2-by-2 pivot block will be used
381 *
382 ABSAKK = ABS( DBLE( A( K, K ) ) )
383 *
384 * IMAX is the row-index of the largest off-diagonal element in
385 * column K, and COLMAX is its absolute value
386 *
387 IF( K.LT.N ) THEN
388 IMAX = K + IZAMAX( N-K, A( K+1, K ), 1 )
389 COLMAX = CABS1( A( IMAX, K ) )
390 ELSE
391 COLMAX = ZERO
392 END IF
393 *
394 IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
395 *
396 * Column K is zero or contains a NaN: set INFO and continue
397 *
398 IF( INFO.EQ.0 )
399 $ INFO = K
400 KP = K
401 A( K, K ) = DBLE( A( K, K ) )
402 ELSE
403 IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
404 *
405 * no interchange, use 1-by-1 pivot block
406 *
407 KP = K
408 ELSE
409 *
410 * JMAX is the column-index of the largest off-diagonal
411 * element in row IMAX, and ROWMAX is its absolute value
412 *
413 JMAX = K - 1 + IZAMAX( IMAX-K, A( IMAX, K ), LDA )
414 ROWMAX = CABS1( A( IMAX, JMAX ) )
415 IF( IMAX.LT.N ) THEN
416 JMAX = IMAX + IZAMAX( N-IMAX, A( IMAX+1, IMAX ), 1 )
417 ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) )
418 END IF
419 *
420 IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
421 *
422 * no interchange, use 1-by-1 pivot block
423 *
424 KP = K
425 ELSE IF( ABS( DBLE( A( IMAX, IMAX ) ) ).GE.ALPHA*ROWMAX )
426 $ THEN
427 *
428 * interchange rows and columns K and IMAX, use 1-by-1
429 * pivot block
430 *
431 KP = IMAX
432 ELSE
433 *
434 * interchange rows and columns K+1 and IMAX, use 2-by-2
435 * pivot block
436 *
437 KP = IMAX
438 KSTEP = 2
439 END IF
440 END IF
441 *
442 KK = K + KSTEP - 1
443 IF( KP.NE.KK ) THEN
444 *
445 * Interchange rows and columns KK and KP in the trailing
446 * submatrix A(k:n,k:n)
447 *
448 IF( KP.LT.N )
449 $ CALL ZSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
450 DO 60 J = KK + 1, KP - 1
451 T = DCONJG( A( J, KK ) )
452 A( J, KK ) = DCONJG( A( KP, J ) )
453 A( KP, J ) = T
454 60 CONTINUE
455 A( KP, KK ) = DCONJG( A( KP, KK ) )
456 R1 = DBLE( A( KK, KK ) )
457 A( KK, KK ) = DBLE( A( KP, KP ) )
458 A( KP, KP ) = R1
459 IF( KSTEP.EQ.2 ) THEN
460 A( K, K ) = DBLE( A( K, K ) )
461 T = A( K+1, K )
462 A( K+1, K ) = A( KP, K )
463 A( KP, K ) = T
464 END IF
465 ELSE
466 A( K, K ) = DBLE( A( K, K ) )
467 IF( KSTEP.EQ.2 )
468 $ A( K+1, K+1 ) = DBLE( A( K+1, K+1 ) )
469 END IF
470 *
471 * Update the trailing submatrix
472 *
473 IF( KSTEP.EQ.1 ) THEN
474 *
475 * 1-by-1 pivot block D(k): column k now holds
476 *
477 * W(k) = L(k)*D(k)
478 *
479 * where L(k) is the k-th column of L
480 *
481 IF( K.LT.N ) THEN
482 *
483 * Perform a rank-1 update of A(k+1:n,k+1:n) as
484 *
485 * A := A - L(k)*D(k)*L(k)**H = A - W(k)*(1/D(k))*W(k)**H
486 *
487 R1 = ONE / DBLE( A( K, K ) )
488 CALL ZHER( UPLO, N-K, -R1, A( K+1, K ), 1,
489 $ A( K+1, K+1 ), LDA )
490 *
491 * Store L(k) in column K
492 *
493 CALL ZDSCAL( N-K, R1, A( K+1, K ), 1 )
494 END IF
495 ELSE
496 *
497 * 2-by-2 pivot block D(k)
498 *
499 IF( K.LT.N-1 ) THEN
500 *
501 * Perform a rank-2 update of A(k+2:n,k+2:n) as
502 *
503 * A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )**H
504 * = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )**H
505 *
506 * where L(k) and L(k+1) are the k-th and (k+1)-th
507 * columns of L
508 *
509 D = DLAPY2( DBLE( A( K+1, K ) ),
510 $ DIMAG( A( K+1, K ) ) )
511 D11 = DBLE( A( K+1, K+1 ) ) / D
512 D22 = DBLE( A( K, K ) ) / D
513 TT = ONE / ( D11*D22-ONE )
514 D21 = A( K+1, K ) / D
515 D = TT / D
516 *
517 DO 80 J = K + 2, N
518 WK = D*( D11*A( J, K )-D21*A( J, K+1 ) )
519 WKP1 = D*( D22*A( J, K+1 )-DCONJG( D21 )*
520 $ A( J, K ) )
521 DO 70 I = J, N
522 A( I, J ) = A( I, J ) - A( I, K )*DCONJG( WK ) -
523 $ A( I, K+1 )*DCONJG( WKP1 )
524 70 CONTINUE
525 A( J, K ) = WK
526 A( J, K+1 ) = WKP1
527 A( J, J ) = DCMPLX( DBLE( A( J, J ) ), 0.0D+0 )
528 80 CONTINUE
529 END IF
530 END IF
531 END IF
532 *
533 * Store details of the interchanges in IPIV
534 *
535 IF( KSTEP.EQ.1 ) THEN
536 IPIV( K ) = KP
537 ELSE
538 IPIV( K ) = -KP
539 IPIV( K+1 ) = -KP
540 END IF
541 *
542 * Increase K and return to the start of the main loop
543 *
544 K = K + KSTEP
545 GO TO 50
546 *
547 END IF
548 *
549 90 CONTINUE
550 RETURN
551 *
552 * End of ZHETF2
553 *
554 END
2 *
3 * -- LAPACK routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
7 *
8 * .. Scalar Arguments ..
9 CHARACTER UPLO
10 INTEGER INFO, LDA, N
11 * ..
12 * .. Array Arguments ..
13 INTEGER IPIV( * )
14 COMPLEX*16 A( LDA, * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZHETF2 computes the factorization of a complex Hermitian matrix A
21 * using the Bunch-Kaufman diagonal pivoting method:
22 *
23 * A = U*D*U**H or A = L*D*L**H
24 *
25 * where U (or L) is a product of permutation and unit upper (lower)
26 * triangular matrices, U**H is the conjugate transpose of U, and D is
27 * Hermitian and block diagonal with 1-by-1 and 2-by-2 diagonal blocks.
28 *
29 * This is the unblocked version of the algorithm, calling Level 2 BLAS.
30 *
31 * Arguments
32 * =========
33 *
34 * UPLO (input) CHARACTER*1
35 * Specifies whether the upper or lower triangular part of the
36 * Hermitian matrix A is stored:
37 * = 'U': Upper triangular
38 * = 'L': Lower triangular
39 *
40 * N (input) INTEGER
41 * The order of the matrix A. N >= 0.
42 *
43 * A (input/output) COMPLEX*16 array, dimension (LDA,N)
44 * On entry, the Hermitian matrix A. If UPLO = 'U', the leading
45 * n-by-n upper triangular part of A contains the upper
46 * triangular part of the matrix A, and the strictly lower
47 * triangular part of A is not referenced. If UPLO = 'L', the
48 * leading n-by-n lower triangular part of A contains the lower
49 * triangular part of the matrix A, and the strictly upper
50 * triangular part of A is not referenced.
51 *
52 * On exit, the block diagonal matrix D and the multipliers used
53 * to obtain the factor U or L (see below for further details).
54 *
55 * LDA (input) INTEGER
56 * The leading dimension of the array A. LDA >= max(1,N).
57 *
58 * IPIV (output) INTEGER array, dimension (N)
59 * Details of the interchanges and the block structure of D.
60 * If IPIV(k) > 0, then rows and columns k and IPIV(k) were
61 * interchanged and D(k,k) is a 1-by-1 diagonal block.
62 * If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
63 * columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
64 * is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) =
65 * IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
66 * interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
67 *
68 * INFO (output) INTEGER
69 * = 0: successful exit
70 * < 0: if INFO = -k, the k-th argument had an illegal value
71 * > 0: if INFO = k, D(k,k) is exactly zero. The factorization
72 * has been completed, but the block diagonal matrix D is
73 * exactly singular, and division by zero will occur if it
74 * is used to solve a system of equations.
75 *
76 * Further Details
77 * ===============
78 *
79 * 09-29-06 - patch from
80 * Bobby Cheng, MathWorks
81 *
82 * Replace l.210 and l.393
83 * IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
84 * by
85 * IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
86 *
87 * 01-01-96 - Based on modifications by
88 * J. Lewis, Boeing Computer Services Company
89 * A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
90 *
91 * If UPLO = 'U', then A = U*D*U**H, where
92 * U = P(n)*U(n)* ... *P(k)U(k)* ...,
93 * i.e., U is a product of terms P(k)*U(k), where k decreases from n to
94 * 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
95 * and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
96 * defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
97 * that if the diagonal block D(k) is of order s (s = 1 or 2), then
98 *
99 * ( I v 0 ) k-s
100 * U(k) = ( 0 I 0 ) s
101 * ( 0 0 I ) n-k
102 * k-s s n-k
103 *
104 * If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
105 * If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
106 * and A(k,k), and v overwrites A(1:k-2,k-1:k).
107 *
108 * If UPLO = 'L', then A = L*D*L**H, where
109 * L = P(1)*L(1)* ... *P(k)*L(k)* ...,
110 * i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
111 * n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
112 * and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
113 * defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
114 * that if the diagonal block D(k) is of order s (s = 1 or 2), then
115 *
116 * ( I 0 0 ) k-1
117 * L(k) = ( 0 I 0 ) s
118 * ( 0 v I ) n-k-s+1
119 * k-1 s n-k-s+1
120 *
121 * If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
122 * If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
123 * and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
124 *
125 * =====================================================================
126 *
127 * .. Parameters ..
128 DOUBLE PRECISION ZERO, ONE
129 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
130 DOUBLE PRECISION EIGHT, SEVTEN
131 PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
132 * ..
133 * .. Local Scalars ..
134 LOGICAL UPPER
135 INTEGER I, IMAX, J, JMAX, K, KK, KP, KSTEP
136 DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D, D11, D22, R1, ROWMAX,
137 $ TT
138 COMPLEX*16 D12, D21, T, WK, WKM1, WKP1, ZDUM
139 * ..
140 * .. External Functions ..
141 LOGICAL LSAME, DISNAN
142 INTEGER IZAMAX
143 DOUBLE PRECISION DLAPY2
144 EXTERNAL LSAME, IZAMAX, DLAPY2, DISNAN
145 * ..
146 * .. External Subroutines ..
147 EXTERNAL XERBLA, ZDSCAL, ZHER, ZSWAP
148 * ..
149 * .. Intrinsic Functions ..
150 INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, SQRT
151 * ..
152 * .. Statement Functions ..
153 DOUBLE PRECISION CABS1
154 * ..
155 * .. Statement Function definitions ..
156 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
157 * ..
158 * .. Executable Statements ..
159 *
160 * Test the input parameters.
161 *
162 INFO = 0
163 UPPER = LSAME( UPLO, 'U' )
164 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
165 INFO = -1
166 ELSE IF( N.LT.0 ) THEN
167 INFO = -2
168 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
169 INFO = -4
170 END IF
171 IF( INFO.NE.0 ) THEN
172 CALL XERBLA( 'ZHETF2', -INFO )
173 RETURN
174 END IF
175 *
176 * Initialize ALPHA for use in choosing pivot block size.
177 *
178 ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
179 *
180 IF( UPPER ) THEN
181 *
182 * Factorize A as U*D*U**H using the upper triangle of A
183 *
184 * K is the main loop index, decreasing from N to 1 in steps of
185 * 1 or 2
186 *
187 K = N
188 10 CONTINUE
189 *
190 * If K < 1, exit from loop
191 *
192 IF( K.LT.1 )
193 $ GO TO 90
194 KSTEP = 1
195 *
196 * Determine rows and columns to be interchanged and whether
197 * a 1-by-1 or 2-by-2 pivot block will be used
198 *
199 ABSAKK = ABS( DBLE( A( K, K ) ) )
200 *
201 * IMAX is the row-index of the largest off-diagonal element in
202 * column K, and COLMAX is its absolute value
203 *
204 IF( K.GT.1 ) THEN
205 IMAX = IZAMAX( K-1, A( 1, K ), 1 )
206 COLMAX = CABS1( A( IMAX, K ) )
207 ELSE
208 COLMAX = ZERO
209 END IF
210 *
211 IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
212 *
213 * Column K is zero or contains a NaN: set INFO and continue
214 *
215 IF( INFO.EQ.0 )
216 $ INFO = K
217 KP = K
218 A( K, K ) = DBLE( A( K, K ) )
219 ELSE
220 IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
221 *
222 * no interchange, use 1-by-1 pivot block
223 *
224 KP = K
225 ELSE
226 *
227 * JMAX is the column-index of the largest off-diagonal
228 * element in row IMAX, and ROWMAX is its absolute value
229 *
230 JMAX = IMAX + IZAMAX( K-IMAX, A( IMAX, IMAX+1 ), LDA )
231 ROWMAX = CABS1( A( IMAX, JMAX ) )
232 IF( IMAX.GT.1 ) THEN
233 JMAX = IZAMAX( IMAX-1, A( 1, IMAX ), 1 )
234 ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) )
235 END IF
236 *
237 IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
238 *
239 * no interchange, use 1-by-1 pivot block
240 *
241 KP = K
242 ELSE IF( ABS( DBLE( A( IMAX, IMAX ) ) ).GE.ALPHA*ROWMAX )
243 $ THEN
244 *
245 * interchange rows and columns K and IMAX, use 1-by-1
246 * pivot block
247 *
248 KP = IMAX
249 ELSE
250 *
251 * interchange rows and columns K-1 and IMAX, use 2-by-2
252 * pivot block
253 *
254 KP = IMAX
255 KSTEP = 2
256 END IF
257 END IF
258 *
259 KK = K - KSTEP + 1
260 IF( KP.NE.KK ) THEN
261 *
262 * Interchange rows and columns KK and KP in the leading
263 * submatrix A(1:k,1:k)
264 *
265 CALL ZSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
266 DO 20 J = KP + 1, KK - 1
267 T = DCONJG( A( J, KK ) )
268 A( J, KK ) = DCONJG( A( KP, J ) )
269 A( KP, J ) = T
270 20 CONTINUE
271 A( KP, KK ) = DCONJG( A( KP, KK ) )
272 R1 = DBLE( A( KK, KK ) )
273 A( KK, KK ) = DBLE( A( KP, KP ) )
274 A( KP, KP ) = R1
275 IF( KSTEP.EQ.2 ) THEN
276 A( K, K ) = DBLE( A( K, K ) )
277 T = A( K-1, K )
278 A( K-1, K ) = A( KP, K )
279 A( KP, K ) = T
280 END IF
281 ELSE
282 A( K, K ) = DBLE( A( K, K ) )
283 IF( KSTEP.EQ.2 )
284 $ A( K-1, K-1 ) = DBLE( A( K-1, K-1 ) )
285 END IF
286 *
287 * Update the leading submatrix
288 *
289 IF( KSTEP.EQ.1 ) THEN
290 *
291 * 1-by-1 pivot block D(k): column k now holds
292 *
293 * W(k) = U(k)*D(k)
294 *
295 * where U(k) is the k-th column of U
296 *
297 * Perform a rank-1 update of A(1:k-1,1:k-1) as
298 *
299 * A := A - U(k)*D(k)*U(k)**H = A - W(k)*1/D(k)*W(k)**H
300 *
301 R1 = ONE / DBLE( A( K, K ) )
302 CALL ZHER( UPLO, K-1, -R1, A( 1, K ), 1, A, LDA )
303 *
304 * Store U(k) in column k
305 *
306 CALL ZDSCAL( K-1, R1, A( 1, K ), 1 )
307 ELSE
308 *
309 * 2-by-2 pivot block D(k): columns k and k-1 now hold
310 *
311 * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
312 *
313 * where U(k) and U(k-1) are the k-th and (k-1)-th columns
314 * of U
315 *
316 * Perform a rank-2 update of A(1:k-2,1:k-2) as
317 *
318 * A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**H
319 * = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )**H
320 *
321 IF( K.GT.2 ) THEN
322 *
323 D = DLAPY2( DBLE( A( K-1, K ) ),
324 $ DIMAG( A( K-1, K ) ) )
325 D22 = DBLE( A( K-1, K-1 ) ) / D
326 D11 = DBLE( A( K, K ) ) / D
327 TT = ONE / ( D11*D22-ONE )
328 D12 = A( K-1, K ) / D
329 D = TT / D
330 *
331 DO 40 J = K - 2, 1, -1
332 WKM1 = D*( D11*A( J, K-1 )-DCONJG( D12 )*
333 $ A( J, K ) )
334 WK = D*( D22*A( J, K )-D12*A( J, K-1 ) )
335 DO 30 I = J, 1, -1
336 A( I, J ) = A( I, J ) - A( I, K )*DCONJG( WK ) -
337 $ A( I, K-1 )*DCONJG( WKM1 )
338 30 CONTINUE
339 A( J, K ) = WK
340 A( J, K-1 ) = WKM1
341 A( J, J ) = DCMPLX( DBLE( A( J, J ) ), 0.0D+0 )
342 40 CONTINUE
343 *
344 END IF
345 *
346 END IF
347 END IF
348 *
349 * Store details of the interchanges in IPIV
350 *
351 IF( KSTEP.EQ.1 ) THEN
352 IPIV( K ) = KP
353 ELSE
354 IPIV( K ) = -KP
355 IPIV( K-1 ) = -KP
356 END IF
357 *
358 * Decrease K and return to the start of the main loop
359 *
360 K = K - KSTEP
361 GO TO 10
362 *
363 ELSE
364 *
365 * Factorize A as L*D*L**H using the lower triangle of A
366 *
367 * K is the main loop index, increasing from 1 to N in steps of
368 * 1 or 2
369 *
370 K = 1
371 50 CONTINUE
372 *
373 * If K > N, exit from loop
374 *
375 IF( K.GT.N )
376 $ GO TO 90
377 KSTEP = 1
378 *
379 * Determine rows and columns to be interchanged and whether
380 * a 1-by-1 or 2-by-2 pivot block will be used
381 *
382 ABSAKK = ABS( DBLE( A( K, K ) ) )
383 *
384 * IMAX is the row-index of the largest off-diagonal element in
385 * column K, and COLMAX is its absolute value
386 *
387 IF( K.LT.N ) THEN
388 IMAX = K + IZAMAX( N-K, A( K+1, K ), 1 )
389 COLMAX = CABS1( A( IMAX, K ) )
390 ELSE
391 COLMAX = ZERO
392 END IF
393 *
394 IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
395 *
396 * Column K is zero or contains a NaN: set INFO and continue
397 *
398 IF( INFO.EQ.0 )
399 $ INFO = K
400 KP = K
401 A( K, K ) = DBLE( A( K, K ) )
402 ELSE
403 IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
404 *
405 * no interchange, use 1-by-1 pivot block
406 *
407 KP = K
408 ELSE
409 *
410 * JMAX is the column-index of the largest off-diagonal
411 * element in row IMAX, and ROWMAX is its absolute value
412 *
413 JMAX = K - 1 + IZAMAX( IMAX-K, A( IMAX, K ), LDA )
414 ROWMAX = CABS1( A( IMAX, JMAX ) )
415 IF( IMAX.LT.N ) THEN
416 JMAX = IMAX + IZAMAX( N-IMAX, A( IMAX+1, IMAX ), 1 )
417 ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) )
418 END IF
419 *
420 IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
421 *
422 * no interchange, use 1-by-1 pivot block
423 *
424 KP = K
425 ELSE IF( ABS( DBLE( A( IMAX, IMAX ) ) ).GE.ALPHA*ROWMAX )
426 $ THEN
427 *
428 * interchange rows and columns K and IMAX, use 1-by-1
429 * pivot block
430 *
431 KP = IMAX
432 ELSE
433 *
434 * interchange rows and columns K+1 and IMAX, use 2-by-2
435 * pivot block
436 *
437 KP = IMAX
438 KSTEP = 2
439 END IF
440 END IF
441 *
442 KK = K + KSTEP - 1
443 IF( KP.NE.KK ) THEN
444 *
445 * Interchange rows and columns KK and KP in the trailing
446 * submatrix A(k:n,k:n)
447 *
448 IF( KP.LT.N )
449 $ CALL ZSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
450 DO 60 J = KK + 1, KP - 1
451 T = DCONJG( A( J, KK ) )
452 A( J, KK ) = DCONJG( A( KP, J ) )
453 A( KP, J ) = T
454 60 CONTINUE
455 A( KP, KK ) = DCONJG( A( KP, KK ) )
456 R1 = DBLE( A( KK, KK ) )
457 A( KK, KK ) = DBLE( A( KP, KP ) )
458 A( KP, KP ) = R1
459 IF( KSTEP.EQ.2 ) THEN
460 A( K, K ) = DBLE( A( K, K ) )
461 T = A( K+1, K )
462 A( K+1, K ) = A( KP, K )
463 A( KP, K ) = T
464 END IF
465 ELSE
466 A( K, K ) = DBLE( A( K, K ) )
467 IF( KSTEP.EQ.2 )
468 $ A( K+1, K+1 ) = DBLE( A( K+1, K+1 ) )
469 END IF
470 *
471 * Update the trailing submatrix
472 *
473 IF( KSTEP.EQ.1 ) THEN
474 *
475 * 1-by-1 pivot block D(k): column k now holds
476 *
477 * W(k) = L(k)*D(k)
478 *
479 * where L(k) is the k-th column of L
480 *
481 IF( K.LT.N ) THEN
482 *
483 * Perform a rank-1 update of A(k+1:n,k+1:n) as
484 *
485 * A := A - L(k)*D(k)*L(k)**H = A - W(k)*(1/D(k))*W(k)**H
486 *
487 R1 = ONE / DBLE( A( K, K ) )
488 CALL ZHER( UPLO, N-K, -R1, A( K+1, K ), 1,
489 $ A( K+1, K+1 ), LDA )
490 *
491 * Store L(k) in column K
492 *
493 CALL ZDSCAL( N-K, R1, A( K+1, K ), 1 )
494 END IF
495 ELSE
496 *
497 * 2-by-2 pivot block D(k)
498 *
499 IF( K.LT.N-1 ) THEN
500 *
501 * Perform a rank-2 update of A(k+2:n,k+2:n) as
502 *
503 * A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )**H
504 * = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )**H
505 *
506 * where L(k) and L(k+1) are the k-th and (k+1)-th
507 * columns of L
508 *
509 D = DLAPY2( DBLE( A( K+1, K ) ),
510 $ DIMAG( A( K+1, K ) ) )
511 D11 = DBLE( A( K+1, K+1 ) ) / D
512 D22 = DBLE( A( K, K ) ) / D
513 TT = ONE / ( D11*D22-ONE )
514 D21 = A( K+1, K ) / D
515 D = TT / D
516 *
517 DO 80 J = K + 2, N
518 WK = D*( D11*A( J, K )-D21*A( J, K+1 ) )
519 WKP1 = D*( D22*A( J, K+1 )-DCONJG( D21 )*
520 $ A( J, K ) )
521 DO 70 I = J, N
522 A( I, J ) = A( I, J ) - A( I, K )*DCONJG( WK ) -
523 $ A( I, K+1 )*DCONJG( WKP1 )
524 70 CONTINUE
525 A( J, K ) = WK
526 A( J, K+1 ) = WKP1
527 A( J, J ) = DCMPLX( DBLE( A( J, J ) ), 0.0D+0 )
528 80 CONTINUE
529 END IF
530 END IF
531 END IF
532 *
533 * Store details of the interchanges in IPIV
534 *
535 IF( KSTEP.EQ.1 ) THEN
536 IPIV( K ) = KP
537 ELSE
538 IPIV( K ) = -KP
539 IPIV( K+1 ) = -KP
540 END IF
541 *
542 * Increase K and return to the start of the main loop
543 *
544 K = K + KSTEP
545 GO TO 50
546 *
547 END IF
548 *
549 90 CONTINUE
550 RETURN
551 *
552 * End of ZHETF2
553 *
554 END