LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ dlasyf_rook()

 subroutine dlasyf_rook ( character uplo, integer n, integer nb, integer kb, double precision, dimension( lda, * ) a, integer lda, integer, dimension( * ) ipiv, double precision, dimension( ldw, * ) w, integer ldw, integer info )

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

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

Purpose:
DLASYF_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.

DLASYF_ROOK is an auxiliary routine called by DSYTRF_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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 dlasyf_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 DOUBLE PRECISION A( LDA, * ), W( LDW, * )
196* ..
197*
198* =====================================================================
199*
200* .. Parameters ..
201 DOUBLE PRECISION ZERO, ONE
202 parameter( zero = 0.0d+0, one = 1.0d+0 )
203 DOUBLE PRECISION EIGHT, SEVTEN
204 parameter( eight = 8.0d+0, sevten = 17.0d+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 DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22,
212 \$ DTEMP, R1, ROWMAX, T, SFMIN
213* ..
214* .. External Functions ..
215 LOGICAL LSAME
216 INTEGER IDAMAX
217 DOUBLE PRECISION DLAMCH
218 EXTERNAL lsame, idamax, dlamch
219* ..
220* .. External Subroutines ..
221 EXTERNAL dcopy, dgemm, dgemv, dscal, dswap
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 = dlamch( '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 dcopy( k, a( 1, k ), 1, w( 1, kw ), 1 )
264 IF( k.LT.n )
265 \$ CALL dgemv( '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 = idamax( 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 dcopy( 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 dcopy( imax, a( 1, imax ), 1, w( 1, kw-1 ), 1 )
321 CALL dcopy( k-imax, a( imax, imax+1 ), lda,
322 \$ w( imax+1, kw-1 ), 1 )
323*
324 IF( k.LT.n )
325 \$ CALL dgemv( '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 + idamax( 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 = idamax( imax-1, w( 1, kw-1 ), 1 )
343 dtemp = abs( w( itemp, kw-1 ) )
344 IF( dtemp.GT.rowmax ) THEN
345 rowmax = dtemp
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 dcopy( 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 dcopy( 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 dcopy( k-p, a( p+1, k ), 1, a( p, p+1 ), lda )
413 CALL dcopy( 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 dswap( n-k+1, a( k, k ), lda, a( p, k ), lda )
419 CALL dswap( 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 dcopy( k-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
430 \$ lda )
431 CALL dcopy( 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 dswap( n-kk+1, a( kk, kk ), lda, a( kp, kk ), lda )
437 CALL dswap( 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 dcopy( 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 dscal( 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 dgemv( '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 dgemm( '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 dswap( 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 dswap( 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 dcopy( n-k+1, a( k, k ), 1, w( k, k ), 1 )
590 IF( k.GT.1 )
591 \$ CALL dgemv( '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 + idamax( 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 dcopy( 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 dcopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1)
647 CALL dcopy( n-imax+1, a( imax, imax ), 1,
648 \$ w( imax, k+1 ), 1 )
649 IF( k.GT.1 )
650 \$ CALL dgemv( '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 + idamax( 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 + idamax( n-imax, w( imax+1, k+1 ), 1)
667 dtemp = abs( w( itemp, k+1 ) )
668 IF( dtemp.GT.rowmax ) THEN
669 rowmax = dtemp
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 dcopy( 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 dcopy( 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 dcopy( p-k, a( k, k ), 1, a( p, k ), lda )
733 CALL dcopy( 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 dswap( k, a( k, 1 ), lda, a( p, 1 ), lda )
739 CALL dswap( 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 dcopy( kp-k-1, a( k+1, kk ), 1, a( kp, k+1 ), lda )
750 CALL dcopy( 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 dswap( kk, a( kk, 1 ), lda, a( kp, 1 ), lda )
755 CALL dswap( 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 dcopy( 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 dscal( 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 dgemv( '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 dgemm( '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 dswap( j, a( jp2, 1 ), lda, a( jj, 1 ), lda )
874 jj = j + 1
875 IF( jp1.NE.jj .AND. kstep.EQ.2 )
876 \$ CALL dswap( 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 DLASYF_ROOK
888*
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:158
integer function idamax(n, dx, incx)
IDAMAX
Definition idamax.f:71
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82
Here is the call graph for this function:
Here is the caller graph for this function: