LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine zhetf2 ( character  UPLO,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
integer  INFO 
)

ZHETF2 computes the factorization of a complex Hermitian matrix, using the diagonal pivoting method (unblocked algorithm, calling Level 2 BLAS).

Download ZHETF2 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZHETF2 computes the factorization of a complex Hermitian 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, U**H is the conjugate transpose of U, and D is
 Hermitian and block diagonal with 1-by-1 and 2-by-2 diagonal blocks.

 This is the unblocked version of the algorithm, calling Level 2 BLAS.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          Specifies whether the upper or lower triangular part of the
          Hermitian matrix A is stored:
          = 'U':  Upper triangular
          = 'L':  Lower triangular
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading
          n-by-n upper triangular part of A contains the upper
          triangular part of the matrix A, and the strictly lower
          triangular part of A is not referenced.  If UPLO = 'L', the
          leading n-by-n lower triangular part of A contains the lower
          triangular part of the matrix A, and the strictly upper
          triangular part of A is not referenced.

          On exit, the block diagonal matrix D and the multipliers used
          to obtain the factor U or L (see below for further details).
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[out]IPIV
          IPIV is INTEGER array, dimension (N)
          Details of the interchanges and the block structure of D.

          If UPLO = 'U':
             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 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':
             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 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 = -k, the k-th argument had an illegal value
          > 0: if INFO = k, D(k,k) 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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2013
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:
  09-29-06 - patch from
    Bobby Cheng, MathWorks

    Replace l.210 and l.393
         IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
    by
         IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN

  01-01-96 - Based on modifications by
    J. Lewis, Boeing Computer Services Company
    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA

Definition at line 193 of file zhetf2.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: