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

◆ dsytf2_rk()

subroutine dsytf2_rk ( character  uplo,
integer  n,
double precision, dimension( lda, * )  a,
integer  lda,
double precision, dimension( * )  e,
integer, dimension( * )  ipiv,
integer  info 
)

DSYTF2_RK computes the factorization of a real symmetric indefinite matrix using the bounded Bunch-Kaufman (rook) diagonal pivoting method (BLAS2 unblocked algorithm).

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

Purpose:
 DSYTF2_RK computes the factorization of a real symmetric matrix A
 using the bounded Bunch-Kaufman (rook) diagonal pivoting method:

    A = P*U*D*(U**T)*(P**T) or A = P*L*D*(L**T)*(P**T),

 where U (or L) is unit upper (or lower) triangular matrix,
 U**T (or L**T) is the transpose of U (or L), P is a permutation
 matrix, P**T is the transpose of P, and D is symmetric and block
 diagonal with 1-by-1 and 2-by-2 diagonal blocks.

 This is the unblocked version of the algorithm, calling Level 2 BLAS.
 For more information see Further Details section.
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,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. For more info see Further
          Details section.

          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 matrix A(1:N,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,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 matrix A(1:N,1:N).
                  If -IPIV(k-1) = k-1, no interchange occurred.

            c) In both cases a) and b), 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 matrix A(1:N,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: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 matrix A(1:N,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 matrix A(1:N,1:N).
                  If -IPIV(k+1) = k+1, no interchange occurred.

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

            d) NOTE: Any entry IPIV(k) is always NONZERO on output.
[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.
Further Details:
 TODO: put further details
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

  01-01-96 - Based on modifications by
    J. Lewis, Boeing Computer Services Company
    A. Petitet, Computer Science Dept.,
                Univ. of Tenn., Knoxville abd , USA

Definition at line 240 of file dsytf2_rk.f.

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