LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
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.

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.```
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*
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*
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*
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
SGEMV
Definition sgemv.f:158
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: