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

CLAHEF computes a partial factorization of a complex Hermitian indefinite matrix using the Bunch-Kaufman diagonal pivoting method (blocked algorithm, calling Level 3 BLAS).

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

Purpose: 
 CLAHEF computes a partial factorization of a complex Hermitian
 matrix A using the Bunch-Kaufman 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.

 CLAHEF is an auxiliary routine called by CHETRF. 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 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) = 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':
             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) = 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]W
          W is COMPLEX 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

Definition at line 179 of file clahef.f.

179 *
180 * -- LAPACK computational routine (version 3.5.0) --
181 * -- LAPACK is a software package provided by Univ. of Tennessee, --
182 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
183 * November 2013
184 *
185 * .. Scalar Arguments ..
186  CHARACTER uplo
187  INTEGER info, kb, lda, ldw, n, nb
188 * ..
189 * .. Array Arguments ..
190  INTEGER ipiv( * )
191  COMPLEX a( lda, * ), w( ldw, * )
192 * ..
193 *
194 * =====================================================================
195 *
196 * .. Parameters ..
197  REAL zero, one
198  parameter ( zero = 0.0e+0, one = 1.0e+0 )
199  COMPLEX cone
200  parameter ( cone = ( 1.0e+0, 0.0e+0 ) )
201  REAL eight, sevten
202  parameter ( eight = 8.0e+0, sevten = 17.0e+0 )
203 * ..
204 * .. Local Scalars ..
205  INTEGER imax, j, jb, jj, jmax, jp, k, kk, kkw, kp,
206  $ kstep, kw
207  REAL absakk, alpha, colmax, r1, rowmax, t
208  COMPLEX d11, d21, d22, z
209 * ..
210 * .. External Functions ..
211  LOGICAL lsame
212  INTEGER icamax
213  EXTERNAL lsame, icamax
214 * ..
215 * .. External Subroutines ..
216  EXTERNAL ccopy, cgemm, cgemv, clacgv, csscal, cswap
217 * ..
218 * .. Intrinsic Functions ..
219  INTRINSIC abs, aimag, conjg, max, min, REAL, sqrt
220 * ..
221 * .. Statement Functions ..
222  REAL cabs1
223 * ..
224 * .. Statement Function definitions ..
225  cabs1( z ) = abs( REAL( Z ) ) + abs( aimag( z ) )
226 * ..
227 * .. Executable Statements ..
228 *
229  info = 0
230 *
231 * Initialize ALPHA for use in choosing pivot block size.
232 *
233  alpha = ( one+sqrt( sevten ) ) / eight
234 *
235  IF( lsame( uplo, 'U' ) ) THEN
236 *
237 * Factorize the trailing columns of A using the upper triangle
238 * of A and working backwards, and compute the matrix W = U12*D
239 * for use in updating A11 (note that conjg(W) is actually stored)
240 *
241 * K is the main loop index, decreasing from N in steps of 1 or 2
242 *
243  k = n
244  10 CONTINUE
245 *
246 * KW is the column of W which corresponds to column K of A
247 *
248  kw = nb + k - n
249 *
250 * Exit from loop
251 *
252  IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
253  $ GO TO 30
254 *
255  kstep = 1
256 *
257 * Copy column K of A to column KW of W and update it
258 *
259  CALL ccopy( k-1, a( 1, k ), 1, w( 1, kw ), 1 )
260  w( k, kw ) = REAL( A( K, K ) )
261  IF( k.LT.n ) THEN
262  CALL cgemv( 'No transpose', k, n-k, -cone, a( 1, k+1 ), lda,
263  $ w( k, kw+1 ), ldw, cone, w( 1, kw ), 1 )
264  w( k, kw ) = REAL( W( K, KW ) )
265  END IF
266 *
267 * Determine rows and columns to be interchanged and whether
268 * a 1-by-1 or 2-by-2 pivot block will be used
269 *
270  absakk = abs( REAL( W( K, KW ) ) )
271 *
272 * IMAX is the row-index of the largest off-diagonal element in
273 * column K, and COLMAX is its absolute value.
274 * Determine both COLMAX and IMAX.
275 *
276  IF( k.GT.1 ) THEN
277  imax = icamax( k-1, w( 1, kw ), 1 )
278  colmax = cabs1( w( imax, kw ) )
279  ELSE
280  colmax = zero
281  END IF
282 *
283  IF( max( absakk, colmax ).EQ.zero ) THEN
284 *
285 * Column K is zero or underflow: set INFO and continue
286 *
287  IF( info.EQ.0 )
288  $ info = k
289  kp = k
290  a( k, k ) = REAL( A( K, K ) )
291  ELSE
292 *
293 * ============================================================
294 *
295 * BEGIN pivot search
296 *
297 * Case(1)
298  IF( absakk.GE.alpha*colmax ) THEN
299 *
300 * no interchange, use 1-by-1 pivot block
301 *
302  kp = k
303  ELSE
304 *
305 * BEGIN pivot search along IMAX row
306 *
307 *
308 * Copy column IMAX to column KW-1 of W and update it
309 *
310  CALL ccopy( imax-1, a( 1, imax ), 1, w( 1, kw-1 ), 1 )
311  w( imax, kw-1 ) = REAL( A( IMAX, IMAX ) )
312  CALL ccopy( k-imax, a( imax, imax+1 ), lda,
313  $ w( imax+1, kw-1 ), 1 )
314  CALL clacgv( k-imax, w( imax+1, kw-1 ), 1 )
315  IF( k.LT.n ) THEN
316  CALL cgemv( 'No transpose', k, n-k, -cone,
317  $ a( 1, k+1 ), lda, w( imax, kw+1 ), ldw,
318  $ cone, w( 1, kw-1 ), 1 )
319  w( imax, kw-1 ) = REAL( W( IMAX, KW-1 ) )
320  END IF
321 *
322 * JMAX is the column-index of the largest off-diagonal
323 * element in row IMAX, and ROWMAX is its absolute value.
324 * Determine only ROWMAX.
325 *
326  jmax = imax + icamax( k-imax, w( imax+1, kw-1 ), 1 )
327  rowmax = cabs1( w( jmax, kw-1 ) )
328  IF( imax.GT.1 ) THEN
329  jmax = icamax( imax-1, w( 1, kw-1 ), 1 )
330  rowmax = max( rowmax, cabs1( w( jmax, kw-1 ) ) )
331  END IF
332 *
333 * Case(2)
334  IF( absakk.GE.alpha*colmax*( colmax / rowmax ) ) THEN
335 *
336 * no interchange, use 1-by-1 pivot block
337 *
338  kp = k
339 *
340 * Case(3)
341  ELSE IF( abs( REAL( W( IMAX, KW-1 ) ) ).GE.alpha*rowmax )
342  $ THEN
343 *
344 * interchange rows and columns K and IMAX, use 1-by-1
345 * pivot block
346 *
347  kp = imax
348 *
349 * copy column KW-1 of W to column KW of W
350 *
351  CALL ccopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
352 *
353 * Case(4)
354  ELSE
355 *
356 * interchange rows and columns K-1 and IMAX, use 2-by-2
357 * pivot block
358 *
359  kp = imax
360  kstep = 2
361  END IF
362 *
363 *
364 * END pivot search along IMAX row
365 *
366  END IF
367 *
368 * END pivot search
369 *
370 * ============================================================
371 *
372 * KK is the column of A where pivoting step stopped
373 *
374  kk = k - kstep + 1
375 *
376 * KKW is the column of W which corresponds to column KK of A
377 *
378  kkw = nb + kk - n
379 *
380 * Interchange rows and columns KP and KK.
381 * Updated column KP is already stored in column KKW of W.
382 *
383  IF( kp.NE.kk ) THEN
384 *
385 * Copy non-updated column KK to column KP of submatrix A
386 * at step K. No need to copy element into column K
387 * (or K and K-1 for 2-by-2 pivot) of A, since these columns
388 * will be later overwritten.
389 *
390  a( kp, kp ) = REAL( A( KK, KK ) )
391  CALL ccopy( kk-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
392  $ lda )
393  CALL clacgv( kk-1-kp, a( kp, kp+1 ), lda )
394  IF( kp.GT.1 )
395  $ CALL ccopy( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
396 *
397 * Interchange rows KK and KP in last K+1 to N columns of A
398 * (columns K (or K and K-1 for 2-by-2 pivot) of A will be
399 * later overwritten). Interchange rows KK and KP
400 * in last KKW to NB columns of W.
401 *
402  IF( k.LT.n )
403  $ CALL cswap( n-k, a( kk, k+1 ), lda, a( kp, k+1 ),
404  $ lda )
405  CALL cswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
406  $ ldw )
407  END IF
408 *
409  IF( kstep.EQ.1 ) THEN
410 *
411 * 1-by-1 pivot block D(k): column kw of W now holds
412 *
413 * W(kw) = U(k)*D(k),
414 *
415 * where U(k) is the k-th column of U
416 *
417 * (1) Store subdiag. elements of column U(k)
418 * and 1-by-1 block D(k) in column k of A.
419 * (NOTE: Diagonal element U(k,k) is a UNIT element
420 * and not stored)
421 * A(k,k) := D(k,k) = W(k,kw)
422 * A(1:k-1,k) := U(1:k-1,k) = W(1:k-1,kw)/D(k,k)
423 *
424 * (NOTE: No need to use for Hermitian matrix
425 * A( K, K ) = DBLE( W( K, K) ) to separately copy diagonal
426 * element D(k,k) from W (potentially saves only one load))
427  CALL ccopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
428  IF( k.GT.1 ) THEN
429 *
430 * (NOTE: No need to check if A(k,k) is NOT ZERO,
431 * since that was ensured earlier in pivot search:
432 * case A(k,k) = 0 falls into 2x2 pivot case(4))
433 *
434  r1 = one / REAL( A( K, K ) )
435  CALL csscal( k-1, r1, a( 1, k ), 1 )
436 *
437 * (2) Conjugate column W(kw)
438 *
439  CALL clacgv( k-1, w( 1, kw ), 1 )
440  END IF
441 *
442  ELSE
443 *
444 * 2-by-2 pivot block D(k): columns kw and kw-1 of W now hold
445 *
446 * ( W(kw-1) W(kw) ) = ( U(k-1) U(k) )*D(k)
447 *
448 * where U(k) and U(k-1) are the k-th and (k-1)-th columns
449 * of U
450 *
451 * (1) Store U(1:k-2,k-1) and U(1:k-2,k) and 2-by-2
452 * block D(k-1:k,k-1:k) in columns k-1 and k of A.
453 * (NOTE: 2-by-2 diagonal block U(k-1:k,k-1:k) is a UNIT
454 * block and not stored)
455 * A(k-1:k,k-1:k) := D(k-1:k,k-1:k) = W(k-1:k,kw-1:kw)
456 * A(1:k-2,k-1:k) := U(1:k-2,k:k-1:k) =
457 * = W(1:k-2,kw-1:kw) * ( D(k-1:k,k-1:k)**(-1) )
458 *
459  IF( k.GT.2 ) THEN
460 *
461 * Factor out the columns of the inverse of 2-by-2 pivot
462 * block D, so that each column contains 1, to reduce the
463 * number of FLOPS when we multiply panel
464 * ( W(kw-1) W(kw) ) by this inverse, i.e. by D**(-1).
465 *
466 * D**(-1) = ( d11 cj(d21) )**(-1) =
467 * ( d21 d22 )
468 *
469 * = 1/(d11*d22-|d21|**2) * ( ( d22) (-cj(d21) ) ) =
470 * ( (-d21) ( d11 ) )
471 *
472 * = 1/(|d21|**2) * 1/((d11/cj(d21))*(d22/d21)-1) *
473 *
474 * * ( d21*( d22/d21 ) conj(d21)*( - 1 ) ) =
475 * ( ( -1 ) ( d11/conj(d21) ) )
476 *
477 * = 1/(|d21|**2) * 1/(D22*D11-1) *
478 *
479 * * ( d21*( D11 ) conj(d21)*( -1 ) ) =
480 * ( ( -1 ) ( D22 ) )
481 *
482 * = (1/|d21|**2) * T * ( d21*( D11 ) conj(d21)*( -1 ) ) =
483 * ( ( -1 ) ( D22 ) )
484 *
485 * = ( (T/conj(d21))*( D11 ) (T/d21)*( -1 ) ) =
486 * ( ( -1 ) ( D22 ) )
487 *
488 * = ( conj(D21)*( D11 ) D21*( -1 ) )
489 * ( ( -1 ) ( D22 ) ),
490 *
491 * where D11 = d22/d21,
492 * D22 = d11/conj(d21),
493 * D21 = T/d21,
494 * T = 1/(D22*D11-1).
495 *
496 * (NOTE: No need to check for division by ZERO,
497 * since that was ensured earlier in pivot search:
498 * (a) d21 != 0, since in 2x2 pivot case(4)
499 * |d21| should be larger than |d11| and |d22|;
500 * (b) (D22*D11 - 1) != 0, since from (a),
501 * both |D11| < 1, |D22| < 1, hence |D22*D11| << 1.)
502 *
503  d21 = w( k-1, kw )
504  d11 = w( k, kw ) / conjg( d21 )
505  d22 = w( k-1, kw-1 ) / d21
506  t = one / ( REAL( d11*d22 )-one )
507  d21 = t / d21
508 *
509 * Update elements in columns A(k-1) and A(k) as
510 * dot products of rows of ( W(kw-1) W(kw) ) and columns
511 * of D**(-1)
512 *
513  DO 20 j = 1, k - 2
514  a( j, k-1 ) = d21*( d11*w( j, kw-1 )-w( j, kw ) )
515  a( j, k ) = conjg( d21 )*
516  $ ( d22*w( j, kw )-w( j, kw-1 ) )
517  20 CONTINUE
518  END IF
519 *
520 * Copy D(k) to A
521 *
522  a( k-1, k-1 ) = w( k-1, kw-1 )
523  a( k-1, k ) = w( k-1, kw )
524  a( k, k ) = w( k, kw )
525 *
526 * (2) Conjugate columns W(kw) and W(kw-1)
527 *
528  CALL clacgv( k-1, w( 1, kw ), 1 )
529  CALL clacgv( k-2, w( 1, kw-1 ), 1 )
530 *
531  END IF
532 *
533  END IF
534 *
535 * Store details of the interchanges in IPIV
536 *
537  IF( kstep.EQ.1 ) THEN
538  ipiv( k ) = kp
539  ELSE
540  ipiv( k ) = -kp
541  ipiv( k-1 ) = -kp
542  END IF
543 *
544 * Decrease K and return to the start of the main loop
545 *
546  k = k - kstep
547  GO TO 10
548 *
549  30 CONTINUE
550 *
551 * Update the upper triangle of A11 (= A(1:k,1:k)) as
552 *
553 * A11 := A11 - U12*D*U12**H = A11 - U12*W**H
554 *
555 * computing blocks of NB columns at a time (note that conjg(W) is
556 * actually stored)
557 *
558  DO 50 j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
559  jb = min( nb, k-j+1 )
560 *
561 * Update the upper triangle of the diagonal block
562 *
563  DO 40 jj = j, j + jb - 1
564  a( jj, jj ) = REAL( A( JJ, JJ ) )
565  CALL cgemv( 'No transpose', jj-j+1, n-k, -cone,
566  $ a( j, k+1 ), lda, w( jj, kw+1 ), ldw, cone,
567  $ a( j, jj ), 1 )
568  a( jj, jj ) = REAL( A( JJ, JJ ) )
569  40 CONTINUE
570 *
571 * Update the rectangular superdiagonal block
572 *
573  CALL cgemm( 'No transpose', 'Transpose', j-1, jb, n-k,
574  $ -cone, a( 1, k+1 ), lda, w( j, kw+1 ), ldw,
575  $ cone, a( 1, j ), lda )
576  50 CONTINUE
577 *
578 * Put U12 in standard form by partially undoing the interchanges
579 * in of rows in columns k+1:n looping backwards from k+1 to n
580 *
581  j = k + 1
582  60 CONTINUE
583 *
584 * Undo the interchanges (if any) of rows J and JP
585 * at each step J
586 *
587 * (Here, J is a diagonal index)
588  jj = j
589  jp = ipiv( j )
590  IF( jp.LT.0 ) THEN
591  jp = -jp
592 * (Here, J is a diagonal index)
593  j = j + 1
594  END IF
595 * (NOTE: Here, J is used to determine row length. Length N-J+1
596 * of the rows to swap back doesn't include diagonal element)
597  j = j + 1
598  IF( jp.NE.jj .AND. j.LE.n )
599  $ CALL cswap( n-j+1, a( jp, j ), lda, a( jj, j ), lda )
600  IF( j.LE.n )
601  $ GO TO 60
602 *
603 * Set KB to the number of columns factorized
604 *
605  kb = n - k
606 *
607  ELSE
608 *
609 * Factorize the leading columns of A using the lower triangle
610 * of A and working forwards, and compute the matrix W = L21*D
611 * for use in updating A22 (note that conjg(W) is actually stored)
612 *
613 * K is the main loop index, increasing from 1 in steps of 1 or 2
614 *
615  k = 1
616  70 CONTINUE
617 *
618 * Exit from loop
619 *
620  IF( ( k.GE.nb .AND. nb.LT.n ) .OR. k.GT.n )
621  $ GO TO 90
622 *
623  kstep = 1
624 *
625 * Copy column K of A to column K of W and update it
626 *
627  w( k, k ) = REAL( A( K, K ) )
628  IF( k.LT.n )
629  $ CALL ccopy( n-k, a( k+1, k ), 1, w( k+1, k ), 1 )
630  CALL cgemv( 'No transpose', n-k+1, k-1, -cone, a( k, 1 ), lda,
631  $ w( k, 1 ), ldw, cone, w( k, k ), 1 )
632  w( k, k ) = REAL( W( K, K ) )
633 *
634 * Determine rows and columns to be interchanged and whether
635 * a 1-by-1 or 2-by-2 pivot block will be used
636 *
637  absakk = abs( REAL( W( K, K ) ) )
638 *
639 * IMAX is the row-index of the largest off-diagonal element in
640 * column K, and COLMAX is its absolute value.
641 * Determine both COLMAX and IMAX.
642 *
643  IF( k.LT.n ) THEN
644  imax = k + icamax( n-k, w( k+1, k ), 1 )
645  colmax = cabs1( w( imax, k ) )
646  ELSE
647  colmax = zero
648  END IF
649 *
650  IF( max( absakk, colmax ).EQ.zero ) THEN
651 *
652 * Column K is zero or underflow: set INFO and continue
653 *
654  IF( info.EQ.0 )
655  $ info = k
656  kp = k
657  a( k, k ) = REAL( A( K, K ) )
658  ELSE
659 *
660 * ============================================================
661 *
662 * BEGIN pivot search
663 *
664 * Case(1)
665  IF( absakk.GE.alpha*colmax ) THEN
666 *
667 * no interchange, use 1-by-1 pivot block
668 *
669  kp = k
670  ELSE
671 *
672 * BEGIN pivot search along IMAX row
673 *
674 *
675 * Copy column IMAX to column K+1 of W and update it
676 *
677  CALL ccopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1 )
678  CALL clacgv( imax-k, w( k, k+1 ), 1 )
679  w( imax, k+1 ) = REAL( A( IMAX, IMAX ) )
680  IF( imax.LT.n )
681  $ CALL ccopy( n-imax, a( imax+1, imax ), 1,
682  $ w( imax+1, k+1 ), 1 )
683  CALL cgemv( 'No transpose', n-k+1, k-1, -cone, a( k, 1 ),
684  $ lda, w( imax, 1 ), ldw, cone, w( k, k+1 ),
685  $ 1 )
686  w( imax, k+1 ) = REAL( W( IMAX, K+1 ) )
687 *
688 * JMAX is the column-index of the largest off-diagonal
689 * element in row IMAX, and ROWMAX is its absolute value.
690 * Determine only ROWMAX.
691 *
692  jmax = k - 1 + icamax( imax-k, w( k, k+1 ), 1 )
693  rowmax = cabs1( w( jmax, k+1 ) )
694  IF( imax.LT.n ) THEN
695  jmax = imax + icamax( n-imax, w( imax+1, k+1 ), 1 )
696  rowmax = max( rowmax, cabs1( w( jmax, k+1 ) ) )
697  END IF
698 *
699 * Case(2)
700  IF( absakk.GE.alpha*colmax*( colmax / rowmax ) ) THEN
701 *
702 * no interchange, use 1-by-1 pivot block
703 *
704  kp = k
705 *
706 * Case(3)
707  ELSE IF( abs( REAL( W( IMAX, K+1 ) ) ).GE.alpha*rowmax )
708  $ THEN
709 *
710 * interchange rows and columns K and IMAX, use 1-by-1
711 * pivot block
712 *
713  kp = imax
714 *
715 * copy column K+1 of W to column K of W
716 *
717  CALL ccopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
718 *
719 * Case(4)
720  ELSE
721 *
722 * interchange rows and columns K+1 and IMAX, use 2-by-2
723 * pivot block
724 *
725  kp = imax
726  kstep = 2
727  END IF
728 *
729 *
730 * END pivot search along IMAX row
731 *
732  END IF
733 *
734 * END pivot search
735 *
736 * ============================================================
737 *
738 * KK is the column of A where pivoting step stopped
739 *
740  kk = k + kstep - 1
741 *
742 * Interchange rows and columns KP and KK.
743 * Updated column KP is already stored in column KK of W.
744 *
745  IF( kp.NE.kk ) THEN
746 *
747 * Copy non-updated column KK to column KP of submatrix A
748 * at step K. No need to copy element into column K
749 * (or K and K+1 for 2-by-2 pivot) of A, since these columns
750 * will be later overwritten.
751 *
752  a( kp, kp ) = REAL( A( KK, KK ) )
753  CALL ccopy( kp-kk-1, a( kk+1, kk ), 1, a( kp, kk+1 ),
754  $ lda )
755  CALL clacgv( kp-kk-1, a( kp, kk+1 ), lda )
756  IF( kp.LT.n )
757  $ CALL ccopy( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
758 *
759 * Interchange rows KK and KP in first K-1 columns of A
760 * (columns K (or K and K+1 for 2-by-2 pivot) of A will be
761 * later overwritten). Interchange rows KK and KP
762 * in first KK columns of W.
763 *
764  IF( k.GT.1 )
765  $ CALL cswap( k-1, a( kk, 1 ), lda, a( kp, 1 ), lda )
766  CALL cswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
767  END IF
768 *
769  IF( kstep.EQ.1 ) THEN
770 *
771 * 1-by-1 pivot block D(k): column k of W now holds
772 *
773 * W(k) = L(k)*D(k),
774 *
775 * where L(k) is the k-th column of L
776 *
777 * (1) Store subdiag. elements of column L(k)
778 * and 1-by-1 block D(k) in column k of A.
779 * (NOTE: Diagonal element L(k,k) is a UNIT element
780 * and not stored)
781 * A(k,k) := D(k,k) = W(k,k)
782 * A(k+1:N,k) := L(k+1:N,k) = W(k+1:N,k)/D(k,k)
783 *
784 * (NOTE: No need to use for Hermitian matrix
785 * A( K, K ) = DBLE( W( K, K) ) to separately copy diagonal
786 * element D(k,k) from W (potentially saves only one load))
787  CALL ccopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
788  IF( k.LT.n ) THEN
789 *
790 * (NOTE: No need to check if A(k,k) is NOT ZERO,
791 * since that was ensured earlier in pivot search:
792 * case A(k,k) = 0 falls into 2x2 pivot case(4))
793 *
794  r1 = one / REAL( A( K, K ) )
795  CALL csscal( n-k, r1, a( k+1, k ), 1 )
796 *
797 * (2) Conjugate column W(k)
798 *
799  CALL clacgv( n-k, w( k+1, k ), 1 )
800  END IF
801 *
802  ELSE
803 *
804 * 2-by-2 pivot block D(k): columns k and k+1 of W now hold
805 *
806 * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
807 *
808 * where L(k) and L(k+1) are the k-th and (k+1)-th columns
809 * of L
810 *
811 * (1) Store L(k+2:N,k) and L(k+2:N,k+1) and 2-by-2
812 * block D(k:k+1,k:k+1) in columns k and k+1 of A.
813 * (NOTE: 2-by-2 diagonal block L(k:k+1,k:k+1) is a UNIT
814 * block and not stored)
815 * A(k:k+1,k:k+1) := D(k:k+1,k:k+1) = W(k:k+1,k:k+1)
816 * A(k+2:N,k:k+1) := L(k+2:N,k:k+1) =
817 * = W(k+2:N,k:k+1) * ( D(k:k+1,k:k+1)**(-1) )
818 *
819  IF( k.LT.n-1 ) THEN
820 *
821 * Factor out the columns of the inverse of 2-by-2 pivot
822 * block D, so that each column contains 1, to reduce the
823 * number of FLOPS when we multiply panel
824 * ( W(kw-1) W(kw) ) by this inverse, i.e. by D**(-1).
825 *
826 * D**(-1) = ( d11 cj(d21) )**(-1) =
827 * ( d21 d22 )
828 *
829 * = 1/(d11*d22-|d21|**2) * ( ( d22) (-cj(d21) ) ) =
830 * ( (-d21) ( d11 ) )
831 *
832 * = 1/(|d21|**2) * 1/((d11/cj(d21))*(d22/d21)-1) *
833 *
834 * * ( d21*( d22/d21 ) conj(d21)*( - 1 ) ) =
835 * ( ( -1 ) ( d11/conj(d21) ) )
836 *
837 * = 1/(|d21|**2) * 1/(D22*D11-1) *
838 *
839 * * ( d21*( D11 ) conj(d21)*( -1 ) ) =
840 * ( ( -1 ) ( D22 ) )
841 *
842 * = (1/|d21|**2) * T * ( d21*( D11 ) conj(d21)*( -1 ) ) =
843 * ( ( -1 ) ( D22 ) )
844 *
845 * = ( (T/conj(d21))*( D11 ) (T/d21)*( -1 ) ) =
846 * ( ( -1 ) ( D22 ) )
847 *
848 * = ( conj(D21)*( D11 ) D21*( -1 ) )
849 * ( ( -1 ) ( D22 ) )
850 *
851 * where D11 = d22/d21,
852 * D22 = d11/conj(d21),
853 * D21 = T/d21,
854 * T = 1/(D22*D11-1).
855 *
856 * (NOTE: No need to check for division by ZERO,
857 * since that was ensured earlier in pivot search:
858 * (a) d21 != 0, since in 2x2 pivot case(4)
859 * |d21| should be larger than |d11| and |d22|;
860 * (b) (D22*D11 - 1) != 0, since from (a),
861 * both |D11| < 1, |D22| < 1, hence |D22*D11| << 1.)
862 *
863  d21 = w( k+1, k )
864  d11 = w( k+1, k+1 ) / d21
865  d22 = w( k, k ) / conjg( d21 )
866  t = one / ( REAL( d11*d22 )-one )
867  d21 = t / d21
868 *
869 * Update elements in columns A(k) and A(k+1) as
870 * dot products of rows of ( W(k) W(k+1) ) and columns
871 * of D**(-1)
872 *
873  DO 80 j = k + 2, n
874  a( j, k ) = conjg( d21 )*
875  $ ( d11*w( j, k )-w( j, k+1 ) )
876  a( j, k+1 ) = d21*( d22*w( j, k+1 )-w( j, k ) )
877  80 CONTINUE
878  END IF
879 *
880 * Copy D(k) to A
881 *
882  a( k, k ) = w( k, k )
883  a( k+1, k ) = w( k+1, k )
884  a( k+1, k+1 ) = w( k+1, k+1 )
885 *
886 * (2) Conjugate columns W(k) and W(k+1)
887 *
888  CALL clacgv( n-k, w( k+1, k ), 1 )
889  CALL clacgv( n-k-1, w( k+2, k+1 ), 1 )
890 *
891  END IF
892 *
893  END IF
894 *
895 * Store details of the interchanges in IPIV
896 *
897  IF( kstep.EQ.1 ) THEN
898  ipiv( k ) = kp
899  ELSE
900  ipiv( k ) = -kp
901  ipiv( k+1 ) = -kp
902  END IF
903 *
904 * Increase K and return to the start of the main loop
905 *
906  k = k + kstep
907  GO TO 70
908 *
909  90 CONTINUE
910 *
911 * Update the lower triangle of A22 (= A(k:n,k:n)) as
912 *
913 * A22 := A22 - L21*D*L21**H = A22 - L21*W**H
914 *
915 * computing blocks of NB columns at a time (note that conjg(W) is
916 * actually stored)
917 *
918  DO 110 j = k, n, nb
919  jb = min( nb, n-j+1 )
920 *
921 * Update the lower triangle of the diagonal block
922 *
923  DO 100 jj = j, j + jb - 1
924  a( jj, jj ) = REAL( A( JJ, JJ ) )
925  CALL cgemv( 'No transpose', j+jb-jj, k-1, -cone,
926  $ a( jj, 1 ), lda, w( jj, 1 ), ldw, cone,
927  $ a( jj, jj ), 1 )
928  a( jj, jj ) = REAL( A( JJ, JJ ) )
929  100 CONTINUE
930 *
931 * Update the rectangular subdiagonal block
932 *
933  IF( j+jb.LE.n )
934  $ CALL cgemm( 'No transpose', 'Transpose', n-j-jb+1, jb,
935  $ k-1, -cone, a( j+jb, 1 ), lda, w( j, 1 ),
936  $ ldw, cone, a( j+jb, j ), lda )
937  110 CONTINUE
938 *
939 * Put L21 in standard form by partially undoing the interchanges
940 * of rows in columns 1:k-1 looping backwards from k-1 to 1
941 *
942  j = k - 1
943  120 CONTINUE
944 *
945 * Undo the interchanges (if any) of rows J and JP
946 * at each step J
947 *
948 * (Here, J is a diagonal index)
949  jj = j
950  jp = ipiv( j )
951  IF( jp.LT.0 ) THEN
952  jp = -jp
953 * (Here, J is a diagonal index)
954  j = j - 1
955  END IF
956 * (NOTE: Here, J is used to determine row length. Length J
957 * of the rows to swap back doesn't include diagonal element)
958  j = j - 1
959  IF( jp.NE.jj .AND. j.GE.1 )
960  $ CALL cswap( j, a( jp, 1 ), lda, a( jj, 1 ), lda )
961  IF( j.GE.1 )
962  $ GO TO 120
963 *
964 * Set KB to the number of columns factorized
965 *
966  kb = k - 1
967 *
968  END IF
969  RETURN
970 *
971 * End of CLAHEF
972 *
cgemv
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:160
icamax
integer function icamax(N, CX, INCX)
ICAMAX
Definition: icamax.f:53
ccopy
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:52
clacgv
subroutine clacgv(N, X, INCX)
CLACGV conjugates a complex vector.
Definition: clacgv.f:76
cswap
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:52
cgemm
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:189
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
csscal
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:54

Here is the call graph for this function:

Here is the caller graph for this function: