LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ 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 ("rook") 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 193 of file ssytf2_rook.f.

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