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