LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ 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.
Date
December 2016
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 243 of file dsytf2_rk.f.

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