LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ zhptrf()

 subroutine zhptrf ( character UPLO, integer N, complex*16, dimension( * ) AP, integer, dimension( * ) IPIV, integer INFO )

ZHPTRF

Purpose:
``` ZHPTRF computes the factorization of a complex Hermitian packed
matrix A using the Bunch-Kaufman diagonal pivoting method:

A = U*D*U**H  or  A = L*D*L**H

where U (or L) is a product of permutation and unit upper (lower)
triangular matrices, and D is Hermitian and block diagonal with
1-by-1 and 2-by-2 diagonal blocks.```
Parameters
 [in] UPLO ``` UPLO is CHARACTER*1 = 'U': Upper triangle of A is stored; = 'L': Lower triangle of A is stored.``` [in] N ``` N is INTEGER The order of the matrix A. N >= 0.``` [in,out] AP ``` AP is COMPLEX*16 array, dimension (N*(N+1)/2) On entry, the upper or lower triangle of the Hermitian matrix A, packed columnwise in a linear array. The j-th column of A is stored in the array AP as follows: if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. On exit, the block diagonal matrix D and the multipliers used to obtain the factor U or L, stored as a packed triangular matrix overwriting A (see below for further details).``` [out] IPIV ``` IPIV is INTEGER array, dimension (N) Details of the interchanges and the block structure of D. If IPIV(k) > 0, then rows and columns k and IPIV(k) were interchanged and D(k,k) is a 1-by-1 diagonal block. If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k) is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) = IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.``` [out] INFO ``` INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: if INFO = i, D(i,i) is exactly zero. The factorization has been completed, but the block diagonal matrix D is exactly singular, and division by zero will occur if it is used to solve a system of equations.```
Further Details:
```  If UPLO = 'U', then A = U*D*U**H, where
U = P(n)*U(n)* ... *P(k)U(k)* ...,
i.e., U is a product of terms P(k)*U(k), where k decreases from n to
1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
that if the diagonal block D(k) is of order s (s = 1 or 2), then

(   I    v    0   )   k-s
U(k) =  (   0    I    0   )   s
(   0    0    I   )   n-k
k-s   s   n-k

If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
and A(k,k), and v overwrites A(1:k-2,k-1:k).

If UPLO = 'L', then A = L*D*L**H, where
L = P(1)*L(1)* ... *P(k)*L(k)* ...,
i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
that if the diagonal block D(k) is of order s (s = 1 or 2), then

(   I    0     0   )  k-1
L(k) =  (   0    I     0   )  s
(   0    v     I   )  n-k-s+1
k-1   s  n-k-s+1

If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).```
Contributors:
J. Lewis, Boeing Computer Services Company

Definition at line 158 of file zhptrf.f.

159*
160* -- LAPACK computational routine --
161* -- LAPACK is a software package provided by Univ. of Tennessee, --
162* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
163*
164* .. Scalar Arguments ..
165 CHARACTER UPLO
166 INTEGER INFO, N
167* ..
168* .. Array Arguments ..
169 INTEGER IPIV( * )
170 COMPLEX*16 AP( * )
171* ..
172*
173* =====================================================================
174*
175* .. Parameters ..
176 DOUBLE PRECISION ZERO, ONE
177 parameter( zero = 0.0d+0, one = 1.0d+0 )
178 DOUBLE PRECISION EIGHT, SEVTEN
179 parameter( eight = 8.0d+0, sevten = 17.0d+0 )
180* ..
181* .. Local Scalars ..
182 LOGICAL UPPER
183 INTEGER I, IMAX, J, JMAX, K, KC, KK, KNC, KP, KPC,
184 \$ KSTEP, KX, NPP
185 DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D, D11, D22, R1, ROWMAX,
186 \$ TT
187 COMPLEX*16 D12, D21, T, WK, WKM1, WKP1, ZDUM
188* ..
189* .. External Functions ..
190 LOGICAL LSAME
191 INTEGER IZAMAX
192 DOUBLE PRECISION DLAPY2
193 EXTERNAL lsame, izamax, dlapy2
194* ..
195* .. External Subroutines ..
196 EXTERNAL xerbla, zdscal, zhpr, zswap
197* ..
198* .. Intrinsic Functions ..
199 INTRINSIC abs, dble, dcmplx, dconjg, dimag, max, sqrt
200* ..
201* .. Statement Functions ..
202 DOUBLE PRECISION CABS1
203* ..
204* .. Statement Function definitions ..
205 cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
206* ..
207* .. Executable Statements ..
208*
209* Test the input parameters.
210*
211 info = 0
212 upper = lsame( uplo, 'U' )
213 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
214 info = -1
215 ELSE IF( n.LT.0 ) THEN
216 info = -2
217 END IF
218 IF( info.NE.0 ) THEN
219 CALL xerbla( 'ZHPTRF', -info )
220 RETURN
221 END IF
222*
223* Initialize ALPHA for use in choosing pivot block size.
224*
225 alpha = ( one+sqrt( sevten ) ) / eight
226*
227 IF( upper ) THEN
228*
229* Factorize A as U*D*U**H using the upper triangle of A
230*
231* K is the main loop index, decreasing from N to 1 in steps of
232* 1 or 2
233*
234 k = n
235 kc = ( n-1 )*n / 2 + 1
236 10 CONTINUE
237 knc = kc
238*
239* If K < 1, exit from loop
240*
241 IF( k.LT.1 )
242 \$ GO TO 110
243 kstep = 1
244*
245* Determine rows and columns to be interchanged and whether
246* a 1-by-1 or 2-by-2 pivot block will be used
247*
248 absakk = abs( dble( ap( kc+k-1 ) ) )
249*
250* IMAX is the row-index of the largest off-diagonal element in
251* column K, and COLMAX is its absolute value
252*
253 IF( k.GT.1 ) THEN
254 imax = izamax( k-1, ap( kc ), 1 )
255 colmax = cabs1( ap( kc+imax-1 ) )
256 ELSE
257 colmax = zero
258 END IF
259*
260 IF( max( absakk, colmax ).EQ.zero ) THEN
261*
262* Column K is zero: set INFO and continue
263*
264 IF( info.EQ.0 )
265 \$ info = k
266 kp = k
267 ap( kc+k-1 ) = dble( ap( kc+k-1 ) )
268 ELSE
269 IF( absakk.GE.alpha*colmax ) THEN
270*
271* no interchange, use 1-by-1 pivot block
272*
273 kp = k
274 ELSE
275*
276* JMAX is the column-index of the largest off-diagonal
277* element in row IMAX, and ROWMAX is its absolute value
278*
279 rowmax = zero
280 jmax = imax
281 kx = imax*( imax+1 ) / 2 + imax
282 DO 20 j = imax + 1, k
283 IF( cabs1( ap( kx ) ).GT.rowmax ) THEN
284 rowmax = cabs1( ap( kx ) )
285 jmax = j
286 END IF
287 kx = kx + j
288 20 CONTINUE
289 kpc = ( imax-1 )*imax / 2 + 1
290 IF( imax.GT.1 ) THEN
291 jmax = izamax( imax-1, ap( kpc ), 1 )
292 rowmax = max( rowmax, cabs1( ap( kpc+jmax-1 ) ) )
293 END IF
294*
295 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) ) THEN
296*
297* no interchange, use 1-by-1 pivot block
298*
299 kp = k
300 ELSE IF( abs( dble( ap( kpc+imax-1 ) ) ).GE.alpha*
301 \$ rowmax ) THEN
302*
303* interchange rows and columns K and IMAX, use 1-by-1
304* pivot block
305*
306 kp = imax
307 ELSE
308*
309* interchange rows and columns K-1 and IMAX, use 2-by-2
310* pivot block
311*
312 kp = imax
313 kstep = 2
314 END IF
315 END IF
316*
317 kk = k - kstep + 1
318 IF( kstep.EQ.2 )
319 \$ knc = knc - k + 1
320 IF( kp.NE.kk ) THEN
321*
322* Interchange rows and columns KK and KP in the leading
323* submatrix A(1:k,1:k)
324*
325 CALL zswap( kp-1, ap( knc ), 1, ap( kpc ), 1 )
326 kx = kpc + kp - 1
327 DO 30 j = kp + 1, kk - 1
328 kx = kx + j - 1
329 t = dconjg( ap( knc+j-1 ) )
330 ap( knc+j-1 ) = dconjg( ap( kx ) )
331 ap( kx ) = t
332 30 CONTINUE
333 ap( kx+kk-1 ) = dconjg( ap( kx+kk-1 ) )
334 r1 = dble( ap( knc+kk-1 ) )
335 ap( knc+kk-1 ) = dble( ap( kpc+kp-1 ) )
336 ap( kpc+kp-1 ) = r1
337 IF( kstep.EQ.2 ) THEN
338 ap( kc+k-1 ) = dble( ap( kc+k-1 ) )
339 t = ap( kc+k-2 )
340 ap( kc+k-2 ) = ap( kc+kp-1 )
341 ap( kc+kp-1 ) = t
342 END IF
343 ELSE
344 ap( kc+k-1 ) = dble( ap( kc+k-1 ) )
345 IF( kstep.EQ.2 )
346 \$ ap( kc-1 ) = dble( ap( kc-1 ) )
347 END IF
348*
350*
351 IF( kstep.EQ.1 ) THEN
352*
353* 1-by-1 pivot block D(k): column k now holds
354*
355* W(k) = U(k)*D(k)
356*
357* where U(k) is the k-th column of U
358*
359* Perform a rank-1 update of A(1:k-1,1:k-1) as
360*
361* A := A - U(k)*D(k)*U(k)**H = A - W(k)*1/D(k)*W(k)**H
362*
363 r1 = one / dble( ap( kc+k-1 ) )
364 CALL zhpr( uplo, k-1, -r1, ap( kc ), 1, ap )
365*
366* Store U(k) in column k
367*
368 CALL zdscal( k-1, r1, ap( kc ), 1 )
369 ELSE
370*
371* 2-by-2 pivot block D(k): columns k and k-1 now hold
372*
373* ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
374*
375* where U(k) and U(k-1) are the k-th and (k-1)-th columns
376* of U
377*
378* Perform a rank-2 update of A(1:k-2,1:k-2) as
379*
380* A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**H
381* = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )**H
382*
383 IF( k.GT.2 ) THEN
384*
385 d = dlapy2( dble( ap( k-1+( k-1 )*k / 2 ) ),
386 \$ dimag( ap( k-1+( k-1 )*k / 2 ) ) )
387 d22 = dble( ap( k-1+( k-2 )*( k-1 ) / 2 ) ) / d
388 d11 = dble( ap( k+( k-1 )*k / 2 ) ) / d
389 tt = one / ( d11*d22-one )
390 d12 = ap( k-1+( k-1 )*k / 2 ) / d
391 d = tt / d
392*
393 DO 50 j = k - 2, 1, -1
394 wkm1 = d*( d11*ap( j+( k-2 )*( k-1 ) / 2 )-
395 \$ dconjg( d12 )*ap( j+( k-1 )*k / 2 ) )
396 wk = d*( d22*ap( j+( k-1 )*k / 2 )-d12*
397 \$ ap( j+( k-2 )*( k-1 ) / 2 ) )
398 DO 40 i = j, 1, -1
399 ap( i+( j-1 )*j / 2 ) = ap( i+( j-1 )*j / 2 ) -
400 \$ ap( i+( k-1 )*k / 2 )*dconjg( wk ) -
401 \$ ap( i+( k-2 )*( k-1 ) / 2 )*dconjg( wkm1 )
402 40 CONTINUE
403 ap( j+( k-1 )*k / 2 ) = wk
404 ap( j+( k-2 )*( k-1 ) / 2 ) = wkm1
405 ap( j+( j-1 )*j / 2 ) = dcmplx( dble( ap( j+( j-
406 \$ 1 )*j / 2 ) ), 0.0d+0 )
407 50 CONTINUE
408*
409 END IF
410*
411 END IF
412 END IF
413*
414* Store details of the interchanges in IPIV
415*
416 IF( kstep.EQ.1 ) THEN
417 ipiv( k ) = kp
418 ELSE
419 ipiv( k ) = -kp
420 ipiv( k-1 ) = -kp
421 END IF
422*
423* Decrease K and return to the start of the main loop
424*
425 k = k - kstep
426 kc = knc - k
427 GO TO 10
428*
429 ELSE
430*
431* Factorize A as L*D*L**H using the lower triangle of A
432*
433* K is the main loop index, increasing from 1 to N in steps of
434* 1 or 2
435*
436 k = 1
437 kc = 1
438 npp = n*( n+1 ) / 2
439 60 CONTINUE
440 knc = kc
441*
442* If K > N, exit from loop
443*
444 IF( k.GT.n )
445 \$ GO TO 110
446 kstep = 1
447*
448* Determine rows and columns to be interchanged and whether
449* a 1-by-1 or 2-by-2 pivot block will be used
450*
451 absakk = abs( dble( ap( kc ) ) )
452*
453* IMAX is the row-index of the largest off-diagonal element in
454* column K, and COLMAX is its absolute value
455*
456 IF( k.LT.n ) THEN
457 imax = k + izamax( n-k, ap( kc+1 ), 1 )
458 colmax = cabs1( ap( kc+imax-k ) )
459 ELSE
460 colmax = zero
461 END IF
462*
463 IF( max( absakk, colmax ).EQ.zero ) THEN
464*
465* Column K is zero: set INFO and continue
466*
467 IF( info.EQ.0 )
468 \$ info = k
469 kp = k
470 ap( kc ) = dble( ap( kc ) )
471 ELSE
472 IF( absakk.GE.alpha*colmax ) THEN
473*
474* no interchange, use 1-by-1 pivot block
475*
476 kp = k
477 ELSE
478*
479* JMAX is the column-index of the largest off-diagonal
480* element in row IMAX, and ROWMAX is its absolute value
481*
482 rowmax = zero
483 kx = kc + imax - k
484 DO 70 j = k, imax - 1
485 IF( cabs1( ap( kx ) ).GT.rowmax ) THEN
486 rowmax = cabs1( ap( kx ) )
487 jmax = j
488 END IF
489 kx = kx + n - j
490 70 CONTINUE
491 kpc = npp - ( n-imax+1 )*( n-imax+2 ) / 2 + 1
492 IF( imax.LT.n ) THEN
493 jmax = imax + izamax( n-imax, ap( kpc+1 ), 1 )
494 rowmax = max( rowmax, cabs1( ap( kpc+jmax-imax ) ) )
495 END IF
496*
497 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) ) THEN
498*
499* no interchange, use 1-by-1 pivot block
500*
501 kp = k
502 ELSE IF( abs( dble( ap( kpc ) ) ).GE.alpha*rowmax ) THEN
503*
504* interchange rows and columns K and IMAX, use 1-by-1
505* pivot block
506*
507 kp = imax
508 ELSE
509*
510* interchange rows and columns K+1 and IMAX, use 2-by-2
511* pivot block
512*
513 kp = imax
514 kstep = 2
515 END IF
516 END IF
517*
518 kk = k + kstep - 1
519 IF( kstep.EQ.2 )
520 \$ knc = knc + n - k + 1
521 IF( kp.NE.kk ) THEN
522*
523* Interchange rows and columns KK and KP in the trailing
524* submatrix A(k:n,k:n)
525*
526 IF( kp.LT.n )
527 \$ CALL zswap( n-kp, ap( knc+kp-kk+1 ), 1, ap( kpc+1 ),
528 \$ 1 )
529 kx = knc + kp - kk
530 DO 80 j = kk + 1, kp - 1
531 kx = kx + n - j + 1
532 t = dconjg( ap( knc+j-kk ) )
533 ap( knc+j-kk ) = dconjg( ap( kx ) )
534 ap( kx ) = t
535 80 CONTINUE
536 ap( knc+kp-kk ) = dconjg( ap( knc+kp-kk ) )
537 r1 = dble( ap( knc ) )
538 ap( knc ) = dble( ap( kpc ) )
539 ap( kpc ) = r1
540 IF( kstep.EQ.2 ) THEN
541 ap( kc ) = dble( ap( kc ) )
542 t = ap( kc+1 )
543 ap( kc+1 ) = ap( kc+kp-k )
544 ap( kc+kp-k ) = t
545 END IF
546 ELSE
547 ap( kc ) = dble( ap( kc ) )
548 IF( kstep.EQ.2 )
549 \$ ap( knc ) = dble( ap( knc ) )
550 END IF
551*
552* Update the trailing submatrix
553*
554 IF( kstep.EQ.1 ) THEN
555*
556* 1-by-1 pivot block D(k): column k now holds
557*
558* W(k) = L(k)*D(k)
559*
560* where L(k) is the k-th column of L
561*
562 IF( k.LT.n ) THEN
563*
564* Perform a rank-1 update of A(k+1:n,k+1:n) as
565*
566* A := A - L(k)*D(k)*L(k)**H = A - W(k)*(1/D(k))*W(k)**H
567*
568 r1 = one / dble( ap( kc ) )
569 CALL zhpr( uplo, n-k, -r1, ap( kc+1 ), 1,
570 \$ ap( kc+n-k+1 ) )
571*
572* Store L(k) in column K
573*
574 CALL zdscal( n-k, r1, ap( kc+1 ), 1 )
575 END IF
576 ELSE
577*
578* 2-by-2 pivot block D(k): columns K and K+1 now hold
579*
580* ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
581*
582* where L(k) and L(k+1) are the k-th and (k+1)-th columns
583* of L
584*
585 IF( k.LT.n-1 ) THEN
586*
587* Perform a rank-2 update of A(k+2:n,k+2:n) as
588*
589* A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )**H
590* = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )**H
591*
592* where L(k) and L(k+1) are the k-th and (k+1)-th
593* columns of L
594*
595 d = dlapy2( dble( ap( k+1+( k-1 )*( 2*n-k ) / 2 ) ),
596 \$ dimag( ap( k+1+( k-1 )*( 2*n-k ) / 2 ) ) )
597 d11 = dble( ap( k+1+k*( 2*n-k-1 ) / 2 ) ) / d
598 d22 = dble( ap( k+( k-1 )*( 2*n-k ) / 2 ) ) / d
599 tt = one / ( d11*d22-one )
600 d21 = ap( k+1+( k-1 )*( 2*n-k ) / 2 ) / d
601 d = tt / d
602*
603 DO 100 j = k + 2, n
604 wk = d*( d11*ap( j+( k-1 )*( 2*n-k ) / 2 )-d21*
605 \$ ap( j+k*( 2*n-k-1 ) / 2 ) )
606 wkp1 = d*( d22*ap( j+k*( 2*n-k-1 ) / 2 )-
607 \$ dconjg( d21 )*ap( j+( k-1 )*( 2*n-k ) /
608 \$ 2 ) )
609 DO 90 i = j, n
610 ap( i+( j-1 )*( 2*n-j ) / 2 ) = ap( i+( j-1 )*
611 \$ ( 2*n-j ) / 2 ) - ap( i+( k-1 )*( 2*n-k ) /
612 \$ 2 )*dconjg( wk ) - ap( i+k*( 2*n-k-1 ) / 2 )*
613 \$ dconjg( wkp1 )
614 90 CONTINUE
615 ap( j+( k-1 )*( 2*n-k ) / 2 ) = wk
616 ap( j+k*( 2*n-k-1 ) / 2 ) = wkp1
617 ap( j+( j-1 )*( 2*n-j ) / 2 )
618 \$ = dcmplx( dble( ap( j+( j-1 )*( 2*n-j ) / 2 ) ),
619 \$ 0.0d+0 )
620 100 CONTINUE
621 END IF
622 END IF
623 END IF
624*
625* Store details of the interchanges in IPIV
626*
627 IF( kstep.EQ.1 ) THEN
628 ipiv( k ) = kp
629 ELSE
630 ipiv( k ) = -kp
631 ipiv( k+1 ) = -kp
632 END IF
633*
634* Increase K and return to the start of the main loop
635*
636 k = k + kstep
637 kc = knc + n - k + 2
638 GO TO 60
639*
640 END IF
641*
642 110 CONTINUE
643 RETURN
644*
645* End of ZHPTRF
646*
double precision function dlapy2(X, Y)
DLAPY2 returns sqrt(x2+y2).
Definition: dlapy2.f:63
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:71
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:81
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.f:78
subroutine zhpr(UPLO, N, ALPHA, X, INCX, AP)
ZHPR
Definition: zhpr.f:130
Here is the call graph for this function:
Here is the caller graph for this function: