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

◆ ssytf2_rook()

subroutine ssytf2_rook ( character uplo,
integer n,
real, dimension( lda, * ) a,
integer lda,
integer, dimension( * ) ipiv,
integer info )

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

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

Purpose:
!>
!> SSYTF2_ROOK computes the factorization of a real symmetric matrix A
!> using the bounded Bunch-Kaufman () diagonal pivoting method:
!>
!>    A = U*D*U**T  or  A = L*D*L**T
!>
!> where U (or L) is a product of permutation and unit upper (lower)
!> triangular matrices, U**T is the transpose of U, 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.
!> 
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, the block diagonal matrix D and the multipliers used
!>          to obtain the factor U or L (see below for further details).
!> 
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!> 
[out]IPIV
!>          IPIV is INTEGER array, dimension (N)
!>          Details of the interchanges and the block structure of D.
!>
!>          If UPLO = 'U':
!>             If IPIV(k) > 0, then rows and columns k and IPIV(k)
!>             were interchanged and D(k,k) is a 1-by-1 diagonal block.
!>
!>             If IPIV(k) < 0 and IPIV(k-1) < 0, then rows and
!>             columns k and -IPIV(k) were interchanged and rows and
!>             columns k-1 and -IPIV(k-1) were inerchaged,
!>             D(k-1:k,k-1:k) is a 2-by-2 diagonal block.
!>
!>          If UPLO = 'L':
!>             If IPIV(k) > 0, then rows and columns k and IPIV(k)
!>             were interchanged and D(k,k) is a 1-by-1 diagonal block.
!>
!>             If IPIV(k) < 0 and IPIV(k+1) < 0, then rows and
!>             columns k and -IPIV(k) were interchanged and rows and
!>             columns k+1 and -IPIV(k+1) were inerchaged,
!>             D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
!> 
[out]INFO
!>          INFO is INTEGER
!>          = 0: successful exit
!>          < 0: if INFO = -k, the k-th argument had an illegal value
!>          > 0: if INFO = k, D(k,k) is exactly zero.  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.
!> 
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  If UPLO = 'U', then A = U*D*U**T, where
!>     U = P(n)*U(n)* ... *P(k)U(k)* ...,
!>  i.e., U is a product of terms P(k)*U(k), where k decreases from n to
!>  1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
!>  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
!>  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
!>  that if the diagonal block D(k) is of order s (s = 1 or 2), then
!>
!>             (   I    v    0   )   k-s
!>     U(k) =  (   0    I    0   )   s
!>             (   0    0    I   )   n-k
!>                k-s   s   n-k
!>
!>  If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
!>  If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
!>  and A(k,k), and v overwrites A(1:k-2,k-1:k).
!>
!>  If UPLO = 'L', then A = L*D*L**T, where
!>     L = P(1)*L(1)* ... *P(k)*L(k)* ...,
!>  i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
!>  n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
!>  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
!>  defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
!>  that if the diagonal block D(k) is of order s (s = 1 or 2), then
!>
!>             (   I    0     0   )  k-1
!>     L(k) =  (   0    I     0   )  s
!>             (   0    v     I   )  n-k-s+1
!>                k-1   s  n-k-s+1
!>
!>  If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
!>  If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
!>  and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
!> 
Contributors:
!>
!>  November 2013,     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 191 of file ssytf2_rook.f.

192*
193* -- LAPACK computational routine --
194* -- LAPACK is a software package provided by Univ. of Tennessee, --
195* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
196*
197* .. Scalar Arguments ..
198 CHARACTER UPLO
199 INTEGER INFO, LDA, N
200* ..
201* .. Array Arguments ..
202 INTEGER IPIV( * )
203 REAL A( LDA, * )
204* ..
205*
206* =====================================================================
207*
208* .. Parameters ..
209 REAL ZERO, ONE
210 parameter( zero = 0.0e+0, one = 1.0e+0 )
211 REAL EIGHT, SEVTEN
212 parameter( eight = 8.0e+0, sevten = 17.0e+0 )
213* ..
214* .. Local Scalars ..
215 LOGICAL UPPER, DONE
216 INTEGER I, IMAX, J, JMAX, ITEMP, K, KK, KP, KSTEP,
217 $ P, II
218 REAL ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22,
219 $ ROWMAX, STEMP, T, WK, WKM1, WKP1, SFMIN
220* ..
221* .. External Functions ..
222 LOGICAL LSAME
223 INTEGER ISAMAX
224 REAL SLAMCH
225 EXTERNAL lsame, isamax, slamch
226* ..
227* .. External Subroutines ..
228 EXTERNAL sscal, sswap, ssyr, xerbla
229* ..
230* .. Intrinsic Functions ..
231 INTRINSIC abs, max, sqrt
232* ..
233* .. Executable Statements ..
234*
235* Test the input parameters.
236*
237 info = 0
238 upper = lsame( uplo, 'U' )
239 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
240 info = -1
241 ELSE IF( n.LT.0 ) THEN
242 info = -2
243 ELSE IF( lda.LT.max( 1, n ) ) THEN
244 info = -4
245 END IF
246 IF( info.NE.0 ) THEN
247 CALL xerbla( 'SSYTF2_ROOK', -info )
248 RETURN
249 END IF
250*
251* Initialize ALPHA for use in choosing pivot block size.
252*
253 alpha = ( one+sqrt( sevten ) ) / eight
254*
255* Compute machine safe minimum
256*
257 sfmin = slamch( 'S' )
258*
259 IF( upper ) THEN
260*
261* Factorize A as U*D*U**T using the upper triangle of A
262*
263* K is the main loop index, decreasing from N to 1 in steps of
264* 1 or 2
265*
266 k = n
267 10 CONTINUE
268*
269* If K < 1, exit from loop
270*
271 IF( k.LT.1 )
272 $ GO TO 70
273 kstep = 1
274 p = k
275*
276* Determine rows and columns to be interchanged and whether
277* a 1-by-1 or 2-by-2 pivot block will be used
278*
279 absakk = abs( a( k, k ) )
280*
281* IMAX is the row-index of the largest off-diagonal element in
282* column K, and COLMAX is its absolute value.
283* Determine both COLMAX and IMAX.
284*
285 IF( k.GT.1 ) THEN
286 imax = isamax( k-1, a( 1, k ), 1 )
287 colmax = abs( a( imax, k ) )
288 ELSE
289 colmax = zero
290 END IF
291*
292 IF( (max( absakk, colmax ).EQ.zero) ) THEN
293*
294* Column K is zero or underflow: set INFO and continue
295*
296 IF( info.EQ.0 )
297 $ info = k
298 kp = k
299 ELSE
300*
301* Test for interchange
302*
303* Equivalent to testing for (used to handle NaN and Inf)
304* ABSAKK.GE.ALPHA*COLMAX
305*
306 IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
307*
308* no interchange,
309* use 1-by-1 pivot block
310*
311 kp = k
312 ELSE
313*
314 done = .false.
315*
316* Loop until pivot found
317*
318 12 CONTINUE
319*
320* Begin pivot search loop body
321*
322* JMAX is the column-index of the largest off-diagonal
323* element in row IMAX, and ROWMAX is its absolute value.
324* Determine both ROWMAX and JMAX.
325*
326 IF( imax.NE.k ) THEN
327 jmax = imax + isamax( k-imax, a( imax, imax+1 ),
328 $ lda )
329 rowmax = abs( a( imax, jmax ) )
330 ELSE
331 rowmax = zero
332 END IF
333*
334 IF( imax.GT.1 ) THEN
335 itemp = isamax( imax-1, a( 1, imax ), 1 )
336 stemp = abs( a( itemp, imax ) )
337 IF( stemp.GT.rowmax ) THEN
338 rowmax = stemp
339 jmax = itemp
340 END IF
341 END IF
342*
343* Equivalent to testing for (used to handle NaN and Inf)
344* ABS( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX
345*
346 IF( .NOT.( abs( a( imax, imax ) ).LT.alpha*rowmax ) )
347 $ THEN
348*
349* interchange rows and columns K and IMAX,
350* use 1-by-1 pivot block
351*
352 kp = imax
353 done = .true.
354*
355* Equivalent to testing for ROWMAX .EQ. COLMAX,
356* used to handle NaN and Inf
357*
358 ELSE IF( ( p.EQ.jmax ).OR.( rowmax.LE.colmax ) ) THEN
359*
360* interchange rows and columns K+1 and IMAX,
361* use 2-by-2 pivot block
362*
363 kp = imax
364 kstep = 2
365 done = .true.
366 ELSE
367*
368* Pivot NOT found, set variables and repeat
369*
370 p = imax
371 colmax = rowmax
372 imax = jmax
373 END IF
374*
375* End pivot search loop body
376*
377 IF( .NOT. done ) GOTO 12
378*
379 END IF
380*
381* Swap TWO rows and TWO columns
382*
383* First swap
384*
385 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
386*
387* Interchange rows and column K and P in the leading
388* submatrix A(1:k,1:k) if we have a 2-by-2 pivot
389*
390 IF( p.GT.1 )
391 $ CALL sswap( p-1, a( 1, k ), 1, a( 1, p ), 1 )
392 IF( p.LT.(k-1) )
393 $ CALL sswap( k-p-1, a( p+1, k ), 1, a( p, p+1 ),
394 $ lda )
395 t = a( k, k )
396 a( k, k ) = a( p, p )
397 a( p, p ) = t
398 END IF
399*
400* Second swap
401*
402 kk = k - kstep + 1
403 IF( kp.NE.kk ) THEN
404*
405* Interchange rows and columns KK and KP in the leading
406* submatrix A(1:k,1:k)
407*
408 IF( kp.GT.1 )
409 $ CALL sswap( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
410 IF( ( kk.GT.1 ) .AND. ( kp.LT.(kk-1) ) )
411 $ CALL sswap( kk-kp-1, a( kp+1, kk ), 1, a( kp,
412 $ kp+1 ),
413 $ lda )
414 t = a( kk, kk )
415 a( kk, kk ) = a( kp, kp )
416 a( kp, kp ) = t
417 IF( kstep.EQ.2 ) THEN
418 t = a( k-1, k )
419 a( k-1, k ) = a( kp, k )
420 a( kp, k ) = t
421 END IF
422 END IF
423*
424* Update the leading submatrix
425*
426 IF( kstep.EQ.1 ) THEN
427*
428* 1-by-1 pivot block D(k): column k now holds
429*
430* W(k) = U(k)*D(k)
431*
432* where U(k) is the k-th column of U
433*
434 IF( k.GT.1 ) THEN
435*
436* Perform a rank-1 update of A(1:k-1,1:k-1) and
437* store U(k) in column k
438*
439 IF( abs( a( k, k ) ).GE.sfmin ) THEN
440*
441* Perform a rank-1 update of A(1:k-1,1:k-1) as
442* A := A - U(k)*D(k)*U(k)**T
443* = A - W(k)*1/D(k)*W(k)**T
444*
445 d11 = one / a( k, k )
446 CALL ssyr( uplo, k-1, -d11, a( 1, k ), 1, a,
447 $ lda )
448*
449* Store U(k) in column k
450*
451 CALL sscal( k-1, d11, a( 1, k ), 1 )
452 ELSE
453*
454* Store L(k) in column K
455*
456 d11 = a( k, k )
457 DO 16 ii = 1, k - 1
458 a( ii, k ) = a( ii, k ) / d11
459 16 CONTINUE
460*
461* Perform a rank-1 update of A(k+1:n,k+1:n) as
462* A := A - U(k)*D(k)*U(k)**T
463* = A - W(k)*(1/D(k))*W(k)**T
464* = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
465*
466 CALL ssyr( uplo, k-1, -d11, a( 1, k ), 1, a,
467 $ lda )
468 END IF
469 END IF
470*
471 ELSE
472*
473* 2-by-2 pivot block D(k): columns k and k-1 now hold
474*
475* ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
476*
477* where U(k) and U(k-1) are the k-th and (k-1)-th columns
478* of U
479*
480* Perform a rank-2 update of A(1:k-2,1:k-2) as
481*
482* A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T
483* = A - ( ( A(k-1)A(k) )*inv(D(k)) ) * ( A(k-1)A(k) )**T
484*
485* and store L(k) and L(k+1) in columns k and k+1
486*
487 IF( k.GT.2 ) THEN
488*
489 d12 = a( k-1, k )
490 d22 = a( k-1, k-1 ) / d12
491 d11 = a( k, k ) / d12
492 t = one / ( d11*d22-one )
493*
494 DO 30 j = k - 2, 1, -1
495*
496 wkm1 = t*( d11*a( j, k-1 )-a( j, k ) )
497 wk = t*( d22*a( j, k )-a( j, k-1 ) )
498*
499 DO 20 i = j, 1, -1
500 a( i, j ) = a( i, j ) - (a( i, k ) / d12 )*wk -
501 $ ( a( i, k-1 ) / d12 )*wkm1
502 20 CONTINUE
503*
504* Store U(k) and U(k-1) in cols k and k-1 for row J
505*
506 a( j, k ) = wk / d12
507 a( j, k-1 ) = wkm1 / d12
508*
509 30 CONTINUE
510*
511 END IF
512*
513 END IF
514 END IF
515*
516* Store details of the interchanges in IPIV
517*
518 IF( kstep.EQ.1 ) THEN
519 ipiv( k ) = kp
520 ELSE
521 ipiv( k ) = -p
522 ipiv( k-1 ) = -kp
523 END IF
524*
525* Decrease K and return to the start of the main loop
526*
527 k = k - kstep
528 GO TO 10
529*
530 ELSE
531*
532* Factorize A as L*D*L**T using the lower triangle of A
533*
534* K is the main loop index, increasing from 1 to N in steps of
535* 1 or 2
536*
537 k = 1
538 40 CONTINUE
539*
540* If K > N, exit from loop
541*
542 IF( k.GT.n )
543 $ GO TO 70
544 kstep = 1
545 p = k
546*
547* Determine rows and columns to be interchanged and whether
548* a 1-by-1 or 2-by-2 pivot block will be used
549*
550 absakk = abs( a( k, k ) )
551*
552* IMAX is the row-index of the largest off-diagonal element in
553* column K, and COLMAX is its absolute value.
554* Determine both COLMAX and IMAX.
555*
556 IF( k.LT.n ) THEN
557 imax = k + isamax( n-k, a( k+1, k ), 1 )
558 colmax = abs( a( imax, k ) )
559 ELSE
560 colmax = zero
561 END IF
562*
563 IF( ( max( absakk, colmax ).EQ.zero ) ) THEN
564*
565* Column K is zero or underflow: set INFO and continue
566*
567 IF( info.EQ.0 )
568 $ info = k
569 kp = k
570 ELSE
571*
572* Test for interchange
573*
574* Equivalent to testing for (used to handle NaN and Inf)
575* ABSAKK.GE.ALPHA*COLMAX
576*
577 IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
578*
579* no interchange, use 1-by-1 pivot block
580*
581 kp = k
582 ELSE
583*
584 done = .false.
585*
586* Loop until pivot found
587*
588 42 CONTINUE
589*
590* Begin pivot search loop body
591*
592* JMAX is the column-index of the largest off-diagonal
593* element in row IMAX, and ROWMAX is its absolute value.
594* Determine both ROWMAX and JMAX.
595*
596 IF( imax.NE.k ) THEN
597 jmax = k - 1 + isamax( imax-k, a( imax, k ),
598 $ lda )
599 rowmax = abs( a( imax, jmax ) )
600 ELSE
601 rowmax = zero
602 END IF
603*
604 IF( imax.LT.n ) THEN
605 itemp = imax + isamax( n-imax, a( imax+1,
606 $ imax ),
607 $ 1 )
608 stemp = abs( a( itemp, imax ) )
609 IF( stemp.GT.rowmax ) THEN
610 rowmax = stemp
611 jmax = itemp
612 END IF
613 END IF
614*
615* Equivalent to testing for (used to handle NaN and Inf)
616* ABS( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX
617*
618 IF( .NOT.( abs( a( imax, imax ) ).LT.alpha*rowmax ) )
619 $ THEN
620*
621* interchange rows and columns K and IMAX,
622* use 1-by-1 pivot block
623*
624 kp = imax
625 done = .true.
626*
627* Equivalent to testing for ROWMAX .EQ. COLMAX,
628* used to handle NaN and Inf
629*
630 ELSE IF( ( p.EQ.jmax ).OR.( rowmax.LE.colmax ) ) THEN
631*
632* interchange rows and columns K+1 and IMAX,
633* use 2-by-2 pivot block
634*
635 kp = imax
636 kstep = 2
637 done = .true.
638 ELSE
639*
640* Pivot NOT found, set variables and repeat
641*
642 p = imax
643 colmax = rowmax
644 imax = jmax
645 END IF
646*
647* End pivot search loop body
648*
649 IF( .NOT. done ) GOTO 42
650*
651 END IF
652*
653* Swap TWO rows and TWO columns
654*
655* First swap
656*
657 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
658*
659* Interchange rows and column K and P in the trailing
660* submatrix A(k:n,k:n) if we have a 2-by-2 pivot
661*
662 IF( p.LT.n )
663 $ CALL sswap( n-p, a( p+1, k ), 1, a( p+1, p ), 1 )
664 IF( p.GT.(k+1) )
665 $ CALL sswap( p-k-1, a( k+1, k ), 1, a( p, k+1 ),
666 $ lda )
667 t = a( k, k )
668 a( k, k ) = a( p, p )
669 a( p, p ) = t
670 END IF
671*
672* Second swap
673*
674 kk = k + kstep - 1
675 IF( kp.NE.kk ) THEN
676*
677* Interchange rows and columns KK and KP in the trailing
678* submatrix A(k:n,k:n)
679*
680 IF( kp.LT.n )
681 $ CALL sswap( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ),
682 $ 1 )
683 IF( ( kk.LT.n ) .AND. ( kp.GT.(kk+1) ) )
684 $ CALL sswap( kp-kk-1, a( kk+1, kk ), 1, a( kp,
685 $ kk+1 ),
686 $ lda )
687 t = a( kk, kk )
688 a( kk, kk ) = a( kp, kp )
689 a( kp, kp ) = t
690 IF( kstep.EQ.2 ) THEN
691 t = a( k+1, k )
692 a( k+1, k ) = a( kp, k )
693 a( kp, k ) = t
694 END IF
695 END IF
696*
697* Update the trailing submatrix
698*
699 IF( kstep.EQ.1 ) THEN
700*
701* 1-by-1 pivot block D(k): column k now holds
702*
703* W(k) = L(k)*D(k)
704*
705* where L(k) is the k-th column of L
706*
707 IF( k.LT.n ) THEN
708*
709* Perform a rank-1 update of A(k+1:n,k+1:n) and
710* store L(k) in column k
711*
712 IF( abs( a( k, k ) ).GE.sfmin ) THEN
713*
714* Perform a rank-1 update of A(k+1:n,k+1:n) as
715* A := A - L(k)*D(k)*L(k)**T
716* = A - W(k)*(1/D(k))*W(k)**T
717*
718 d11 = one / a( k, k )
719 CALL ssyr( uplo, n-k, -d11, a( k+1, k ), 1,
720 $ a( k+1, k+1 ), lda )
721*
722* Store L(k) in column k
723*
724 CALL sscal( n-k, d11, a( k+1, k ), 1 )
725 ELSE
726*
727* Store L(k) in column k
728*
729 d11 = a( k, k )
730 DO 46 ii = k + 1, n
731 a( ii, k ) = a( ii, k ) / d11
732 46 CONTINUE
733*
734* Perform a rank-1 update of A(k+1:n,k+1:n) as
735* A := A - L(k)*D(k)*L(k)**T
736* = A - W(k)*(1/D(k))*W(k)**T
737* = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
738*
739 CALL ssyr( uplo, n-k, -d11, a( k+1, k ), 1,
740 $ a( k+1, k+1 ), lda )
741 END IF
742 END IF
743*
744 ELSE
745*
746* 2-by-2 pivot block D(k): columns k and k+1 now hold
747*
748* ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
749*
750* where L(k) and L(k+1) are the k-th and (k+1)-th columns
751* of L
752*
753*
754* Perform a rank-2 update of A(k+2:n,k+2:n) as
755*
756* A := A - ( L(k) L(k+1) ) * D(k) * ( L(k) L(k+1) )**T
757* = A - ( ( A(k)A(k+1) )*inv(D(k) ) * ( A(k)A(k+1) )**T
758*
759* and store L(k) and L(k+1) in columns k and k+1
760*
761 IF( k.LT.n-1 ) THEN
762*
763 d21 = a( k+1, k )
764 d11 = a( k+1, k+1 ) / d21
765 d22 = a( k, k ) / d21
766 t = one / ( d11*d22-one )
767*
768 DO 60 j = k + 2, n
769*
770* Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J
771*
772 wk = t*( d11*a( j, k )-a( j, k+1 ) )
773 wkp1 = t*( d22*a( j, k+1 )-a( j, k ) )
774*
775* Perform a rank-2 update of A(k+2:n,k+2:n)
776*
777 DO 50 i = j, n
778 a( i, j ) = a( i, j ) - ( a( i, k ) / d21 )*wk -
779 $ ( a( i, k+1 ) / d21 )*wkp1
780 50 CONTINUE
781*
782* Store L(k) and L(k+1) in cols k and k+1 for row J
783*
784 a( j, k ) = wk / d21
785 a( j, k+1 ) = wkp1 / d21
786*
787 60 CONTINUE
788*
789 END IF
790*
791 END IF
792 END IF
793*
794* Store details of the interchanges in IPIV
795*
796 IF( kstep.EQ.1 ) THEN
797 ipiv( k ) = kp
798 ELSE
799 ipiv( k ) = -p
800 ipiv( k+1 ) = -kp
801 END IF
802*
803* Increase K and return to the start of the main loop
804*
805 k = k + kstep
806 GO TO 40
807*
808 END IF
809*
810 70 CONTINUE
811*
812 RETURN
813*
814* End of SSYTF2_ROOK
815*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine ssyr(uplo, n, alpha, x, incx, a, lda)
SSYR
Definition ssyr.f:132
integer function isamax(n, sx, incx)
ISAMAX
Definition isamax.f:71
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
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: