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

◆ 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: