LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ ssytf2_rk()

subroutine ssytf2_rk ( character  UPLO,
integer  N,
real, dimension( lda, * )  A,
integer  LDA,
real, dimension( * )  E,
integer, dimension( * )  IPIV,
integer  INFO 
)

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

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

Purpose:
 SSYTF2_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 REAL 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 REAL 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 ssytf2_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  REAL A( LDA, * ), E( * )
253 * ..
254 *
255 * =====================================================================
256 *
257 * .. Parameters ..
258  REAL ZERO, ONE
259  parameter( zero = 0.0e+0, one = 1.0e+0 )
260  REAL EIGHT, SEVTEN
261  parameter( eight = 8.0e+0, sevten = 17.0e+0 )
262 * ..
263 * .. Local Scalars ..
264  LOGICAL UPPER, DONE
265  INTEGER I, IMAX, J, JMAX, ITEMP, K, KK, KP, KSTEP,
266  $ P, II
267  REAL ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22,
268  $ ROWMAX, STEMP, T, WK, WKM1, WKP1, SFMIN
269 * ..
270 * .. External Functions ..
271  LOGICAL LSAME
272  INTEGER ISAMAX
273  REAL SLAMCH
274  EXTERNAL lsame, isamax, slamch
275 * ..
276 * .. External Subroutines ..
277  EXTERNAL sscal, sswap, ssyr, 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( 'SSYTF2_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 = slamch( '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 = isamax( 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 + isamax( 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 = isamax( imax-1, a( 1, imax ), 1 )
396  stemp = abs( a( itemp, imax ) )
397  IF( stemp.GT.rowmax ) THEN
398  rowmax = stemp
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 sswap( p-1, a( 1, k ), 1, a( 1, p ), 1 )
452  IF( p.LT.(k-1) )
453  $ CALL sswap( 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 sswap( 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 sswap( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
477  IF( ( kk.GT.1 ) .AND. ( kp.LT.(kk-1) ) )
478  $ CALL sswap( 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 sswap( 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 ssyr( uplo, k-1, -d11, a( 1, k ), 1, a, lda )
521 *
522 * Store U(k) in column k
523 *
524  CALL sscal( 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 ssyr( 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 + isamax( 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 + isamax( 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 + isamax( n-imax, a( imax+1, imax ),
705  $ 1 )
706  stemp = abs( a( itemp, imax ) )
707  IF( stemp.GT.rowmax ) THEN
708  rowmax = stemp
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 sswap( n-p, a( p+1, k ), 1, a( p+1, p ), 1 )
762  IF( p.GT.(k+1) )
763  $ CALL sswap( 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 sswap( 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 sswap( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
786  IF( ( kk.LT.n ) .AND. ( kp.GT.(kk+1) ) )
787  $ CALL sswap( 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 sswap( 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 ssyr( 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 sscal( 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 ssyr( 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 SSYTF2_RK
939 *
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:71
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:82
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
subroutine ssyr(UPLO, N, ALPHA, X, INCX, A, LDA)
SSYR
Definition: ssyr.f:132
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: