LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ clasyf_rk()

subroutine clasyf_rk ( character  UPLO,
integer  N,
integer  NB,
integer  KB,
complex, dimension( lda, * )  A,
integer  LDA,
complex, dimension( * )  E,
integer, dimension( * )  IPIV,
complex, dimension( ldw, * )  W,
integer  LDW,
integer  INFO 
)

CLASYF_RK computes a partial factorization of a complex symmetric indefinite matrix using bounded Bunch-Kaufman (rook) diagonal pivoting method.

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

Purpose:
 CLASYF_RK computes a partial factorization of a complex symmetric
 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**T U22**T )

 A  =  ( L11  0 ) (  D   0  ) ( L11**T L21**T )  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.

 CLASYF_RK is an auxiliary routine called by CSYTRF_RK. 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
          symmetric 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 symmetric 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, contains:
            a) ONLY diagonal elements of the symmetric block diagonal
               matrix D on the diagonal of A, i.e. D(k,k) = A(k,k);
               (superdiagonal (or subdiagonal) elements of D
                are stored on exit in array E), and
            b) If UPLO = 'U': factor U in the superdiagonal part of A.
               If UPLO = 'L': factor L in the subdiagonal part of A.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[out]E
          E is COMPLEX array, dimension (N)
          On exit, contains the superdiagonal (or subdiagonal)
          elements of the symmetric block diagonal matrix D
          with 1-by-1 or 2-by-2 diagonal blocks, where
          If UPLO = 'U': E(i) = D(i-1,i), i=2:N, E(1) is set to 0;
          If UPLO = 'L': E(i) = D(i+1,i), i=1:N-1, E(N) is set to 0.

          NOTE: For 1-by-1 diagonal block D(k), where
          1 <= k <= N, the element E(k) is set to 0 in both
          UPLO = 'U' or UPLO = 'L' cases.
[out]IPIV
          IPIV is INTEGER array, dimension (N)
          IPIV describes the permutation matrix P in the factorization
          of matrix A as follows. The absolute value of IPIV(k)
          represents the index of row and column that were
          interchanged with the k-th row and column. The value of UPLO
          describes the order in which the interchanges were applied.
          Also, the sign of IPIV represents the block structure of
          the symmetric block diagonal matrix D with 1-by-1 or 2-by-2
          diagonal blocks which correspond to 1 or 2 interchanges
          at each factorization step.

          If UPLO = 'U',
          ( in factorization order, k decreases from N to 1 ):
            a) A single positive entry IPIV(k) > 0 means:
               D(k,k) is a 1-by-1 diagonal block.
               If IPIV(k) != k, rows and columns k and IPIV(k) were
               interchanged in the submatrix A(1:N,N-KB+1:N);
               If IPIV(k) = k, no interchange occurred.


            b) A pair of consecutive negative entries
               IPIV(k) < 0 and IPIV(k-1) < 0 means:
               D(k-1:k,k-1:k) is a 2-by-2 diagonal block.
               (NOTE: negative entries in IPIV appear ONLY in pairs).
               1) If -IPIV(k) != k, rows and columns
                  k and -IPIV(k) were interchanged
                  in the matrix A(1:N,N-KB+1:N).
                  If -IPIV(k) = k, no interchange occurred.
               2) If -IPIV(k-1) != k-1, rows and columns
                  k-1 and -IPIV(k-1) were interchanged
                  in the submatrix A(1:N,N-KB+1:N).
                  If -IPIV(k-1) = k-1, no interchange occurred.

            c) In both cases a) and b) is always ABS( IPIV(k) ) <= k.

            d) NOTE: Any entry IPIV(k) is always NONZERO on output.

          If UPLO = 'L',
          ( in factorization order, k increases from 1 to N ):
            a) A single positive entry IPIV(k) > 0 means:
               D(k,k) is a 1-by-1 diagonal block.
               If IPIV(k) != k, rows and columns k and IPIV(k) were
               interchanged in the submatrix A(1:N,1:KB).
               If IPIV(k) = k, no interchange occurred.

            b) A pair of consecutive negative entries
               IPIV(k) < 0 and IPIV(k+1) < 0 means:
               D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
               (NOTE: negative entries in IPIV appear ONLY in pairs).
               1) If -IPIV(k) != k, rows and columns
                  k and -IPIV(k) were interchanged
                  in the submatrix A(1:N,1:KB).
                  If -IPIV(k) = k, no interchange occurred.
               2) If -IPIV(k+1) != k+1, rows and columns
                  k-1 and -IPIV(k-1) were interchanged
                  in the submatrix A(1:N,1:KB).
                  If -IPIV(k+1) = k+1, no interchange occurred.

            c) In both cases a) and b) is always ABS( IPIV(k) ) >= k.

            d) NOTE: Any entry IPIV(k) is always NONZERO on output.
[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, the k-th argument had an illegal value

          > 0: If INFO = k, the matrix A is singular, because:
                 If UPLO = 'U': column k in the upper
                 triangular part of A contains all zeros.
                 If UPLO = 'L': column k in the lower
                 triangular part of A contains all zeros.

               Therefore D(k,k) is exactly zero, and superdiagonal
               elements of column k of U (or subdiagonal elements of
               column k of L ) are all zeros. 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.

               NOTE: INFO only stores the first occurrence of
               a singularity, any subsequent occurrence of singularity
               is not stored in INFO even though the factorization
               always completes.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
December 2016
Contributors:
  December 2016,  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 264 of file clasyf_rk.f.

264 *
265 * -- LAPACK computational routine (version 3.7.0) --
266 * -- LAPACK is a software package provided by Univ. of Tennessee, --
267 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
268 * December 2016
269 *
270 * .. Scalar Arguments ..
271  CHARACTER uplo
272  INTEGER info, kb, lda, ldw, n, nb
273 * ..
274 * .. Array Arguments ..
275  INTEGER ipiv( * )
276  COMPLEX a( lda, * ), e( * ), w( ldw, * )
277 * ..
278 *
279 * =====================================================================
280 *
281 * .. Parameters ..
282  REAL zero, one
283  parameter( zero = 0.0e+0, one = 1.0e+0 )
284  REAL eight, sevten
285  parameter( eight = 8.0e+0, sevten = 17.0e+0 )
286  COMPLEX cone, czero
287  parameter( cone = ( 1.0e+0, 0.0e+0 ),
288  $ czero = ( 0.0e+0, 0.0e+0 ) )
289 * ..
290 * .. Local Scalars ..
291  LOGICAL done
292  INTEGER imax, itemp, j, jb, jj, jmax, k, kk, kw, kkw,
293  $ kp, kstep, p, ii
294  REAL absakk, alpha, colmax, rowmax, sfmin, stemp
295  COMPLEX d11, d12, d21, d22, r1, t, z
296 * ..
297 * .. External Functions ..
298  LOGICAL lsame
299  INTEGER icamax
300  REAL slamch
301  EXTERNAL lsame, icamax, slamch
302 * ..
303 * .. External Subroutines ..
304  EXTERNAL ccopy, cgemm, cgemv, cscal, cswap
305 * ..
306 * .. Intrinsic Functions ..
307  INTRINSIC abs, aimag, max, min, REAL, sqrt
308 * ..
309 * .. Statement Functions ..
310  REAL cabs1
311 * ..
312 * .. Statement Function definitions ..
313  cabs1( z ) = abs( REAL( Z ) ) + abs( aimag( z ) )
314 * ..
315 * .. Executable Statements ..
316 *
317  info = 0
318 *
319 * Initialize ALPHA for use in choosing pivot block size.
320 *
321  alpha = ( one+sqrt( sevten ) ) / eight
322 *
323 * Compute machine safe minimum
324 *
325  sfmin = slamch( 'S' )
326 *
327  IF( lsame( uplo, 'U' ) ) THEN
328 *
329 * Factorize the trailing columns of A using the upper triangle
330 * of A and working backwards, and compute the matrix W = U12*D
331 * for use in updating A11
332 *
333 * Initilize the first entry of array E, where superdiagonal
334 * elements of D are stored
335 *
336  e( 1 ) = czero
337 *
338 * K is the main loop index, decreasing from N in steps of 1 or 2
339 *
340  k = n
341  10 CONTINUE
342 *
343 * KW is the column of W which corresponds to column K of A
344 *
345  kw = nb + k - n
346 *
347 * Exit from loop
348 *
349  IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
350  $ GO TO 30
351 *
352  kstep = 1
353  p = k
354 *
355 * Copy column K of A to column KW of W and update it
356 *
357  CALL ccopy( k, a( 1, k ), 1, w( 1, kw ), 1 )
358  IF( k.LT.n )
359  $ CALL cgemv( 'No transpose', k, n-k, -cone, a( 1, k+1 ),
360  $ lda, w( k, kw+1 ), ldw, cone, w( 1, kw ), 1 )
361 *
362 * Determine rows and columns to be interchanged and whether
363 * a 1-by-1 or 2-by-2 pivot block will be used
364 *
365  absakk = cabs1( w( k, kw ) )
366 *
367 * IMAX is the row-index of the largest off-diagonal element in
368 * column K, and COLMAX is its absolute value.
369 * Determine both COLMAX and IMAX.
370 *
371  IF( k.GT.1 ) THEN
372  imax = icamax( k-1, w( 1, kw ), 1 )
373  colmax = cabs1( w( imax, kw ) )
374  ELSE
375  colmax = zero
376  END IF
377 *
378  IF( max( absakk, colmax ).EQ.zero ) THEN
379 *
380 * Column K is zero or underflow: set INFO and continue
381 *
382  IF( info.EQ.0 )
383  $ info = k
384  kp = k
385  CALL ccopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
386 *
387 * Set E( K ) to zero
388 *
389  IF( k.GT.1 )
390  $ e( k ) = czero
391 *
392  ELSE
393 *
394 * ============================================================
395 *
396 * Test for interchange
397 *
398 * Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
399 * (used to handle NaN and Inf)
400 *
401  IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
402 *
403 * no interchange, use 1-by-1 pivot block
404 *
405  kp = k
406 *
407  ELSE
408 *
409  done = .false.
410 *
411 * Loop until pivot found
412 *
413  12 CONTINUE
414 *
415 * Begin pivot search loop body
416 *
417 *
418 * Copy column IMAX to column KW-1 of W and update it
419 *
420  CALL ccopy( imax, a( 1, imax ), 1, w( 1, kw-1 ), 1 )
421  CALL ccopy( k-imax, a( imax, imax+1 ), lda,
422  $ w( imax+1, kw-1 ), 1 )
423 *
424  IF( k.LT.n )
425  $ CALL cgemv( 'No transpose', k, n-k, -cone,
426  $ a( 1, k+1 ), lda, w( imax, kw+1 ), ldw,
427  $ cone, w( 1, kw-1 ), 1 )
428 *
429 * JMAX is the column-index of the largest off-diagonal
430 * element in row IMAX, and ROWMAX is its absolute value.
431 * Determine both ROWMAX and JMAX.
432 *
433  IF( imax.NE.k ) THEN
434  jmax = imax + icamax( k-imax, w( imax+1, kw-1 ),
435  $ 1 )
436  rowmax = cabs1( w( jmax, kw-1 ) )
437  ELSE
438  rowmax = zero
439  END IF
440 *
441  IF( imax.GT.1 ) THEN
442  itemp = icamax( imax-1, w( 1, kw-1 ), 1 )
443  stemp = cabs1( w( itemp, kw-1 ) )
444  IF( stemp.GT.rowmax ) THEN
445  rowmax = stemp
446  jmax = itemp
447  END IF
448  END IF
449 *
450 * Equivalent to testing for
451 * CABS1( W( IMAX, KW-1 ) ).GE.ALPHA*ROWMAX
452 * (used to handle NaN and Inf)
453 *
454  IF( .NOT.(cabs1( w( imax, kw-1 ) ).LT.alpha*rowmax ) )
455  $ THEN
456 *
457 * interchange rows and columns K and IMAX,
458 * use 1-by-1 pivot block
459 *
460  kp = imax
461 *
462 * copy column KW-1 of W to column KW of W
463 *
464  CALL ccopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
465 *
466  done = .true.
467 *
468 * Equivalent to testing for ROWMAX.EQ.COLMAX,
469 * (used to handle NaN and Inf)
470 *
471  ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
472  $ THEN
473 *
474 * interchange rows and columns K-1 and IMAX,
475 * use 2-by-2 pivot block
476 *
477  kp = imax
478  kstep = 2
479  done = .true.
480  ELSE
481 *
482 * Pivot not found: set params and repeat
483 *
484  p = imax
485  colmax = rowmax
486  imax = jmax
487 *
488 * Copy updated JMAXth (next IMAXth) column to Kth of W
489 *
490  CALL ccopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
491 *
492  END IF
493 *
494 * End pivot search loop body
495 *
496  IF( .NOT. done ) GOTO 12
497 *
498  END IF
499 *
500 * ============================================================
501 *
502  kk = k - kstep + 1
503 *
504 * KKW is the column of W which corresponds to column KK of A
505 *
506  kkw = nb + kk - n
507 *
508  IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
509 *
510 * Copy non-updated column K to column P
511 *
512  CALL ccopy( k-p, a( p+1, k ), 1, a( p, p+1 ), lda )
513  CALL ccopy( p, a( 1, k ), 1, a( 1, p ), 1 )
514 *
515 * Interchange rows K and P in last N-K+1 columns of A
516 * and last N-K+2 columns of W
517 *
518  CALL cswap( n-k+1, a( k, k ), lda, a( p, k ), lda )
519  CALL cswap( n-kk+1, w( k, kkw ), ldw, w( p, kkw ), ldw )
520  END IF
521 *
522 * Updated column KP is already stored in column KKW of W
523 *
524  IF( kp.NE.kk ) THEN
525 *
526 * Copy non-updated column KK to column KP
527 *
528  a( kp, k ) = a( kk, k )
529  CALL ccopy( k-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
530  $ lda )
531  CALL ccopy( kp, a( 1, kk ), 1, a( 1, kp ), 1 )
532 *
533 * Interchange rows KK and KP in last N-KK+1 columns
534 * of A and W
535 *
536  CALL cswap( n-kk+1, a( kk, kk ), lda, a( kp, kk ), lda )
537  CALL cswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
538  $ ldw )
539  END IF
540 *
541  IF( kstep.EQ.1 ) THEN
542 *
543 * 1-by-1 pivot block D(k): column KW of W now holds
544 *
545 * W(k) = U(k)*D(k)
546 *
547 * where U(k) is the k-th column of U
548 *
549 * Store U(k) in column k of A
550 *
551  CALL ccopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
552  IF( k.GT.1 ) THEN
553  IF( cabs1( a( k, k ) ).GE.sfmin ) THEN
554  r1 = cone / a( k, k )
555  CALL cscal( k-1, r1, a( 1, k ), 1 )
556  ELSE IF( a( k, k ).NE.czero ) THEN
557  DO 14 ii = 1, k - 1
558  a( ii, k ) = a( ii, k ) / a( k, k )
559  14 CONTINUE
560  END IF
561 *
562 * Store the superdiagonal element of D in array E
563 *
564  e( k ) = czero
565 *
566  END IF
567 *
568  ELSE
569 *
570 * 2-by-2 pivot block D(k): columns KW and KW-1 of W now
571 * hold
572 *
573 * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
574 *
575 * where U(k) and U(k-1) are the k-th and (k-1)-th columns
576 * of U
577 *
578  IF( k.GT.2 ) THEN
579 *
580 * Store U(k) and U(k-1) in columns k and k-1 of A
581 *
582  d12 = w( k-1, kw )
583  d11 = w( k, kw ) / d12
584  d22 = w( k-1, kw-1 ) / d12
585  t = cone / ( d11*d22-cone )
586  DO 20 j = 1, k - 2
587  a( j, k-1 ) = t*( (d11*w( j, kw-1 )-w( j, kw ) ) /
588  $ d12 )
589  a( j, k ) = t*( ( d22*w( j, kw )-w( j, kw-1 ) ) /
590  $ d12 )
591  20 CONTINUE
592  END IF
593 *
594 * Copy diagonal elements of D(K) to A,
595 * copy superdiagonal element of D(K) to E(K) and
596 * ZERO out superdiagonal entry of A
597 *
598  a( k-1, k-1 ) = w( k-1, kw-1 )
599  a( k-1, k ) = czero
600  a( k, k ) = w( k, kw )
601  e( k ) = w( k-1, kw )
602  e( k-1 ) = czero
603 *
604  END IF
605 *
606 * End column K is nonsingular
607 *
608  END IF
609 *
610 * Store details of the interchanges in IPIV
611 *
612  IF( kstep.EQ.1 ) THEN
613  ipiv( k ) = kp
614  ELSE
615  ipiv( k ) = -p
616  ipiv( k-1 ) = -kp
617  END IF
618 *
619 * Decrease K and return to the start of the main loop
620 *
621  k = k - kstep
622  GO TO 10
623 *
624  30 CONTINUE
625 *
626 * Update the upper triangle of A11 (= A(1:k,1:k)) as
627 *
628 * A11 := A11 - U12*D*U12**T = A11 - U12*W**T
629 *
630 * computing blocks of NB columns at a time
631 *
632  DO 50 j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
633  jb = min( nb, k-j+1 )
634 *
635 * Update the upper triangle of the diagonal block
636 *
637  DO 40 jj = j, j + jb - 1
638  CALL cgemv( 'No transpose', jj-j+1, n-k, -cone,
639  $ a( j, k+1 ), lda, w( jj, kw+1 ), ldw, cone,
640  $ a( j, jj ), 1 )
641  40 CONTINUE
642 *
643 * Update the rectangular superdiagonal block
644 *
645  IF( j.GE.2 )
646  $ CALL cgemm( 'No transpose', 'Transpose', j-1, jb,
647  $ n-k, -cone, a( 1, k+1 ), lda, w( j, kw+1 ),
648  $ ldw, cone, a( 1, j ), lda )
649  50 CONTINUE
650 *
651 * Set KB to the number of columns factorized
652 *
653  kb = n - k
654 *
655  ELSE
656 *
657 * Factorize the leading columns of A using the lower triangle
658 * of A and working forwards, and compute the matrix W = L21*D
659 * for use in updating A22
660 *
661 * Initilize the unused last entry of the subdiagonal array E.
662 *
663  e( n ) = czero
664 *
665 * K is the main loop index, increasing from 1 in steps of 1 or 2
666 *
667  k = 1
668  70 CONTINUE
669 *
670 * Exit from loop
671 *
672  IF( ( k.GE.nb .AND. nb.LT.n ) .OR. k.GT.n )
673  $ GO TO 90
674 *
675  kstep = 1
676  p = k
677 *
678 * Copy column K of A to column K of W and update it
679 *
680  CALL ccopy( n-k+1, a( k, k ), 1, w( k, k ), 1 )
681  IF( k.GT.1 )
682  $ CALL cgemv( 'No transpose', n-k+1, k-1, -cone, a( k, 1 ),
683  $ lda, w( k, 1 ), ldw, cone, w( k, k ), 1 )
684 *
685 * Determine rows and columns to be interchanged and whether
686 * a 1-by-1 or 2-by-2 pivot block will be used
687 *
688  absakk = cabs1( w( k, k ) )
689 *
690 * IMAX is the row-index of the largest off-diagonal element in
691 * column K, and COLMAX is its absolute value.
692 * Determine both COLMAX and IMAX.
693 *
694  IF( k.LT.n ) THEN
695  imax = k + icamax( n-k, w( k+1, k ), 1 )
696  colmax = cabs1( w( imax, k ) )
697  ELSE
698  colmax = zero
699  END IF
700 *
701  IF( max( absakk, colmax ).EQ.zero ) THEN
702 *
703 * Column K is zero or underflow: set INFO and continue
704 *
705  IF( info.EQ.0 )
706  $ info = k
707  kp = k
708  CALL ccopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
709 *
710 * Set E( K ) to zero
711 *
712  IF( k.LT.n )
713  $ e( k ) = czero
714 *
715  ELSE
716 *
717 * ============================================================
718 *
719 * Test for interchange
720 *
721 * Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
722 * (used to handle NaN and Inf)
723 *
724  IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
725 *
726 * no interchange, use 1-by-1 pivot block
727 *
728  kp = k
729 *
730  ELSE
731 *
732  done = .false.
733 *
734 * Loop until pivot found
735 *
736  72 CONTINUE
737 *
738 * Begin pivot search loop body
739 *
740 *
741 * Copy column IMAX to column K+1 of W and update it
742 *
743  CALL ccopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1)
744  CALL ccopy( n-imax+1, a( imax, imax ), 1,
745  $ w( imax, k+1 ), 1 )
746  IF( k.GT.1 )
747  $ CALL cgemv( 'No transpose', n-k+1, k-1, -cone,
748  $ a( k, 1 ), lda, w( imax, 1 ), ldw,
749  $ cone, w( k, k+1 ), 1 )
750 *
751 * JMAX is the column-index of the largest off-diagonal
752 * element in row IMAX, and ROWMAX is its absolute value.
753 * Determine both ROWMAX and JMAX.
754 *
755  IF( imax.NE.k ) THEN
756  jmax = k - 1 + icamax( imax-k, w( k, k+1 ), 1 )
757  rowmax = cabs1( w( jmax, k+1 ) )
758  ELSE
759  rowmax = zero
760  END IF
761 *
762  IF( imax.LT.n ) THEN
763  itemp = imax + icamax( n-imax, w( imax+1, k+1 ), 1)
764  stemp = cabs1( w( itemp, k+1 ) )
765  IF( stemp.GT.rowmax ) THEN
766  rowmax = stemp
767  jmax = itemp
768  END IF
769  END IF
770 *
771 * Equivalent to testing for
772 * CABS1( W( IMAX, K+1 ) ).GE.ALPHA*ROWMAX
773 * (used to handle NaN and Inf)
774 *
775  IF( .NOT.( cabs1( w( imax, k+1 ) ).LT.alpha*rowmax ) )
776  $ THEN
777 *
778 * interchange rows and columns K and IMAX,
779 * use 1-by-1 pivot block
780 *
781  kp = imax
782 *
783 * copy column K+1 of W to column K of W
784 *
785  CALL ccopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
786 *
787  done = .true.
788 *
789 * Equivalent to testing for ROWMAX.EQ.COLMAX,
790 * (used to handle NaN and Inf)
791 *
792  ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
793  $ THEN
794 *
795 * interchange rows and columns K+1 and IMAX,
796 * use 2-by-2 pivot block
797 *
798  kp = imax
799  kstep = 2
800  done = .true.
801  ELSE
802 *
803 * Pivot not found: set params and repeat
804 *
805  p = imax
806  colmax = rowmax
807  imax = jmax
808 *
809 * Copy updated JMAXth (next IMAXth) column to Kth of W
810 *
811  CALL ccopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
812 *
813  END IF
814 *
815 * End pivot search loop body
816 *
817  IF( .NOT. done ) GOTO 72
818 *
819  END IF
820 *
821 * ============================================================
822 *
823  kk = k + kstep - 1
824 *
825  IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
826 *
827 * Copy non-updated column K to column P
828 *
829  CALL ccopy( p-k, a( k, k ), 1, a( p, k ), lda )
830  CALL ccopy( n-p+1, a( p, k ), 1, a( p, p ), 1 )
831 *
832 * Interchange rows K and P in first K columns of A
833 * and first K+1 columns of W
834 *
835  CALL cswap( k, a( k, 1 ), lda, a( p, 1 ), lda )
836  CALL cswap( kk, w( k, 1 ), ldw, w( p, 1 ), ldw )
837  END IF
838 *
839 * Updated column KP is already stored in column KK of W
840 *
841  IF( kp.NE.kk ) THEN
842 *
843 * Copy non-updated column KK to column KP
844 *
845  a( kp, k ) = a( kk, k )
846  CALL ccopy( kp-k-1, a( k+1, kk ), 1, a( kp, k+1 ), lda )
847  CALL ccopy( n-kp+1, a( kp, kk ), 1, a( kp, kp ), 1 )
848 *
849 * Interchange rows KK and KP in first KK columns of A and W
850 *
851  CALL cswap( kk, a( kk, 1 ), lda, a( kp, 1 ), lda )
852  CALL cswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
853  END IF
854 *
855  IF( kstep.EQ.1 ) THEN
856 *
857 * 1-by-1 pivot block D(k): column k of W now holds
858 *
859 * W(k) = L(k)*D(k)
860 *
861 * where L(k) is the k-th column of L
862 *
863 * Store L(k) in column k of A
864 *
865  CALL ccopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
866  IF( k.LT.n ) THEN
867  IF( cabs1( a( k, k ) ).GE.sfmin ) THEN
868  r1 = cone / a( k, k )
869  CALL cscal( n-k, r1, a( k+1, k ), 1 )
870  ELSE IF( a( k, k ).NE.czero ) THEN
871  DO 74 ii = k + 1, n
872  a( ii, k ) = a( ii, k ) / a( k, k )
873  74 CONTINUE
874  END IF
875 *
876 * Store the subdiagonal element of D in array E
877 *
878  e( k ) = czero
879 *
880  END IF
881 *
882  ELSE
883 *
884 * 2-by-2 pivot block D(k): columns k and k+1 of W now hold
885 *
886 * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
887 *
888 * where L(k) and L(k+1) are the k-th and (k+1)-th columns
889 * of L
890 *
891  IF( k.LT.n-1 ) THEN
892 *
893 * Store L(k) and L(k+1) in columns k and k+1 of A
894 *
895  d21 = w( k+1, k )
896  d11 = w( k+1, k+1 ) / d21
897  d22 = w( k, k ) / d21
898  t = cone / ( d11*d22-cone )
899  DO 80 j = k + 2, n
900  a( j, k ) = t*( ( d11*w( j, k )-w( j, k+1 ) ) /
901  $ d21 )
902  a( j, k+1 ) = t*( ( d22*w( j, k+1 )-w( j, k ) ) /
903  $ d21 )
904  80 CONTINUE
905  END IF
906 *
907 * Copy diagonal elements of D(K) to A,
908 * copy subdiagonal element of D(K) to E(K) and
909 * ZERO out subdiagonal entry of A
910 *
911  a( k, k ) = w( k, k )
912  a( k+1, k ) = czero
913  a( k+1, k+1 ) = w( k+1, k+1 )
914  e( k ) = w( k+1, k )
915  e( k+1 ) = czero
916 *
917  END IF
918 *
919 * End column K is nonsingular
920 *
921  END IF
922 *
923 * Store details of the interchanges in IPIV
924 *
925  IF( kstep.EQ.1 ) THEN
926  ipiv( k ) = kp
927  ELSE
928  ipiv( k ) = -p
929  ipiv( k+1 ) = -kp
930  END IF
931 *
932 * Increase K and return to the start of the main loop
933 *
934  k = k + kstep
935  GO TO 70
936 *
937  90 CONTINUE
938 *
939 * Update the lower triangle of A22 (= A(k:n,k:n)) as
940 *
941 * A22 := A22 - L21*D*L21**T = A22 - L21*W**T
942 *
943 * computing blocks of NB columns at a time
944 *
945  DO 110 j = k, n, nb
946  jb = min( nb, n-j+1 )
947 *
948 * Update the lower triangle of the diagonal block
949 *
950  DO 100 jj = j, j + jb - 1
951  CALL cgemv( 'No transpose', j+jb-jj, k-1, -cone,
952  $ a( jj, 1 ), lda, w( jj, 1 ), ldw, cone,
953  $ a( jj, jj ), 1 )
954  100 CONTINUE
955 *
956 * Update the rectangular subdiagonal block
957 *
958  IF( j+jb.LE.n )
959  $ CALL cgemm( 'No transpose', 'Transpose', n-j-jb+1, jb,
960  $ k-1, -cone, a( j+jb, 1 ), lda, w( j, 1 ),
961  $ ldw, cone, a( j+jb, j ), lda )
962  110 CONTINUE
963 *
964 * Set KB to the number of columns factorized
965 *
966  kb = k - 1
967 *
968  END IF
969 *
970  RETURN
971 *
972 * End of CLASYF_RK
973 *
integer function icamax(N, CX, INCX)
ICAMAX
Definition: icamax.f:73
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:80
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:160
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:83
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:83
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:189
Here is the call graph for this function:
Here is the caller graph for this function: