LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ slasyf_rook()

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.
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 182 of file slasyf_rook.f.

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