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

◆ dlasyf_rk()

subroutine dlasyf_rk ( character  uplo,
integer  n,
integer  nb,
integer  kb,
double precision, dimension( lda, * )  a,
integer  lda,
double precision, dimension( * )  e,
integer, dimension( * )  ipiv,
double precision, dimension( ldw, * )  w,
integer  ldw,
integer  info 
)

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

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

Purpose:
 DLASYF_RK computes a partial factorization of a real 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.

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