LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine slasyf_rook ( character  UPLO,
integer  N,
integer  NB,
integer  KB,
real, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
real, dimension( ldw, * )  W,
integer  LDW,
integer  INFO 
)

SLASYF_ROOK computes a partial factorization of a real symmetric matrix using the bounded Bunch-Kaufman ("rook") diagonal pivoting method.

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

Purpose: 
 SLASYF_ROOK computes a partial factorization of a real symmetric
 matrix A using the bounded Bunch-Kaufman ("rook") diagonal
 pivoting method. The partial factorization has the form:

 A  =  ( I  U12 ) ( A11  0  ) (  I       0    )  if UPLO = 'U', or:
       ( 0  U22 ) (  0   D  ) ( U12**T U22**T )

 A  =  ( L11  0 ) (  D   0  ) ( L11**T L21**T )  if UPLO = 'L'
       ( L21  I ) (  0  A22 ) (  0       I    )

 where the order of D is at most NB. The actual order is returned in
 the argument KB, and is either NB or NB-1, or N if N <= NB.

 SLASYF_ROOK is an auxiliary routine called by SSYTRF_ROOK. It uses
 blocked code (calling Level 3 BLAS) to update the submatrix
 A11 (if UPLO = 'U') or A22 (if UPLO = 'L').
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]NB
          NB is INTEGER
          The maximum number of columns of the matrix A that should be
          factored.  NB should be at least 2 to allow for 2-by-2 pivot
          blocks.
[out]KB
          KB is INTEGER
          The number of columns of A that were actually factored.
          KB is either NB-1 or NB, or N if N <= NB.
[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, A contains details of the partial factorization.
[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':
             Only the last KB elements of IPIV are set.

             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':
             Only the first KB elements of IPIV are set.

             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]W
          W is REAL array, dimension (LDW,NB)
[in]LDW
          LDW is INTEGER
          The leading dimension of the array W.  LDW >= max(1,N).
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          > 0: if INFO = k, D(k,k) is exactly zero.  The factorization
               has been completed, but the block diagonal matrix D is
               exactly singular.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2013
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

Definition at line 186 of file slasyf_rook.f.

186 *
187 * -- LAPACK computational routine (version 3.5.0) --
188 * -- LAPACK is a software package provided by Univ. of Tennessee, --
189 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
190 * November 2013
191 *
192 * .. Scalar Arguments ..
193  CHARACTER uplo
194  INTEGER info, kb, lda, ldw, n, nb
195 * ..
196 * .. Array Arguments ..
197  INTEGER ipiv( * )
198  REAL a( lda, * ), w( ldw, * )
199 * ..
200 *
201 * =====================================================================
202 *
203 * .. Parameters ..
204  REAL zero, one
205  parameter ( zero = 0.0e+0, one = 1.0e+0 )
206  REAL eight, sevten
207  parameter ( eight = 8.0e+0, sevten = 17.0e+0 )
208 * ..
209 * .. Local Scalars ..
210  LOGICAL done
211  INTEGER imax, itemp, j, jb, jj, jmax, jp1, jp2, k, kk,
212  $ kw, kkw, kp, kstep, p, ii
213 
214  REAL absakk, alpha, colmax, d11, d12, d21, d22,
215  $ stemp, r1, rowmax, t, sfmin
216 * ..
217 * .. External Functions ..
218  LOGICAL lsame
219  INTEGER isamax
220  REAL slamch
221  EXTERNAL lsame, isamax, slamch
222 * ..
223 * .. External Subroutines ..
224  EXTERNAL scopy, sgemm, sgemv, sscal, sswap
225 * ..
226 * .. Intrinsic Functions ..
227  INTRINSIC abs, max, min, sqrt
228 * ..
229 * .. Executable Statements ..
230 *
231  info = 0
232 *
233 * Initialize ALPHA for use in choosing pivot block size.
234 *
235  alpha = ( one+sqrt( sevten ) ) / eight
236 *
237 * Compute machine safe minimum
238 *
239  sfmin = slamch( 'S' )
240 *
241  IF( lsame( uplo, 'U' ) ) THEN
242 *
243 * Factorize the trailing columns of A using the upper triangle
244 * of A and working backwards, and compute the matrix W = U12*D
245 * for use in updating A11
246 *
247 * K is the main loop index, decreasing from N in steps of 1 or 2
248 *
249  k = n
250  10 CONTINUE
251 *
252 * KW is the column of W which corresponds to column K of A
253 *
254  kw = nb + k - n
255 *
256 * Exit from loop
257 *
258  IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
259  $ GO TO 30
260 *
261  kstep = 1
262  p = k
263 *
264 * Copy column K of A to column KW of W and update it
265 *
266  CALL scopy( k, a( 1, k ), 1, w( 1, kw ), 1 )
267  IF( k.LT.n )
268  $ CALL sgemv( 'No transpose', k, n-k, -one, a( 1, k+1 ),
269  $ lda, w( k, kw+1 ), ldw, one, w( 1, kw ), 1 )
270 *
271 * Determine rows and columns to be interchanged and whether
272 * a 1-by-1 or 2-by-2 pivot block will be used
273 *
274  absakk = abs( w( k, kw ) )
275 *
276 * IMAX is the row-index of the largest off-diagonal element in
277 * column K, and COLMAX is its absolute value.
278 * Determine both COLMAX and IMAX.
279 *
280  IF( k.GT.1 ) THEN
281  imax = isamax( k-1, w( 1, kw ), 1 )
282  colmax = abs( w( imax, kw ) )
283  ELSE
284  colmax = zero
285  END IF
286 *
287  IF( max( absakk, colmax ).EQ.zero ) THEN
288 *
289 * Column K is zero or underflow: set INFO and continue
290 *
291  IF( info.EQ.0 )
292  $ info = k
293  kp = k
294  CALL scopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
295  ELSE
296 *
297 * ============================================================
298 *
299 * Test for interchange
300 *
301 * Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
302 * (used to handle NaN and Inf)
303 *
304  IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
305 *
306 * no interchange, use 1-by-1 pivot block
307 *
308  kp = k
309 *
310  ELSE
311 *
312  done = .false.
313 *
314 * Loop until pivot found
315 *
316  12 CONTINUE
317 *
318 * Begin pivot search loop body
319 *
320 *
321 * Copy column IMAX to column KW-1 of W and update it
322 *
323  CALL scopy( imax, a( 1, imax ), 1, w( 1, kw-1 ), 1 )
324  CALL scopy( k-imax, a( imax, imax+1 ), lda,
325  $ w( imax+1, kw-1 ), 1 )
326 *
327  IF( k.LT.n )
328  $ CALL sgemv( 'No transpose', k, n-k, -one,
329  $ a( 1, k+1 ), lda, w( imax, kw+1 ), ldw,
330  $ one, w( 1, kw-1 ), 1 )
331 *
332 * JMAX is the column-index of the largest off-diagonal
333 * element in row IMAX, and ROWMAX is its absolute value.
334 * Determine both ROWMAX and JMAX.
335 *
336  IF( imax.NE.k ) THEN
337  jmax = imax + isamax( k-imax, w( imax+1, kw-1 ),
338  $ 1 )
339  rowmax = abs( w( jmax, kw-1 ) )
340  ELSE
341  rowmax = zero
342  END IF
343 *
344  IF( imax.GT.1 ) THEN
345  itemp = isamax( imax-1, w( 1, kw-1 ), 1 )
346  stemp = abs( w( itemp, kw-1 ) )
347  IF( stemp.GT.rowmax ) THEN
348  rowmax = stemp
349  jmax = itemp
350  END IF
351  END IF
352 *
353 * Equivalent to testing for
354 * ABS( W( IMAX, KW-1 ) ).GE.ALPHA*ROWMAX
355 * (used to handle NaN and Inf)
356 *
357  IF( .NOT.(abs( w( imax, kw-1 ) ).LT.alpha*rowmax ) )
358  $ THEN
359 *
360 * interchange rows and columns K and IMAX,
361 * use 1-by-1 pivot block
362 *
363  kp = imax
364 *
365 * copy column KW-1 of W to column KW of W
366 *
367  CALL scopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
368 *
369  done = .true.
370 *
371 * Equivalent to testing for ROWMAX.EQ.COLMAX,
372 * (used to handle NaN and Inf)
373 *
374  ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
375  $ THEN
376 *
377 * interchange rows and columns K-1 and IMAX,
378 * use 2-by-2 pivot block
379 *
380  kp = imax
381  kstep = 2
382  done = .true.
383  ELSE
384 *
385 * Pivot not found: set params and repeat
386 *
387  p = imax
388  colmax = rowmax
389  imax = jmax
390 *
391 * Copy updated JMAXth (next IMAXth) column to Kth of W
392 *
393  CALL scopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
394 *
395  END IF
396 *
397 * End pivot search loop body
398 *
399  IF( .NOT. done ) GOTO 12
400 *
401  END IF
402 *
403 * ============================================================
404 *
405  kk = k - kstep + 1
406 *
407 * KKW is the column of W which corresponds to column KK of A
408 *
409  kkw = nb + kk - n
410 *
411  IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
412 *
413 * Copy non-updated column K to column P
414 *
415  CALL scopy( k-p, a( p+1, k ), 1, a( p, p+1 ), lda )
416  CALL scopy( p, a( 1, k ), 1, a( 1, p ), 1 )
417 *
418 * Interchange rows K and P in last N-K+1 columns of A
419 * and last N-K+2 columns of W
420 *
421  CALL sswap( n-k+1, a( k, k ), lda, a( p, k ), lda )
422  CALL sswap( n-kk+1, w( k, kkw ), ldw, w( p, kkw ), ldw )
423  END IF
424 *
425 * Updated column KP is already stored in column KKW of W
426 *
427  IF( kp.NE.kk ) THEN
428 *
429 * Copy non-updated column KK to column KP
430 *
431  a( kp, k ) = a( kk, k )
432  CALL scopy( k-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
433  $ lda )
434  CALL scopy( kp, a( 1, kk ), 1, a( 1, kp ), 1 )
435 *
436 * Interchange rows KK and KP in last N-KK+1 columns
437 * of A and W
438 *
439  CALL sswap( n-kk+1, a( kk, kk ), lda, a( kp, kk ), lda )
440  CALL sswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
441  $ ldw )
442  END IF
443 *
444  IF( kstep.EQ.1 ) THEN
445 *
446 * 1-by-1 pivot block D(k): column KW of W now holds
447 *
448 * W(k) = U(k)*D(k)
449 *
450 * where U(k) is the k-th column of U
451 *
452 * Store U(k) in column k of A
453 *
454  CALL scopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
455  IF( k.GT.1 ) THEN
456  IF( abs( a( k, k ) ).GE.sfmin ) THEN
457  r1 = one / a( k, k )
458  CALL sscal( k-1, r1, a( 1, k ), 1 )
459  ELSE IF( a( k, k ).NE.zero ) THEN
460  DO 14 ii = 1, k - 1
461  a( ii, k ) = a( ii, k ) / a( k, k )
462  14 CONTINUE
463  END IF
464  END IF
465 *
466  ELSE
467 *
468 * 2-by-2 pivot block D(k): columns KW and KW-1 of W now
469 * hold
470 *
471 * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
472 *
473 * where U(k) and U(k-1) are the k-th and (k-1)-th columns
474 * of U
475 *
476  IF( k.GT.2 ) THEN
477 *
478 * Store U(k) and U(k-1) in columns k and k-1 of A
479 *
480  d12 = w( k-1, kw )
481  d11 = w( k, kw ) / d12
482  d22 = w( k-1, kw-1 ) / d12
483  t = one / ( d11*d22-one )
484  DO 20 j = 1, k - 2
485  a( j, k-1 ) = t*( (d11*w( j, kw-1 )-w( j, kw ) ) /
486  $ d12 )
487  a( j, k ) = t*( ( d22*w( j, kw )-w( j, kw-1 ) ) /
488  $ d12 )
489  20 CONTINUE
490  END IF
491 *
492 * Copy D(k) to A
493 *
494  a( k-1, k-1 ) = w( k-1, kw-1 )
495  a( k-1, k ) = w( k-1, kw )
496  a( k, k ) = w( k, kw )
497  END IF
498  END IF
499 *
500 * Store details of the interchanges in IPIV
501 *
502  IF( kstep.EQ.1 ) THEN
503  ipiv( k ) = kp
504  ELSE
505  ipiv( k ) = -p
506  ipiv( k-1 ) = -kp
507  END IF
508 *
509 * Decrease K and return to the start of the main loop
510 *
511  k = k - kstep
512  GO TO 10
513 *
514  30 CONTINUE
515 *
516 * Update the upper triangle of A11 (= A(1:k,1:k)) as
517 *
518 * A11 := A11 - U12*D*U12**T = A11 - U12*W**T
519 *
520 * computing blocks of NB columns at a time
521 *
522  DO 50 j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
523  jb = min( nb, k-j+1 )
524 *
525 * Update the upper triangle of the diagonal block
526 *
527  DO 40 jj = j, j + jb - 1
528  CALL sgemv( 'No transpose', jj-j+1, n-k, -one,
529  $ a( j, k+1 ), lda, w( jj, kw+1 ), ldw, one,
530  $ a( j, jj ), 1 )
531  40 CONTINUE
532 *
533 * Update the rectangular superdiagonal block
534 *
535  IF( j.GE.2 )
536  $ CALL sgemm( 'No transpose', 'Transpose', j-1, jb,
537  $ n-k, -one, a( 1, k+1 ), lda, w( j, kw+1 ), ldw,
538  $ one, a( 1, j ), lda )
539  50 CONTINUE
540 *
541 * Put U12 in standard form by partially undoing the interchanges
542 * in columns k+1:n
543 *
544  j = k + 1
545  60 CONTINUE
546 *
547  kstep = 1
548  jp1 = 1
549  jj = j
550  jp2 = ipiv( j )
551  IF( jp2.LT.0 ) THEN
552  jp2 = -jp2
553  j = j + 1
554  jp1 = -ipiv( j )
555  kstep = 2
556  END IF
557 *
558  j = j + 1
559  IF( jp2.NE.jj .AND. j.LE.n )
560  $ CALL sswap( n-j+1, a( jp2, j ), lda, a( jj, j ), lda )
561  jj = j - 1
562  IF( jp1.NE.jj .AND. kstep.EQ.2 )
563  $ CALL sswap( n-j+1, a( jp1, j ), lda, a( jj, j ), lda )
564  IF( j.LE.n )
565  $ GO TO 60
566 *
567 * Set KB to the number of columns factorized
568 *
569  kb = n - k
570 *
571  ELSE
572 *
573 * Factorize the leading columns of A using the lower triangle
574 * of A and working forwards, and compute the matrix W = L21*D
575 * for use in updating A22
576 *
577 * K is the main loop index, increasing from 1 in steps of 1 or 2
578 *
579  k = 1
580  70 CONTINUE
581 *
582 * Exit from loop
583 *
584  IF( ( k.GE.nb .AND. nb.LT.n ) .OR. k.GT.n )
585  $ GO TO 90
586 *
587  kstep = 1
588  p = k
589 *
590 * Copy column K of A to column K of W and update it
591 *
592  CALL scopy( n-k+1, a( k, k ), 1, w( k, k ), 1 )
593  IF( k.GT.1 )
594  $ CALL sgemv( 'No transpose', n-k+1, k-1, -one, a( k, 1 ),
595  $ lda, w( k, 1 ), ldw, one, w( k, k ), 1 )
596 *
597 * Determine rows and columns to be interchanged and whether
598 * a 1-by-1 or 2-by-2 pivot block will be used
599 *
600  absakk = abs( w( k, k ) )
601 *
602 * IMAX is the row-index of the largest off-diagonal element in
603 * column K, and COLMAX is its absolute value.
604 * Determine both COLMAX and IMAX.
605 *
606  IF( k.LT.n ) THEN
607  imax = k + isamax( n-k, w( k+1, k ), 1 )
608  colmax = abs( w( imax, k ) )
609  ELSE
610  colmax = zero
611  END IF
612 *
613  IF( max( absakk, colmax ).EQ.zero ) THEN
614 *
615 * Column K is zero or underflow: set INFO and continue
616 *
617  IF( info.EQ.0 )
618  $ info = k
619  kp = k
620  CALL scopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
621  ELSE
622 *
623 * ============================================================
624 *
625 * Test for interchange
626 *
627 * Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
628 * (used to handle NaN and Inf)
629 *
630  IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
631 *
632 * no interchange, use 1-by-1 pivot block
633 *
634  kp = k
635 *
636  ELSE
637 *
638  done = .false.
639 *
640 * Loop until pivot found
641 *
642  72 CONTINUE
643 *
644 * Begin pivot search loop body
645 *
646 *
647 * Copy column IMAX to column K+1 of W and update it
648 *
649  CALL scopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1)
650  CALL scopy( n-imax+1, a( imax, imax ), 1,
651  $ w( imax, k+1 ), 1 )
652  IF( k.GT.1 )
653  $ CALL sgemv( 'No transpose', n-k+1, k-1, -one,
654  $ a( k, 1 ), lda, w( imax, 1 ), ldw,
655  $ one, w( k, k+1 ), 1 )
656 *
657 * JMAX is the column-index of the largest off-diagonal
658 * element in row IMAX, and ROWMAX is its absolute value.
659 * Determine both ROWMAX and JMAX.
660 *
661  IF( imax.NE.k ) THEN
662  jmax = k - 1 + isamax( imax-k, w( k, k+1 ), 1 )
663  rowmax = abs( w( jmax, k+1 ) )
664  ELSE
665  rowmax = zero
666  END IF
667 *
668  IF( imax.LT.n ) THEN
669  itemp = imax + isamax( n-imax, w( imax+1, k+1 ), 1)
670  stemp = abs( w( itemp, k+1 ) )
671  IF( stemp.GT.rowmax ) THEN
672  rowmax = stemp
673  jmax = itemp
674  END IF
675  END IF
676 *
677 * Equivalent to testing for
678 * ABS( W( IMAX, K+1 ) ).GE.ALPHA*ROWMAX
679 * (used to handle NaN and Inf)
680 *
681  IF( .NOT.( abs( w( imax, k+1 ) ).LT.alpha*rowmax ) )
682  $ THEN
683 *
684 * interchange rows and columns K and IMAX,
685 * use 1-by-1 pivot block
686 *
687  kp = imax
688 *
689 * copy column K+1 of W to column K of W
690 *
691  CALL scopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
692 *
693  done = .true.
694 *
695 * Equivalent to testing for ROWMAX.EQ.COLMAX,
696 * (used to handle NaN and Inf)
697 *
698  ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
699  $ THEN
700 *
701 * interchange rows and columns K+1 and IMAX,
702 * use 2-by-2 pivot block
703 *
704  kp = imax
705  kstep = 2
706  done = .true.
707  ELSE
708 *
709 * Pivot not found: set params and repeat
710 *
711  p = imax
712  colmax = rowmax
713  imax = jmax
714 *
715 * Copy updated JMAXth (next IMAXth) column to Kth of W
716 *
717  CALL scopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
718 *
719  END IF
720 *
721 * End pivot search loop body
722 *
723  IF( .NOT. done ) GOTO 72
724 *
725  END IF
726 *
727 * ============================================================
728 *
729  kk = k + kstep - 1
730 *
731  IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
732 *
733 * Copy non-updated column K to column P
734 *
735  CALL scopy( p-k, a( k, k ), 1, a( p, k ), lda )
736  CALL scopy( n-p+1, a( p, k ), 1, a( p, p ), 1 )
737 *
738 * Interchange rows K and P in first K columns of A
739 * and first K+1 columns of W
740 *
741  CALL sswap( k, a( k, 1 ), lda, a( p, 1 ), lda )
742  CALL sswap( kk, w( k, 1 ), ldw, w( p, 1 ), ldw )
743  END IF
744 *
745 * Updated column KP is already stored in column KK of W
746 *
747  IF( kp.NE.kk ) THEN
748 *
749 * Copy non-updated column KK to column KP
750 *
751  a( kp, k ) = a( kk, k )
752  CALL scopy( kp-k-1, a( k+1, kk ), 1, a( kp, k+1 ), lda )
753  CALL scopy( n-kp+1, a( kp, kk ), 1, a( kp, kp ), 1 )
754 *
755 * Interchange rows KK and KP in first KK columns of A and W
756 *
757  CALL sswap( kk, a( kk, 1 ), lda, a( kp, 1 ), lda )
758  CALL sswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
759  END IF
760 *
761  IF( kstep.EQ.1 ) THEN
762 *
763 * 1-by-1 pivot block D(k): column k of W now holds
764 *
765 * W(k) = L(k)*D(k)
766 *
767 * where L(k) is the k-th column of L
768 *
769 * Store L(k) in column k of A
770 *
771  CALL scopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
772  IF( k.LT.n ) THEN
773  IF( abs( a( k, k ) ).GE.sfmin ) THEN
774  r1 = one / a( k, k )
775  CALL sscal( n-k, r1, a( k+1, k ), 1 )
776  ELSE IF( a( k, k ).NE.zero ) THEN
777  DO 74 ii = k + 1, n
778  a( ii, k ) = a( ii, k ) / a( k, k )
779  74 CONTINUE
780  END IF
781  END IF
782 *
783  ELSE
784 *
785 * 2-by-2 pivot block D(k): columns k and k+1 of W now hold
786 *
787 * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
788 *
789 * where L(k) and L(k+1) are the k-th and (k+1)-th columns
790 * of L
791 *
792  IF( k.LT.n-1 ) THEN
793 *
794 * Store L(k) and L(k+1) in columns k and k+1 of A
795 *
796  d21 = w( k+1, k )
797  d11 = w( k+1, k+1 ) / d21
798  d22 = w( k, k ) / d21
799  t = one / ( d11*d22-one )
800  DO 80 j = k + 2, n
801  a( j, k ) = t*( ( d11*w( j, k )-w( j, k+1 ) ) /
802  $ d21 )
803  a( j, k+1 ) = t*( ( d22*w( j, k+1 )-w( j, k ) ) /
804  $ d21 )
805  80 CONTINUE
806  END IF
807 *
808 * Copy D(k) to A
809 *
810  a( k, k ) = w( k, k )
811  a( k+1, k ) = w( k+1, k )
812  a( k+1, k+1 ) = w( k+1, k+1 )
813  END IF
814  END IF
815 *
816 * Store details of the interchanges in IPIV
817 *
818  IF( kstep.EQ.1 ) THEN
819  ipiv( k ) = kp
820  ELSE
821  ipiv( k ) = -p
822  ipiv( k+1 ) = -kp
823  END IF
824 *
825 * Increase K and return to the start of the main loop
826 *
827  k = k + kstep
828  GO TO 70
829 *
830  90 CONTINUE
831 *
832 * Update the lower triangle of A22 (= A(k:n,k:n)) as
833 *
834 * A22 := A22 - L21*D*L21**T = A22 - L21*W**T
835 *
836 * computing blocks of NB columns at a time
837 *
838  DO 110 j = k, n, nb
839  jb = min( nb, n-j+1 )
840 *
841 * Update the lower triangle of the diagonal block
842 *
843  DO 100 jj = j, j + jb - 1
844  CALL sgemv( 'No transpose', j+jb-jj, k-1, -one,
845  $ a( jj, 1 ), lda, w( jj, 1 ), ldw, one,
846  $ a( jj, jj ), 1 )
847  100 CONTINUE
848 *
849 * Update the rectangular subdiagonal block
850 *
851  IF( j+jb.LE.n )
852  $ CALL sgemm( 'No transpose', 'Transpose', n-j-jb+1, jb,
853  $ k-1, -one, a( j+jb, 1 ), lda, w( j, 1 ), ldw,
854  $ one, a( j+jb, j ), lda )
855  110 CONTINUE
856 *
857 * Put L21 in standard form by partially undoing the interchanges
858 * in columns 1:k-1
859 *
860  j = k - 1
861  120 CONTINUE
862 *
863  kstep = 1
864  jp1 = 1
865  jj = j
866  jp2 = ipiv( j )
867  IF( jp2.LT.0 ) THEN
868  jp2 = -jp2
869  j = j - 1
870  jp1 = -ipiv( j )
871  kstep = 2
872  END IF
873 *
874  j = j - 1
875  IF( jp2.NE.jj .AND. j.GE.1 )
876  $ CALL sswap( j, a( jp2, 1 ), lda, a( jj, 1 ), lda )
877  jj = j + 1
878  IF( jp1.NE.jj .AND. kstep.EQ.2 )
879  $ CALL sswap( j, a( jp1, 1 ), lda, a( jj, 1 ), lda )
880  IF( j.GE.1 )
881  $ GO TO 120
882 *
883 * Set KB to the number of columns factorized
884 *
885  kb = k - 1
886 *
887  END IF
888  RETURN
889 *
890 * End of SLASYF_ROOK
891 *
isamax
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:53
sgemm
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:189
sgemv
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:158
sscal
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:55
sswap
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:53
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
scopy
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53

Here is the call graph for this function:

Here is the caller graph for this function: