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

◆ zsytf2_rook()

subroutine zsytf2_rook ( character uplo,
integer n,
complex*16, dimension( lda, * ) a,
integer lda,
integer, dimension( * ) ipiv,
integer info )

ZSYTF2_ROOK computes the factorization of a complex symmetric indefinite matrix using the bounded Bunch-Kaufman ("rook") diagonal pivoting method (unblocked algorithm).

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

Purpose:
!>
!> ZSYTF2_ROOK computes the factorization of a complex symmetric matrix A
!> using the bounded Bunch-Kaufman () diagonal pivoting method:
!>
!>    A = U*D*U**T  or  A = L*D*L**T
!>
!> where U (or L) is a product of permutation and unit upper (lower)
!> triangular matrices, U**T is the transpose of U, and D is symmetric and
!> block diagonal with 1-by-1 and 2-by-2 diagonal blocks.
!>
!> This is the unblocked version of the algorithm, calling Level 2 BLAS.
!> 
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,out]A
!>          A is COMPLEX*16 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, the block diagonal matrix D and the multipliers used
!>          to obtain the factor U or L (see below for further details).
!> 
[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':
!>             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':
!>             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]INFO
!>          INFO is INTEGER
!>          = 0: successful exit
!>          < 0: if INFO = -k, the k-th argument had an illegal value
!>          > 0: if INFO = k, D(k,k) is exactly zero.  The factorization
!>               has been completed, but the block diagonal matrix D is
!>               exactly singular, and division by zero will occur if it
!>               is used to solve a system of equations.
!> 
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  If UPLO = 'U', then A = U*D*U**T, where
!>     U = P(n)*U(n)* ... *P(k)U(k)* ...,
!>  i.e., U is a product of terms P(k)*U(k), where k decreases from n to
!>  1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
!>  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
!>  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
!>  that if the diagonal block D(k) is of order s (s = 1 or 2), then
!>
!>             (   I    v    0   )   k-s
!>     U(k) =  (   0    I    0   )   s
!>             (   0    0    I   )   n-k
!>                k-s   s   n-k
!>
!>  If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
!>  If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
!>  and A(k,k), and v overwrites A(1:k-2,k-1:k).
!>
!>  If UPLO = 'L', then A = L*D*L**T, where
!>     L = P(1)*L(1)* ... *P(k)*L(k)* ...,
!>  i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
!>  n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
!>  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
!>  defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
!>  that if the diagonal block D(k) is of order s (s = 1 or 2), then
!>
!>             (   I    0     0   )  k-1
!>     L(k) =  (   0    I     0   )  s
!>             (   0    v     I   )  n-k-s+1
!>                k-1   s  n-k-s+1
!>
!>  If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
!>  If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
!>  and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
!> 
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
!>
!>  01-01-96 - Based on modifications by
!>    J. Lewis, Boeing Computer Services Company
!>    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville abd , USA
!> 

Definition at line 191 of file zsytf2_rook.f.

192*
193* -- LAPACK computational routine --
194* -- LAPACK is a software package provided by Univ. of Tennessee, --
195* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
196*
197* .. Scalar Arguments ..
198 CHARACTER UPLO
199 INTEGER INFO, LDA, N
200* ..
201* .. Array Arguments ..
202 INTEGER IPIV( * )
203 COMPLEX*16 A( LDA, * )
204* ..
205*
206* =====================================================================
207*
208* .. Parameters ..
209 DOUBLE PRECISION ZERO, ONE
210 parameter( zero = 0.0d+0, one = 1.0d+0 )
211 DOUBLE PRECISION EIGHT, SEVTEN
212 parameter( eight = 8.0d+0, sevten = 17.0d+0 )
213 COMPLEX*16 CONE
214 parameter( cone = ( 1.0d+0, 0.0d+0 ) )
215* ..
216* .. Local Scalars ..
217 LOGICAL UPPER, DONE
218 INTEGER I, IMAX, J, JMAX, ITEMP, K, KK, KP, KSTEP,
219 $ P, II
220 DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, ROWMAX, DTEMP, SFMIN
221 COMPLEX*16 D11, D12, D21, D22, T, WK, WKM1, WKP1, Z
222* ..
223* .. External Functions ..
224 LOGICAL LSAME
225 INTEGER IZAMAX
226 DOUBLE PRECISION DLAMCH
227 EXTERNAL lsame, izamax, dlamch
228* ..
229* .. External Subroutines ..
230 EXTERNAL zscal, zswap, zsyr, xerbla
231* ..
232* .. Intrinsic Functions ..
233 INTRINSIC abs, max, sqrt, dimag, dble
234* ..
235* .. Statement Functions ..
236 DOUBLE PRECISION CABS1
237* ..
238* .. Statement Function definitions ..
239 cabs1( z ) = abs( dble( z ) ) + abs( dimag( z ) )
240* ..
241* .. Executable Statements ..
242*
243* Test the input parameters.
244*
245 info = 0
246 upper = lsame( uplo, 'U' )
247 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
248 info = -1
249 ELSE IF( n.LT.0 ) THEN
250 info = -2
251 ELSE IF( lda.LT.max( 1, n ) ) THEN
252 info = -4
253 END IF
254 IF( info.NE.0 ) THEN
255 CALL xerbla( 'ZSYTF2_ROOK', -info )
256 RETURN
257 END IF
258*
259* Initialize ALPHA for use in choosing pivot block size.
260*
261 alpha = ( one+sqrt( sevten ) ) / eight
262*
263* Compute machine safe minimum
264*
265 sfmin = dlamch( 'S' )
266*
267 IF( upper ) THEN
268*
269* Factorize A as U*D*U**T using the upper triangle of A
270*
271* K is the main loop index, decreasing from N to 1 in steps of
272* 1 or 2
273*
274 k = n
275 10 CONTINUE
276*
277* If K < 1, exit from loop
278*
279 IF( k.LT.1 )
280 $ GO TO 70
281 kstep = 1
282 p = k
283*
284* Determine rows and columns to be interchanged and whether
285* a 1-by-1 or 2-by-2 pivot block will be used
286*
287 absakk = cabs1( a( k, k ) )
288*
289* IMAX is the row-index of the largest off-diagonal element in
290* column K, and COLMAX is its absolute value.
291* Determine both COLMAX and IMAX.
292*
293 IF( k.GT.1 ) THEN
294 imax = izamax( k-1, a( 1, k ), 1 )
295 colmax = cabs1( a( imax, k ) )
296 ELSE
297 colmax = zero
298 END IF
299*
300 IF( (max( absakk, colmax ).EQ.zero) ) THEN
301*
302* Column K is zero or underflow: set INFO and continue
303*
304 IF( info.EQ.0 )
305 $ info = k
306 kp = k
307 ELSE
308*
309* Test for interchange
310*
311* Equivalent to testing for (used to handle NaN and Inf)
312* ABSAKK.GE.ALPHA*COLMAX
313*
314 IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
315*
316* no interchange,
317* use 1-by-1 pivot block
318*
319 kp = k
320 ELSE
321*
322 done = .false.
323*
324* Loop until pivot found
325*
326 12 CONTINUE
327*
328* Begin pivot search loop body
329*
330* JMAX is the column-index of the largest off-diagonal
331* element in row IMAX, and ROWMAX is its absolute value.
332* Determine both ROWMAX and JMAX.
333*
334 IF( imax.NE.k ) THEN
335 jmax = imax + izamax( k-imax, a( imax, imax+1 ),
336 $ lda )
337 rowmax = cabs1( a( imax, jmax ) )
338 ELSE
339 rowmax = zero
340 END IF
341*
342 IF( imax.GT.1 ) THEN
343 itemp = izamax( imax-1, a( 1, imax ), 1 )
344 dtemp = cabs1( a( itemp, imax ) )
345 IF( dtemp.GT.rowmax ) THEN
346 rowmax = dtemp
347 jmax = itemp
348 END IF
349 END IF
350*
351* Equivalent to testing for (used to handle NaN and Inf)
352* CABS1( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX
353*
354 IF( .NOT.( cabs1(a( imax, imax )).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 done = .true.
362*
363* Equivalent to testing for ROWMAX .EQ. COLMAX,
364* used to handle NaN and Inf
365*
366 ELSE IF( ( p.EQ.jmax ).OR.( rowmax.LE.colmax ) ) THEN
367*
368* interchange rows and columns K+1 and IMAX,
369* use 2-by-2 pivot block
370*
371 kp = imax
372 kstep = 2
373 done = .true.
374 ELSE
375*
376* Pivot NOT found, set variables and repeat
377*
378 p = imax
379 colmax = rowmax
380 imax = jmax
381 END IF
382*
383* End pivot search loop body
384*
385 IF( .NOT. done ) GOTO 12
386*
387 END IF
388*
389* Swap TWO rows and TWO columns
390*
391* First swap
392*
393 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
394*
395* Interchange rows and column K and P in the leading
396* submatrix A(1:k,1:k) if we have a 2-by-2 pivot
397*
398 IF( p.GT.1 )
399 $ CALL zswap( p-1, a( 1, k ), 1, a( 1, p ), 1 )
400 IF( p.LT.(k-1) )
401 $ CALL zswap( k-p-1, a( p+1, k ), 1, a( p, p+1 ),
402 $ lda )
403 t = a( k, k )
404 a( k, k ) = a( p, p )
405 a( p, p ) = t
406 END IF
407*
408* Second swap
409*
410 kk = k - kstep + 1
411 IF( kp.NE.kk ) THEN
412*
413* Interchange rows and columns KK and KP in the leading
414* submatrix A(1:k,1:k)
415*
416 IF( kp.GT.1 )
417 $ CALL zswap( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
418 IF( ( kk.GT.1 ) .AND. ( kp.LT.(kk-1) ) )
419 $ CALL zswap( kk-kp-1, a( kp+1, kk ), 1, a( kp,
420 $ kp+1 ),
421 $ lda )
422 t = a( kk, kk )
423 a( kk, kk ) = a( kp, kp )
424 a( kp, kp ) = t
425 IF( kstep.EQ.2 ) THEN
426 t = a( k-1, k )
427 a( k-1, k ) = a( kp, k )
428 a( kp, k ) = t
429 END IF
430 END IF
431*
432* Update the leading submatrix
433*
434 IF( kstep.EQ.1 ) THEN
435*
436* 1-by-1 pivot block D(k): column k now holds
437*
438* W(k) = U(k)*D(k)
439*
440* where U(k) is the k-th column of U
441*
442 IF( k.GT.1 ) THEN
443*
444* Perform a rank-1 update of A(1:k-1,1:k-1) and
445* store U(k) in column k
446*
447 IF( cabs1( a( k, k ) ).GE.sfmin ) THEN
448*
449* Perform a rank-1 update of A(1:k-1,1:k-1) as
450* A := A - U(k)*D(k)*U(k)**T
451* = A - W(k)*1/D(k)*W(k)**T
452*
453 d11 = cone / a( k, k )
454 CALL zsyr( uplo, k-1, -d11, a( 1, k ), 1, a,
455 $ lda )
456*
457* Store U(k) in column k
458*
459 CALL zscal( k-1, d11, a( 1, k ), 1 )
460 ELSE
461*
462* Store L(k) in column K
463*
464 d11 = a( k, k )
465 DO 16 ii = 1, k - 1
466 a( ii, k ) = a( ii, k ) / d11
467 16 CONTINUE
468*
469* Perform a rank-1 update of A(k+1:n,k+1:n) as
470* A := A - U(k)*D(k)*U(k)**T
471* = A - W(k)*(1/D(k))*W(k)**T
472* = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
473*
474 CALL zsyr( uplo, k-1, -d11, a( 1, k ), 1, a,
475 $ lda )
476 END IF
477 END IF
478*
479 ELSE
480*
481* 2-by-2 pivot block D(k): columns k and k-1 now hold
482*
483* ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
484*
485* where U(k) and U(k-1) are the k-th and (k-1)-th columns
486* of U
487*
488* Perform a rank-2 update of A(1:k-2,1:k-2) as
489*
490* A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T
491* = A - ( ( A(k-1)A(k) )*inv(D(k)) ) * ( A(k-1)A(k) )**T
492*
493* and store L(k) and L(k+1) in columns k and k+1
494*
495 IF( k.GT.2 ) THEN
496*
497 d12 = a( k-1, k )
498 d22 = a( k-1, k-1 ) / d12
499 d11 = a( k, k ) / d12
500 t = cone / ( d11*d22-cone )
501*
502 DO 30 j = k - 2, 1, -1
503*
504 wkm1 = t*( d11*a( j, k-1 )-a( j, k ) )
505 wk = t*( d22*a( j, k )-a( j, k-1 ) )
506*
507 DO 20 i = j, 1, -1
508 a( i, j ) = a( i, j ) - (a( i, k ) / d12 )*wk -
509 $ ( a( i, k-1 ) / d12 )*wkm1
510 20 CONTINUE
511*
512* Store U(k) and U(k-1) in cols k and k-1 for row J
513*
514 a( j, k ) = wk / d12
515 a( j, k-1 ) = wkm1 / d12
516*
517 30 CONTINUE
518*
519 END IF
520*
521 END IF
522 END IF
523*
524* Store details of the interchanges in IPIV
525*
526 IF( kstep.EQ.1 ) THEN
527 ipiv( k ) = kp
528 ELSE
529 ipiv( k ) = -p
530 ipiv( k-1 ) = -kp
531 END IF
532*
533* Decrease K and return to the start of the main loop
534*
535 k = k - kstep
536 GO TO 10
537*
538 ELSE
539*
540* Factorize A as L*D*L**T using the lower triangle of A
541*
542* K is the main loop index, increasing from 1 to N in steps of
543* 1 or 2
544*
545 k = 1
546 40 CONTINUE
547*
548* If K > N, exit from loop
549*
550 IF( k.GT.n )
551 $ GO TO 70
552 kstep = 1
553 p = k
554*
555* Determine rows and columns to be interchanged and whether
556* a 1-by-1 or 2-by-2 pivot block will be used
557*
558 absakk = cabs1( a( k, k ) )
559*
560* IMAX is the row-index of the largest off-diagonal element in
561* column K, and COLMAX is its absolute value.
562* Determine both COLMAX and IMAX.
563*
564 IF( k.LT.n ) THEN
565 imax = k + izamax( n-k, a( k+1, k ), 1 )
566 colmax = cabs1( a( imax, k ) )
567 ELSE
568 colmax = zero
569 END IF
570*
571 IF( ( max( absakk, colmax ).EQ.zero ) ) THEN
572*
573* Column K is zero or underflow: set INFO and continue
574*
575 IF( info.EQ.0 )
576 $ info = k
577 kp = k
578 ELSE
579*
580* Test for interchange
581*
582* Equivalent to testing for (used to handle NaN and Inf)
583* ABSAKK.GE.ALPHA*COLMAX
584*
585 IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
586*
587* no interchange, use 1-by-1 pivot block
588*
589 kp = k
590 ELSE
591*
592 done = .false.
593*
594* Loop until pivot found
595*
596 42 CONTINUE
597*
598* Begin pivot search loop body
599*
600* JMAX is the column-index of the largest off-diagonal
601* element in row IMAX, and ROWMAX is its absolute value.
602* Determine both ROWMAX and JMAX.
603*
604 IF( imax.NE.k ) THEN
605 jmax = k - 1 + izamax( imax-k, a( imax, k ),
606 $ lda )
607 rowmax = cabs1( a( imax, jmax ) )
608 ELSE
609 rowmax = zero
610 END IF
611*
612 IF( imax.LT.n ) THEN
613 itemp = imax + izamax( n-imax, a( imax+1,
614 $ imax ),
615 $ 1 )
616 dtemp = cabs1( a( itemp, imax ) )
617 IF( dtemp.GT.rowmax ) THEN
618 rowmax = dtemp
619 jmax = itemp
620 END IF
621 END IF
622*
623* Equivalent to testing for (used to handle NaN and Inf)
624* CABS1( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX
625*
626 IF( .NOT.( cabs1(a( imax, imax )).LT.alpha*rowmax ) )
627 $ THEN
628*
629* interchange rows and columns K and IMAX,
630* use 1-by-1 pivot block
631*
632 kp = imax
633 done = .true.
634*
635* Equivalent to testing for ROWMAX .EQ. COLMAX,
636* used to handle NaN and Inf
637*
638 ELSE IF( ( p.EQ.jmax ).OR.( rowmax.LE.colmax ) ) THEN
639*
640* interchange rows and columns K+1 and IMAX,
641* use 2-by-2 pivot block
642*
643 kp = imax
644 kstep = 2
645 done = .true.
646 ELSE
647*
648* Pivot NOT found, set variables and repeat
649*
650 p = imax
651 colmax = rowmax
652 imax = jmax
653 END IF
654*
655* End pivot search loop body
656*
657 IF( .NOT. done ) GOTO 42
658*
659 END IF
660*
661* Swap TWO rows and TWO columns
662*
663* First swap
664*
665 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
666*
667* Interchange rows and column K and P in the trailing
668* submatrix A(k:n,k:n) if we have a 2-by-2 pivot
669*
670 IF( p.LT.n )
671 $ CALL zswap( n-p, a( p+1, k ), 1, a( p+1, p ), 1 )
672 IF( p.GT.(k+1) )
673 $ CALL zswap( p-k-1, a( k+1, k ), 1, a( p, k+1 ),
674 $ lda )
675 t = a( k, k )
676 a( k, k ) = a( p, p )
677 a( p, p ) = t
678 END IF
679*
680* Second swap
681*
682 kk = k + kstep - 1
683 IF( kp.NE.kk ) THEN
684*
685* Interchange rows and columns KK and KP in the trailing
686* submatrix A(k:n,k:n)
687*
688 IF( kp.LT.n )
689 $ CALL zswap( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ),
690 $ 1 )
691 IF( ( kk.LT.n ) .AND. ( kp.GT.(kk+1) ) )
692 $ CALL zswap( kp-kk-1, a( kk+1, kk ), 1, a( kp,
693 $ kk+1 ),
694 $ lda )
695 t = a( kk, kk )
696 a( kk, kk ) = a( kp, kp )
697 a( kp, kp ) = t
698 IF( kstep.EQ.2 ) THEN
699 t = a( k+1, k )
700 a( k+1, k ) = a( kp, k )
701 a( kp, k ) = t
702 END IF
703 END IF
704*
705* Update the trailing submatrix
706*
707 IF( kstep.EQ.1 ) THEN
708*
709* 1-by-1 pivot block D(k): column k now holds
710*
711* W(k) = L(k)*D(k)
712*
713* where L(k) is the k-th column of L
714*
715 IF( k.LT.n ) THEN
716*
717* Perform a rank-1 update of A(k+1:n,k+1:n) and
718* store L(k) in column k
719*
720 IF( cabs1( a( k, k ) ).GE.sfmin ) THEN
721*
722* Perform a rank-1 update of A(k+1:n,k+1:n) as
723* A := A - L(k)*D(k)*L(k)**T
724* = A - W(k)*(1/D(k))*W(k)**T
725*
726 d11 = cone / a( k, k )
727 CALL zsyr( uplo, n-k, -d11, a( k+1, k ), 1,
728 $ a( k+1, k+1 ), lda )
729*
730* Store L(k) in column k
731*
732 CALL zscal( n-k, d11, a( k+1, k ), 1 )
733 ELSE
734*
735* Store L(k) in column k
736*
737 d11 = a( k, k )
738 DO 46 ii = k + 1, n
739 a( ii, k ) = a( ii, k ) / d11
740 46 CONTINUE
741*
742* Perform a rank-1 update of A(k+1:n,k+1:n) as
743* A := A - L(k)*D(k)*L(k)**T
744* = A - W(k)*(1/D(k))*W(k)**T
745* = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
746*
747 CALL zsyr( uplo, n-k, -d11, a( k+1, k ), 1,
748 $ a( k+1, k+1 ), lda )
749 END IF
750 END IF
751*
752 ELSE
753*
754* 2-by-2 pivot block D(k): columns k and k+1 now hold
755*
756* ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
757*
758* where L(k) and L(k+1) are the k-th and (k+1)-th columns
759* of L
760*
761*
762* Perform a rank-2 update of A(k+2:n,k+2:n) as
763*
764* A := A - ( L(k) L(k+1) ) * D(k) * ( L(k) L(k+1) )**T
765* = A - ( ( A(k)A(k+1) )*inv(D(k) ) * ( A(k)A(k+1) )**T
766*
767* and store L(k) and L(k+1) in columns k and k+1
768*
769 IF( k.LT.n-1 ) THEN
770*
771 d21 = a( k+1, k )
772 d11 = a( k+1, k+1 ) / d21
773 d22 = a( k, k ) / d21
774 t = cone / ( d11*d22-cone )
775*
776 DO 60 j = k + 2, n
777*
778* Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J
779*
780 wk = t*( d11*a( j, k )-a( j, k+1 ) )
781 wkp1 = t*( d22*a( j, k+1 )-a( j, k ) )
782*
783* Perform a rank-2 update of A(k+2:n,k+2:n)
784*
785 DO 50 i = j, n
786 a( i, j ) = a( i, j ) - ( a( i, k ) / d21 )*wk -
787 $ ( a( i, k+1 ) / d21 )*wkp1
788 50 CONTINUE
789*
790* Store L(k) and L(k+1) in cols k and k+1 for row J
791*
792 a( j, k ) = wk / d21
793 a( j, k+1 ) = wkp1 / d21
794*
795 60 CONTINUE
796*
797 END IF
798*
799 END IF
800 END IF
801*
802* Store details of the interchanges in IPIV
803*
804 IF( kstep.EQ.1 ) THEN
805 ipiv( k ) = kp
806 ELSE
807 ipiv( k ) = -p
808 ipiv( k+1 ) = -kp
809 END IF
810*
811* Increase K and return to the start of the main loop
812*
813 k = k + kstep
814 GO TO 40
815*
816 END IF
817*
818 70 CONTINUE
819*
820 RETURN
821*
822* End of ZSYTF2_ROOK
823*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine zsyr(uplo, n, alpha, x, incx, a, lda)
ZSYR performs the symmetric rank-1 update of a complex symmetric matrix.
Definition zsyr.f:133
integer function izamax(n, zx, incx)
IZAMAX
Definition izamax.f:71
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine zscal(n, za, zx, incx)
ZSCAL
Definition zscal.f:78
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP
Definition zswap.f:81
Here is the call graph for this function:
Here is the caller graph for this function: