LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine zlahef_rook ( character  UPLO,
integer  N,
integer  NB,
integer  KB,
complex*16, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
complex*16, dimension( ldw, * )  W,
integer  LDW,
integer  INFO 
)

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

Purpose: 
 ZLAHEF_ROOK computes a partial factorization of a complex Hermitian
 matrix A using the bounded Bunch-Kaufman ("rook") diagonal pivoting
 method. The partial factorization has the form:

 A  =  ( I  U12 ) ( A11  0  ) (  I      0     )  if UPLO = 'U', or:
       ( 0  U22 ) (  0   D  ) ( U12**H U22**H )

 A  =  ( L11  0 ) (  D   0  ) ( L11**H L21**H )  if UPLO = 'L'
       ( L21  I ) (  0  A22 ) (  0      I     )

 where the order of D is at most NB. The actual order is returned in
 the argument KB, and is either NB or NB-1, or N if N <= NB.
 Note that U**H denotes the conjugate transpose of U.

 ZLAHEF_ROOK is an auxiliary routine called by ZHETRF_ROOK. It uses
 blocked code (calling Level 3 BLAS) to update the submatrix
 A11 (if UPLO = 'U') or A22 (if UPLO = 'L').
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]NB
          NB is INTEGER
          The maximum number of columns of the matrix A that should be
          factored.  NB should be at least 2 to allow for 2-by-2 pivot
          blocks.
[out]KB
          KB is INTEGER
          The number of columns of A that were actually factored.
          KB is either NB-1 or NB, or N if N <= NB.
[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, A contains details of the partial factorization.
[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':
             Only the last KB elements of IPIV are set.

             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) < 0 and IPIV(k-1) < 0, then rows and
             columns k and -IPIV(k) were interchanged and rows and
             columns k-1 and -IPIV(k-1) were inerchaged,
             D(k-1:k,k-1:k) is a 2-by-2 diagonal block.

          If UPLO = 'L':
             Only the first KB elements of IPIV are set.

             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) < 0 and IPIV(k+1) < 0, then rows and
             columns k and -IPIV(k) were interchanged and rows and
             columns k+1 and -IPIV(k+1) were inerchaged,
             D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
[out]W
          W is COMPLEX*16 array, dimension (LDW,NB)
[in]LDW
          LDW is INTEGER
          The leading dimension of the array W.  LDW >= max(1,N).
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          > 0: if INFO = k, D(k,k) is exactly zero.  The factorization
               has been completed, but the block diagonal matrix D is
               exactly singular.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2013
Contributors: 
  November 2013,  Igor Kozachenko,
                  Computer Science Division,
                  University of California, Berkeley

  September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
                  School of Mathematics,
                  University of Manchester

Definition at line 186 of file zlahef_rook.f.

186 *
187 * -- LAPACK computational routine (version 3.5.0) --
188 * -- LAPACK is a software package provided by Univ. of Tennessee, --
189 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
190 * November 2013
191 *
192 * .. Scalar Arguments ..
193  CHARACTER uplo
194  INTEGER info, kb, lda, ldw, n, nb
195 * ..
196 * .. Array Arguments ..
197  INTEGER ipiv( * )
198  COMPLEX*16 a( lda, * ), w( ldw, * )
199 * ..
200 *
201 * =====================================================================
202 *
203 * .. Parameters ..
204  DOUBLE PRECISION zero, one
205  parameter ( zero = 0.0d+0, one = 1.0d+0 )
206  COMPLEX*16 cone
207  parameter ( cone = ( 1.0d+0, 0.0d+0 ) )
208  DOUBLE PRECISION eight, sevten
209  parameter ( eight = 8.0d+0, sevten = 17.0d+0 )
210 * ..
211 * .. Local Scalars ..
212  LOGICAL done
213  INTEGER imax, itemp, ii, j, jb, jj, jmax, jp1, jp2, k,
214  $ kk, kkw, kp, kstep, kw, p
215  DOUBLE PRECISION absakk, alpha, colmax, dtemp, r1, rowmax, t,
216  $ sfmin
217  COMPLEX*16 d11, d21, d22, z
218 * ..
219 * .. External Functions ..
220  LOGICAL lsame
221  INTEGER izamax
222  DOUBLE PRECISION dlamch
223  EXTERNAL lsame, izamax, dlamch
224 * ..
225 * .. External Subroutines ..
226  EXTERNAL zcopy, zdscal, zgemm, zgemv, zlacgv, zswap
227 * ..
228 * .. Intrinsic Functions ..
229  INTRINSIC abs, dble, dconjg, dimag, max, min, sqrt
230 * ..
231 * .. Statement Functions ..
232  DOUBLE PRECISION cabs1
233 * ..
234 * .. Statement Function definitions ..
235  cabs1( z ) = abs( dble( z ) ) + abs( dimag( z ) )
236 * ..
237 * .. Executable Statements ..
238 *
239  info = 0
240 *
241 * Initialize ALPHA for use in choosing pivot block size.
242 *
243  alpha = ( one+sqrt( sevten ) ) / eight
244 *
245 * Compute machine safe minimum
246 *
247  sfmin = dlamch( 'S' )
248 *
249  IF( lsame( uplo, 'U' ) ) THEN
250 *
251 * Factorize the trailing columns of A using the upper triangle
252 * of A and working backwards, and compute the matrix W = U12*D
253 * for use in updating A11 (note that conjg(W) is actually stored)
254 *
255 * K is the main loop index, decreasing from N in steps of 1 or 2
256 *
257  k = n
258  10 CONTINUE
259 *
260 * KW is the column of W which corresponds to column K of A
261 *
262  kw = nb + k - n
263 *
264 * Exit from loop
265 *
266  IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
267  $ GO TO 30
268 *
269  kstep = 1
270  p = k
271 *
272 * Copy column K of A to column KW of W and update it
273 *
274  IF( k.GT.1 )
275  $ CALL zcopy( k-1, a( 1, k ), 1, w( 1, kw ), 1 )
276  w( k, kw ) = dble( a( k, k ) )
277  IF( k.LT.n ) THEN
278  CALL zgemv( 'No transpose', k, n-k, -cone, a( 1, k+1 ), lda,
279  $ w( k, kw+1 ), ldw, cone, w( 1, kw ), 1 )
280  w( k, kw ) = dble( w( k, kw ) )
281  END IF
282 *
283 * Determine rows and columns to be interchanged and whether
284 * a 1-by-1 or 2-by-2 pivot block will be used
285 *
286  absakk = abs( dble( w( k, kw ) ) )
287 *
288 * IMAX is the row-index of the largest off-diagonal element in
289 * column K, and COLMAX is its absolute value.
290 * Determine both COLMAX and IMAX.
291 *
292  IF( k.GT.1 ) THEN
293  imax = izamax( k-1, w( 1, kw ), 1 )
294  colmax = cabs1( w( imax, kw ) )
295  ELSE
296  colmax = zero
297  END IF
298 *
299  IF( max( absakk, colmax ).EQ.zero ) THEN
300 *
301 * Column K is zero or underflow: set INFO and continue
302 *
303  IF( info.EQ.0 )
304  $ info = k
305  kp = k
306  a( k, k ) = dble( w( k, kw ) )
307  IF( k.GT.1 )
308  $ CALL zcopy( k-1, w( 1, kw ), 1, a( 1, k ), 1 )
309  ELSE
310 *
311 * ============================================================
312 *
313 * BEGIN pivot search
314 *
315 * Case(1)
316 * Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
317 * (used to handle NaN and Inf)
318  IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
319 *
320 * no interchange, use 1-by-1 pivot block
321 *
322  kp = k
323 *
324  ELSE
325 *
326 * Lop until pivot found
327 *
328  done = .false.
329 *
330  12 CONTINUE
331 *
332 * BEGIN pivot search loop body
333 *
334 *
335 * Copy column IMAX to column KW-1 of W and update it
336 *
337  IF( imax.GT.1 )
338  $ CALL zcopy( imax-1, a( 1, imax ), 1, w( 1, kw-1 ),
339  $ 1 )
340  w( imax, kw-1 ) = dble( a( imax, imax ) )
341 *
342  CALL zcopy( k-imax, a( imax, imax+1 ), lda,
343  $ w( imax+1, kw-1 ), 1 )
344  CALL zlacgv( k-imax, w( imax+1, kw-1 ), 1 )
345 *
346  IF( k.LT.n ) THEN
347  CALL zgemv( 'No transpose', k, n-k, -cone,
348  $ a( 1, k+1 ), lda, w( imax, kw+1 ), ldw,
349  $ cone, w( 1, kw-1 ), 1 )
350  w( imax, kw-1 ) = dble( w( imax, kw-1 ) )
351  END IF
352 *
353 * JMAX is the column-index of the largest off-diagonal
354 * element in row IMAX, and ROWMAX is its absolute value.
355 * Determine both ROWMAX and JMAX.
356 *
357  IF( imax.NE.k ) THEN
358  jmax = imax + izamax( k-imax, w( imax+1, kw-1 ),
359  $ 1 )
360  rowmax = cabs1( w( jmax, kw-1 ) )
361  ELSE
362  rowmax = zero
363  END IF
364 *
365  IF( imax.GT.1 ) THEN
366  itemp = izamax( imax-1, w( 1, kw-1 ), 1 )
367  dtemp = cabs1( w( itemp, kw-1 ) )
368  IF( dtemp.GT.rowmax ) THEN
369  rowmax = dtemp
370  jmax = itemp
371  END IF
372  END IF
373 *
374 * Case(2)
375 * Equivalent to testing for
376 * ABS( REAL( W( IMAX,KW-1 ) ) ).GE.ALPHA*ROWMAX
377 * (used to handle NaN and Inf)
378 *
379  IF( .NOT.( abs( dble( w( imax,kw-1 ) ) )
380  $ .LT.alpha*rowmax ) ) THEN
381 *
382 * interchange rows and columns K and IMAX,
383 * use 1-by-1 pivot block
384 *
385  kp = imax
386 *
387 * copy column KW-1 of W to column KW of W
388 *
389  CALL zcopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
390 *
391  done = .true.
392 *
393 * Case(3)
394 * Equivalent to testing for ROWMAX.EQ.COLMAX,
395 * (used to handle NaN and Inf)
396 *
397  ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
398  $ THEN
399 *
400 * interchange rows and columns K-1 and IMAX,
401 * use 2-by-2 pivot block
402 *
403  kp = imax
404  kstep = 2
405  done = .true.
406 *
407 * Case(4)
408  ELSE
409 *
410 * Pivot not found: set params and repeat
411 *
412  p = imax
413  colmax = rowmax
414  imax = jmax
415 *
416 * Copy updated JMAXth (next IMAXth) column to Kth of W
417 *
418  CALL zcopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
419 *
420  END IF
421 *
422 *
423 * END pivot search loop body
424 *
425  IF( .NOT.done ) GOTO 12
426 *
427  END IF
428 *
429 * END pivot search
430 *
431 * ============================================================
432 *
433 * KK is the column of A where pivoting step stopped
434 *
435  kk = k - kstep + 1
436 *
437 * KKW is the column of W which corresponds to column KK of A
438 *
439  kkw = nb + kk - n
440 *
441 * Interchange rows and columns P and K.
442 * Updated column P is already stored in column KW of W.
443 *
444  IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
445 *
446 * Copy non-updated column K to column P of submatrix A
447 * at step K. No need to copy element into columns
448 * K and K-1 of A for 2-by-2 pivot, since these columns
449 * will be later overwritten.
450 *
451  a( p, p ) = dble( a( k, k ) )
452  CALL zcopy( k-1-p, a( p+1, k ), 1, a( p, p+1 ),
453  $ lda )
454  CALL zlacgv( k-1-p, a( p, p+1 ), lda )
455  IF( p.GT.1 )
456  $ CALL zcopy( p-1, a( 1, k ), 1, a( 1, p ), 1 )
457 *
458 * Interchange rows K and P in the last K+1 to N columns of A
459 * (columns K and K-1 of A for 2-by-2 pivot will be
460 * later overwritten). Interchange rows K and P
461 * in last KKW to NB columns of W.
462 *
463  IF( k.LT.n )
464  $ CALL zswap( n-k, a( k, k+1 ), lda, a( p, k+1 ),
465  $ lda )
466  CALL zswap( n-kk+1, w( k, kkw ), ldw, w( p, kkw ),
467  $ ldw )
468  END IF
469 *
470 * Interchange rows and columns KP and KK.
471 * Updated column KP is already stored in column KKW of W.
472 *
473  IF( kp.NE.kk ) THEN
474 *
475 * Copy non-updated column KK to column KP of submatrix A
476 * at step K. No need to copy element into column K
477 * (or K and K-1 for 2-by-2 pivot) of A, since these columns
478 * will be later overwritten.
479 *
480  a( kp, kp ) = dble( a( kk, kk ) )
481  CALL zcopy( kk-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
482  $ lda )
483  CALL zlacgv( kk-1-kp, a( kp, kp+1 ), lda )
484  IF( kp.GT.1 )
485  $ CALL zcopy( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
486 *
487 * Interchange rows KK and KP in last K+1 to N columns of A
488 * (columns K (or K and K-1 for 2-by-2 pivot) of A will be
489 * later overwritten). Interchange rows KK and KP
490 * in last KKW to NB columns of W.
491 *
492  IF( k.LT.n )
493  $ CALL zswap( n-k, a( kk, k+1 ), lda, a( kp, k+1 ),
494  $ lda )
495  CALL zswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
496  $ ldw )
497  END IF
498 *
499  IF( kstep.EQ.1 ) THEN
500 *
501 * 1-by-1 pivot block D(k): column kw of W now holds
502 *
503 * W(kw) = U(k)*D(k),
504 *
505 * where U(k) is the k-th column of U
506 *
507 * (1) Store subdiag. elements of column U(k)
508 * and 1-by-1 block D(k) in column k of A.
509 * (NOTE: Diagonal element U(k,k) is a UNIT element
510 * and not stored)
511 * A(k,k) := D(k,k) = W(k,kw)
512 * A(1:k-1,k) := U(1:k-1,k) = W(1:k-1,kw)/D(k,k)
513 *
514 * (NOTE: No need to use for Hermitian matrix
515 * A( K, K ) = REAL( W( K, K) ) to separately copy diagonal
516 * element D(k,k) from W (potentially saves only one load))
517  CALL zcopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
518  IF( k.GT.1 ) THEN
519 *
520 * (NOTE: No need to check if A(k,k) is NOT ZERO,
521 * since that was ensured earlier in pivot search:
522 * case A(k,k) = 0 falls into 2x2 pivot case(3))
523 *
524 * Handle division by a small number
525 *
526  t = dble( a( k, k ) )
527  IF( abs( t ).GE.sfmin ) THEN
528  r1 = one / t
529  CALL zdscal( k-1, r1, a( 1, k ), 1 )
530  ELSE
531  DO 14 ii = 1, k-1
532  a( ii, k ) = a( ii, k ) / t
533  14 CONTINUE
534  END IF
535 *
536 * (2) Conjugate column W(kw)
537 *
538  CALL zlacgv( k-1, w( 1, kw ), 1 )
539  END IF
540 *
541  ELSE
542 *
543 * 2-by-2 pivot block D(k): columns kw and kw-1 of W now hold
544 *
545 * ( W(kw-1) W(kw) ) = ( U(k-1) U(k) )*D(k)
546 *
547 * where U(k) and U(k-1) are the k-th and (k-1)-th columns
548 * of U
549 *
550 * (1) Store U(1:k-2,k-1) and U(1:k-2,k) and 2-by-2
551 * block D(k-1:k,k-1:k) in columns k-1 and k of A.
552 * (NOTE: 2-by-2 diagonal block U(k-1:k,k-1:k) is a UNIT
553 * block and not stored)
554 * A(k-1:k,k-1:k) := D(k-1:k,k-1:k) = W(k-1:k,kw-1:kw)
555 * A(1:k-2,k-1:k) := U(1:k-2,k:k-1:k) =
556 * = W(1:k-2,kw-1:kw) * ( D(k-1:k,k-1:k)**(-1) )
557 *
558  IF( k.GT.2 ) THEN
559 *
560 * Factor out the columns of the inverse of 2-by-2 pivot
561 * block D, so that each column contains 1, to reduce the
562 * number of FLOPS when we multiply panel
563 * ( W(kw-1) W(kw) ) by this inverse, i.e. by D**(-1).
564 *
565 * D**(-1) = ( d11 cj(d21) )**(-1) =
566 * ( d21 d22 )
567 *
568 * = 1/(d11*d22-|d21|**2) * ( ( d22) (-cj(d21) ) ) =
569 * ( (-d21) ( d11 ) )
570 *
571 * = 1/(|d21|**2) * 1/((d11/cj(d21))*(d22/d21)-1) *
572 *
573 * * ( d21*( d22/d21 ) conj(d21)*( - 1 ) ) =
574 * ( ( -1 ) ( d11/conj(d21) ) )
575 *
576 * = 1/(|d21|**2) * 1/(D22*D11-1) *
577 *
578 * * ( d21*( D11 ) conj(d21)*( -1 ) ) =
579 * ( ( -1 ) ( D22 ) )
580 *
581 * = (1/|d21|**2) * T * ( d21*( D11 ) conj(d21)*( -1 ) ) =
582 * ( ( -1 ) ( D22 ) )
583 *
584 * = ( (T/conj(d21))*( D11 ) (T/d21)*( -1 ) ) =
585 * ( ( -1 ) ( D22 ) )
586 *
587 * Handle division by a small number. (NOTE: order of
588 * operations is important)
589 *
590 * = ( T*(( D11 )/conj(D21)) T*(( -1 )/D21 ) )
591 * ( (( -1 ) ) (( D22 ) ) ),
592 *
593 * where D11 = d22/d21,
594 * D22 = d11/conj(d21),
595 * D21 = d21,
596 * T = 1/(D22*D11-1).
597 *
598 * (NOTE: No need to check for division by ZERO,
599 * since that was ensured earlier in pivot search:
600 * (a) d21 != 0 in 2x2 pivot case(4),
601 * since |d21| should be larger than |d11| and |d22|;
602 * (b) (D22*D11 - 1) != 0, since from (a),
603 * both |D11| < 1, |D22| < 1, hence |D22*D11| << 1.)
604 *
605  d21 = w( k-1, kw )
606  d11 = w( k, kw ) / dconjg( d21 )
607  d22 = w( k-1, kw-1 ) / d21
608  t = one / ( dble( d11*d22 )-one )
609 *
610 * Update elements in columns A(k-1) and A(k) as
611 * dot products of rows of ( W(kw-1) W(kw) ) and columns
612 * of D**(-1)
613 *
614  DO 20 j = 1, k - 2
615  a( j, k-1 ) = t*( ( d11*w( j, kw-1 )-w( j, kw ) ) /
616  $ d21 )
617  a( j, k ) = t*( ( d22*w( j, kw )-w( j, kw-1 ) ) /
618  $ dconjg( d21 ) )
619  20 CONTINUE
620  END IF
621 *
622 * Copy D(k) to A
623 *
624  a( k-1, k-1 ) = w( k-1, kw-1 )
625  a( k-1, k ) = w( k-1, kw )
626  a( k, k ) = w( k, kw )
627 *
628 * (2) Conjugate columns W(kw) and W(kw-1)
629 *
630  CALL zlacgv( k-1, w( 1, kw ), 1 )
631  CALL zlacgv( k-2, w( 1, kw-1 ), 1 )
632 *
633  END IF
634 *
635  END IF
636 *
637 * Store details of the interchanges in IPIV
638 *
639  IF( kstep.EQ.1 ) THEN
640  ipiv( k ) = kp
641  ELSE
642  ipiv( k ) = -p
643  ipiv( k-1 ) = -kp
644  END IF
645 *
646 * Decrease K and return to the start of the main loop
647 *
648  k = k - kstep
649  GO TO 10
650 *
651  30 CONTINUE
652 *
653 * Update the upper triangle of A11 (= A(1:k,1:k)) as
654 *
655 * A11 := A11 - U12*D*U12**H = A11 - U12*W**H
656 *
657 * computing blocks of NB columns at a time (note that conjg(W) is
658 * actually stored)
659 *
660  DO 50 j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
661  jb = min( nb, k-j+1 )
662 *
663 * Update the upper triangle of the diagonal block
664 *
665  DO 40 jj = j, j + jb - 1
666  a( jj, jj ) = dble( a( jj, jj ) )
667  CALL zgemv( 'No transpose', jj-j+1, n-k, -cone,
668  $ a( j, k+1 ), lda, w( jj, kw+1 ), ldw, cone,
669  $ a( j, jj ), 1 )
670  a( jj, jj ) = dble( a( jj, jj ) )
671  40 CONTINUE
672 *
673 * Update the rectangular superdiagonal block
674 *
675  IF( j.GE.2 )
676  $ CALL zgemm( 'No transpose', 'Transpose', j-1, jb, n-k,
677  $ -cone, a( 1, k+1 ), lda, w( j, kw+1 ), ldw,
678  $ cone, a( 1, j ), lda )
679  50 CONTINUE
680 *
681 * Put U12 in standard form by partially undoing the interchanges
682 * in of rows in columns k+1:n looping backwards from k+1 to n
683 *
684  j = k + 1
685  60 CONTINUE
686 *
687 * Undo the interchanges (if any) of rows J and JP2
688 * (or J and JP2, and J+1 and JP1) at each step J
689 *
690  kstep = 1
691  jp1 = 1
692 * (Here, J is a diagonal index)
693  jj = j
694  jp2 = ipiv( j )
695  IF( jp2.LT.0 ) THEN
696  jp2 = -jp2
697 * (Here, J is a diagonal index)
698  j = j + 1
699  jp1 = -ipiv( j )
700  kstep = 2
701  END IF
702 * (NOTE: Here, J is used to determine row length. Length N-J+1
703 * of the rows to swap back doesn't include diagonal element)
704  j = j + 1
705  IF( jp2.NE.jj .AND. j.LE.n )
706  $ CALL zswap( n-j+1, a( jp2, j ), lda, a( jj, j ), lda )
707  jj = jj + 1
708  IF( kstep.EQ.2 .AND. jp1.NE.jj .AND. j.LE.n )
709  $ CALL zswap( n-j+1, a( jp1, j ), lda, a( jj, j ), lda )
710  IF( j.LT.n )
711  $ GO TO 60
712 *
713 * Set KB to the number of columns factorized
714 *
715  kb = n - k
716 *
717  ELSE
718 *
719 * Factorize the leading columns of A using the lower triangle
720 * of A and working forwards, and compute the matrix W = L21*D
721 * for use in updating A22 (note that conjg(W) is actually stored)
722 *
723 * K is the main loop index, increasing from 1 in steps of 1 or 2
724 *
725  k = 1
726  70 CONTINUE
727 *
728 * Exit from loop
729 *
730  IF( ( k.GE.nb .AND. nb.LT.n ) .OR. k.GT.n )
731  $ GO TO 90
732 *
733  kstep = 1
734  p = k
735 *
736 * Copy column K of A to column K of W and update column K of W
737 *
738  w( k, k ) = dble( a( k, k ) )
739  IF( k.LT.n )
740  $ CALL zcopy( n-k, a( k+1, k ), 1, w( k+1, k ), 1 )
741  IF( k.GT.1 ) THEN
742  CALL zgemv( 'No transpose', n-k+1, k-1, -cone, a( k, 1 ),
743  $ lda, w( k, 1 ), ldw, cone, w( k, k ), 1 )
744  w( k, k ) = dble( w( k, k ) )
745  END IF
746 *
747 * Determine rows and columns to be interchanged and whether
748 * a 1-by-1 or 2-by-2 pivot block will be used
749 *
750  absakk = abs( dble( w( k, k ) ) )
751 *
752 * IMAX is the row-index of the largest off-diagonal element in
753 * column K, and COLMAX is its absolute value.
754 * Determine both COLMAX and IMAX.
755 *
756  IF( k.LT.n ) THEN
757  imax = k + izamax( n-k, w( k+1, k ), 1 )
758  colmax = cabs1( w( imax, k ) )
759  ELSE
760  colmax = zero
761  END IF
762 *
763  IF( max( absakk, colmax ).EQ.zero ) THEN
764 *
765 * Column K is zero or underflow: set INFO and continue
766 *
767  IF( info.EQ.0 )
768  $ info = k
769  kp = k
770  a( k, k ) = dble( w( k, k ) )
771  IF( k.LT.n )
772  $ CALL zcopy( n-k, w( k+1, k ), 1, a( k+1, k ), 1 )
773  ELSE
774 *
775 * ============================================================
776 *
777 * BEGIN pivot search
778 *
779 * Case(1)
780 * Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
781 * (used to handle NaN and Inf)
782 *
783  IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
784 *
785 * no interchange, use 1-by-1 pivot block
786 *
787  kp = k
788 *
789  ELSE
790 *
791  done = .false.
792 *
793 * Loop until pivot found
794 *
795  72 CONTINUE
796 *
797 * BEGIN pivot search loop body
798 *
799 *
800 * Copy column IMAX to column k+1 of W and update it
801 *
802  CALL zcopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1)
803  CALL zlacgv( imax-k, w( k, k+1 ), 1 )
804  w( imax, k+1 ) = dble( a( imax, imax ) )
805 *
806  IF( imax.LT.n )
807  $ CALL zcopy( n-imax, a( imax+1, imax ), 1,
808  $ w( imax+1, k+1 ), 1 )
809 *
810  IF( k.GT.1 ) THEN
811  CALL zgemv( 'No transpose', n-k+1, k-1, -cone,
812  $ a( k, 1 ), lda, w( imax, 1 ), ldw,
813  $ cone, w( k, k+1 ), 1 )
814  w( imax, k+1 ) = dble( w( imax, k+1 ) )
815  END IF
816 *
817 * JMAX is the column-index of the largest off-diagonal
818 * element in row IMAX, and ROWMAX is its absolute value.
819 * Determine both ROWMAX and JMAX.
820 *
821  IF( imax.NE.k ) THEN
822  jmax = k - 1 + izamax( imax-k, w( k, k+1 ), 1 )
823  rowmax = cabs1( w( jmax, k+1 ) )
824  ELSE
825  rowmax = zero
826  END IF
827 *
828  IF( imax.LT.n ) THEN
829  itemp = imax + izamax( n-imax, w( imax+1, k+1 ), 1)
830  dtemp = cabs1( w( itemp, k+1 ) )
831  IF( dtemp.GT.rowmax ) THEN
832  rowmax = dtemp
833  jmax = itemp
834  END IF
835  END IF
836 *
837 * Case(2)
838 * Equivalent to testing for
839 * ABS( REAL( W( IMAX,K+1 ) ) ).GE.ALPHA*ROWMAX
840 * (used to handle NaN and Inf)
841 *
842  IF( .NOT.( abs( dble( w( imax,k+1 ) ) )
843  $ .LT.alpha*rowmax ) ) THEN
844 *
845 * interchange rows and columns K and IMAX,
846 * use 1-by-1 pivot block
847 *
848  kp = imax
849 *
850 * copy column K+1 of W to column K of W
851 *
852  CALL zcopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
853 *
854  done = .true.
855 *
856 * Case(3)
857 * Equivalent to testing for ROWMAX.EQ.COLMAX,
858 * (used to handle NaN and Inf)
859 *
860  ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
861  $ THEN
862 *
863 * interchange rows and columns K+1 and IMAX,
864 * use 2-by-2 pivot block
865 *
866  kp = imax
867  kstep = 2
868  done = .true.
869 *
870 * Case(4)
871  ELSE
872 *
873 * Pivot not found: set params and repeat
874 *
875  p = imax
876  colmax = rowmax
877  imax = jmax
878 *
879 * Copy updated JMAXth (next IMAXth) column to Kth of W
880 *
881  CALL zcopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
882 *
883  END IF
884 *
885 *
886 * End pivot search loop body
887 *
888  IF( .NOT.done ) GOTO 72
889 *
890  END IF
891 *
892 * END pivot search
893 *
894 * ============================================================
895 *
896 * KK is the column of A where pivoting step stopped
897 *
898  kk = k + kstep - 1
899 *
900 * Interchange rows and columns P and K (only for 2-by-2 pivot).
901 * Updated column P is already stored in column K of W.
902 *
903  IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
904 *
905 * Copy non-updated column KK-1 to column P of submatrix A
906 * at step K. No need to copy element into columns
907 * K and K+1 of A for 2-by-2 pivot, since these columns
908 * will be later overwritten.
909 *
910  a( p, p ) = dble( a( k, k ) )
911  CALL zcopy( p-k-1, a( k+1, k ), 1, a( p, k+1 ), lda )
912  CALL zlacgv( p-k-1, a( p, k+1 ), lda )
913  IF( p.LT.n )
914  $ CALL zcopy( n-p, a( p+1, k ), 1, a( p+1, p ), 1 )
915 *
916 * Interchange rows K and P in first K-1 columns of A
917 * (columns K and K+1 of A for 2-by-2 pivot will be
918 * later overwritten). Interchange rows K and P
919 * in first KK columns of W.
920 *
921  IF( k.GT.1 )
922  $ CALL zswap( k-1, a( k, 1 ), lda, a( p, 1 ), lda )
923  CALL zswap( kk, w( k, 1 ), ldw, w( p, 1 ), ldw )
924  END IF
925 *
926 * Interchange rows and columns KP and KK.
927 * Updated column KP is already stored in column KK of W.
928 *
929  IF( kp.NE.kk ) THEN
930 *
931 * Copy non-updated column KK to column KP of submatrix A
932 * at step K. No need to copy element into column K
933 * (or K and K+1 for 2-by-2 pivot) of A, since these columns
934 * will be later overwritten.
935 *
936  a( kp, kp ) = dble( a( kk, kk ) )
937  CALL zcopy( kp-kk-1, a( kk+1, kk ), 1, a( kp, kk+1 ),
938  $ lda )
939  CALL zlacgv( kp-kk-1, a( kp, kk+1 ), lda )
940  IF( kp.LT.n )
941  $ CALL zcopy( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
942 *
943 * Interchange rows KK and KP in first K-1 columns of A
944 * (column K (or K and K+1 for 2-by-2 pivot) of A will be
945 * later overwritten). Interchange rows KK and KP
946 * in first KK columns of W.
947 *
948  IF( k.GT.1 )
949  $ CALL zswap( k-1, a( kk, 1 ), lda, a( kp, 1 ), lda )
950  CALL zswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
951  END IF
952 *
953  IF( kstep.EQ.1 ) THEN
954 *
955 * 1-by-1 pivot block D(k): column k of W now holds
956 *
957 * W(k) = L(k)*D(k),
958 *
959 * where L(k) is the k-th column of L
960 *
961 * (1) Store subdiag. elements of column L(k)
962 * and 1-by-1 block D(k) in column k of A.
963 * (NOTE: Diagonal element L(k,k) is a UNIT element
964 * and not stored)
965 * A(k,k) := D(k,k) = W(k,k)
966 * A(k+1:N,k) := L(k+1:N,k) = W(k+1:N,k)/D(k,k)
967 *
968 * (NOTE: No need to use for Hermitian matrix
969 * A( K, K ) = REAL( W( K, K) ) to separately copy diagonal
970 * element D(k,k) from W (potentially saves only one load))
971  CALL zcopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
972  IF( k.LT.n ) THEN
973 *
974 * (NOTE: No need to check if A(k,k) is NOT ZERO,
975 * since that was ensured earlier in pivot search:
976 * case A(k,k) = 0 falls into 2x2 pivot case(3))
977 *
978 * Handle division by a small number
979 *
980  t = dble( a( k, k ) )
981  IF( abs( t ).GE.sfmin ) THEN
982  r1 = one / t
983  CALL zdscal( n-k, r1, a( k+1, k ), 1 )
984  ELSE
985  DO 74 ii = k + 1, n
986  a( ii, k ) = a( ii, k ) / t
987  74 CONTINUE
988  END IF
989 *
990 * (2) Conjugate column W(k)
991 *
992  CALL zlacgv( n-k, w( k+1, k ), 1 )
993  END IF
994 *
995  ELSE
996 *
997 * 2-by-2 pivot block D(k): columns k and k+1 of W now hold
998 *
999 * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
1000 *
1001 * where L(k) and L(k+1) are the k-th and (k+1)-th columns
1002 * of L
1003 *
1004 * (1) Store L(k+2:N,k) and L(k+2:N,k+1) and 2-by-2
1005 * block D(k:k+1,k:k+1) in columns k and k+1 of A.
1006 * NOTE: 2-by-2 diagonal block L(k:k+1,k:k+1) is a UNIT
1007 * block and not stored.
1008 * A(k:k+1,k:k+1) := D(k:k+1,k:k+1) = W(k:k+1,k:k+1)
1009 * A(k+2:N,k:k+1) := L(k+2:N,k:k+1) =
1010 * = W(k+2:N,k:k+1) * ( D(k:k+1,k:k+1)**(-1) )
1011 *
1012  IF( k.LT.n-1 ) THEN
1013 *
1014 * Factor out the columns of the inverse of 2-by-2 pivot
1015 * block D, so that each column contains 1, to reduce the
1016 * number of FLOPS when we multiply panel
1017 * ( W(kw-1) W(kw) ) by this inverse, i.e. by D**(-1).
1018 *
1019 * D**(-1) = ( d11 cj(d21) )**(-1) =
1020 * ( d21 d22 )
1021 *
1022 * = 1/(d11*d22-|d21|**2) * ( ( d22) (-cj(d21) ) ) =
1023 * ( (-d21) ( d11 ) )
1024 *
1025 * = 1/(|d21|**2) * 1/((d11/cj(d21))*(d22/d21)-1) *
1026 *
1027 * * ( d21*( d22/d21 ) conj(d21)*( - 1 ) ) =
1028 * ( ( -1 ) ( d11/conj(d21) ) )
1029 *
1030 * = 1/(|d21|**2) * 1/(D22*D11-1) *
1031 *
1032 * * ( d21*( D11 ) conj(d21)*( -1 ) ) =
1033 * ( ( -1 ) ( D22 ) )
1034 *
1035 * = (1/|d21|**2) * T * ( d21*( D11 ) conj(d21)*( -1 ) ) =
1036 * ( ( -1 ) ( D22 ) )
1037 *
1038 * = ( (T/conj(d21))*( D11 ) (T/d21)*( -1 ) ) =
1039 * ( ( -1 ) ( D22 ) )
1040 *
1041 * Handle division by a small number. (NOTE: order of
1042 * operations is important)
1043 *
1044 * = ( T*(( D11 )/conj(D21)) T*(( -1 )/D21 ) )
1045 * ( (( -1 ) ) (( D22 ) ) ),
1046 *
1047 * where D11 = d22/d21,
1048 * D22 = d11/conj(d21),
1049 * D21 = d21,
1050 * T = 1/(D22*D11-1).
1051 *
1052 * (NOTE: No need to check for division by ZERO,
1053 * since that was ensured earlier in pivot search:
1054 * (a) d21 != 0 in 2x2 pivot case(4),
1055 * since |d21| should be larger than |d11| and |d22|;
1056 * (b) (D22*D11 - 1) != 0, since from (a),
1057 * both |D11| < 1, |D22| < 1, hence |D22*D11| << 1.)
1058 *
1059  d21 = w( k+1, k )
1060  d11 = w( k+1, k+1 ) / d21
1061  d22 = w( k, k ) / dconjg( d21 )
1062  t = one / ( dble( d11*d22 )-one )
1063 *
1064 * Update elements in columns A(k) and A(k+1) as
1065 * dot products of rows of ( W(k) W(k+1) ) and columns
1066 * of D**(-1)
1067 *
1068  DO 80 j = k + 2, n
1069  a( j, k ) = t*( ( d11*w( j, k )-w( j, k+1 ) ) /
1070  $ dconjg( d21 ) )
1071  a( j, k+1 ) = t*( ( d22*w( j, k+1 )-w( j, k ) ) /
1072  $ d21 )
1073  80 CONTINUE
1074  END IF
1075 *
1076 * Copy D(k) to A
1077 *
1078  a( k, k ) = w( k, k )
1079  a( k+1, k ) = w( k+1, k )
1080  a( k+1, k+1 ) = w( k+1, k+1 )
1081 *
1082 * (2) Conjugate columns W(k) and W(k+1)
1083 *
1084  CALL zlacgv( n-k, w( k+1, k ), 1 )
1085  CALL zlacgv( n-k-1, w( k+2, k+1 ), 1 )
1086 *
1087  END IF
1088 *
1089  END IF
1090 *
1091 * Store details of the interchanges in IPIV
1092 *
1093  IF( kstep.EQ.1 ) THEN
1094  ipiv( k ) = kp
1095  ELSE
1096  ipiv( k ) = -p
1097  ipiv( k+1 ) = -kp
1098  END IF
1099 *
1100 * Increase K and return to the start of the main loop
1101 *
1102  k = k + kstep
1103  GO TO 70
1104 *
1105  90 CONTINUE
1106 *
1107 * Update the lower triangle of A22 (= A(k:n,k:n)) as
1108 *
1109 * A22 := A22 - L21*D*L21**H = A22 - L21*W**H
1110 *
1111 * computing blocks of NB columns at a time (note that conjg(W) is
1112 * actually stored)
1113 *
1114  DO 110 j = k, n, nb
1115  jb = min( nb, n-j+1 )
1116 *
1117 * Update the lower triangle of the diagonal block
1118 *
1119  DO 100 jj = j, j + jb - 1
1120  a( jj, jj ) = dble( a( jj, jj ) )
1121  CALL zgemv( 'No transpose', j+jb-jj, k-1, -cone,
1122  $ a( jj, 1 ), lda, w( jj, 1 ), ldw, cone,
1123  $ a( jj, jj ), 1 )
1124  a( jj, jj ) = dble( a( jj, jj ) )
1125  100 CONTINUE
1126 *
1127 * Update the rectangular subdiagonal block
1128 *
1129  IF( j+jb.LE.n )
1130  $ CALL zgemm( 'No transpose', 'Transpose', n-j-jb+1, jb,
1131  $ k-1, -cone, a( j+jb, 1 ), lda, w( j, 1 ),
1132  $ ldw, cone, a( j+jb, j ), lda )
1133  110 CONTINUE
1134 *
1135 * Put L21 in standard form by partially undoing the interchanges
1136 * of rows in columns 1:k-1 looping backwards from k-1 to 1
1137 *
1138  j = k - 1
1139  120 CONTINUE
1140 *
1141 * Undo the interchanges (if any) of rows J and JP2
1142 * (or J and JP2, and J-1 and JP1) at each step J
1143 *
1144  kstep = 1
1145  jp1 = 1
1146 * (Here, J is a diagonal index)
1147  jj = j
1148  jp2 = ipiv( j )
1149  IF( jp2.LT.0 ) THEN
1150  jp2 = -jp2
1151 * (Here, J is a diagonal index)
1152  j = j - 1
1153  jp1 = -ipiv( j )
1154  kstep = 2
1155  END IF
1156 * (NOTE: Here, J is used to determine row length. Length J
1157 * of the rows to swap back doesn't include diagonal element)
1158  j = j - 1
1159  IF( jp2.NE.jj .AND. j.GE.1 )
1160  $ CALL zswap( j, a( jp2, 1 ), lda, a( jj, 1 ), lda )
1161  jj = jj -1
1162  IF( kstep.EQ.2 .AND. jp1.NE.jj .AND. j.GE.1 )
1163  $ CALL zswap( j, a( jp1, 1 ), lda, a( jj, 1 ), lda )
1164  IF( j.GT.1 )
1165  $ GO TO 120
1166 *
1167 * Set KB to the number of columns factorized
1168 *
1169  kb = k - 1
1170 *
1171  END IF
1172  RETURN
1173 *
1174 * End of ZLAHEF_ROOK
1175 *
zcopy
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:52
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
zgemv
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:160
zswap
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:52
zgemm
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
izamax
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:53
zdscal
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.f:54
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
zlacgv
subroutine zlacgv(N, X, INCX)
ZLACGV conjugates a complex vector.
Definition: zlacgv.f:76

Here is the call graph for this function:

Here is the caller graph for this function: