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

◆ 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.
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 260 of file clasyf_rk.f.

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