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

◆ 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 238 of file ssytf2_rk.f.

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