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

◆ csytf2_rk()

subroutine csytf2_rk ( character uplo,
integer n,
complex, dimension( lda, * ) a,
integer lda,
complex, dimension( * ) e,
integer, dimension( * ) ipiv,
integer info )

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

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

Purpose:
!> CSYTF2_RK computes the factorization of a complex symmetric matrix A
!> using the bounded Bunch-Kaufman (rook) diagonal pivoting method:
!>
!>    A = P*U*D*(U**T)*(P**T) or A = P*L*D*(L**T)*(P**T),
!>
!> where U (or L) is unit upper (or lower) triangular matrix,
!> U**T (or L**T) is the transpose of U (or L), P is a permutation
!> matrix, P**T is the transpose of P, 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.
!> For more information see Further Details section.
!> 
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 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, contains:
!>            a) ONLY diagonal elements of the symmetric block diagonal
!>               matrix D on the diagonal of A, i.e. D(k,k) = A(k,k);
!>               (superdiagonal (or subdiagonal) elements of D
!>                are stored on exit in array E), and
!>            b) If UPLO = 'U': factor U in the superdiagonal part of A.
!>               If UPLO = 'L': factor L in the subdiagonal part of A.
!> 
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!> 
[out]E
!>          E is COMPLEX array, dimension (N)
!>          On exit, contains the superdiagonal (or subdiagonal)
!>          elements of the symmetric block diagonal matrix D
!>          with 1-by-1 or 2-by-2 diagonal blocks, where
!>          If UPLO = 'U': E(i) = D(i-1,i), i=2:N, E(1) is set to 0;
!>          If UPLO = 'L': E(i) = D(i+1,i), i=1:N-1, E(N) is set to 0.
!>
!>          NOTE: For 1-by-1 diagonal block D(k), where
!>          1 <= k <= N, the element E(k) is set to 0 in both
!>          UPLO = 'U' or UPLO = 'L' cases.
!> 
[out]IPIV
!>          IPIV is INTEGER array, dimension (N)
!>          IPIV describes the permutation matrix P in the factorization
!>          of matrix A as follows. The absolute value of IPIV(k)
!>          represents the index of row and column that were
!>          interchanged with the k-th row and column. The value of UPLO
!>          describes the order in which the interchanges were applied.
!>          Also, the sign of IPIV represents the block structure of
!>          the symmetric block diagonal matrix D with 1-by-1 or 2-by-2
!>          diagonal blocks which correspond to 1 or 2 interchanges
!>          at each factorization step. For more info see Further
!>          Details section.
!>
!>          If UPLO = 'U',
!>          ( in factorization order, k decreases from N to 1 ):
!>            a) A single positive entry IPIV(k) > 0 means:
!>               D(k,k) is a 1-by-1 diagonal block.
!>               If IPIV(k) != k, rows and columns k and IPIV(k) were
!>               interchanged in the matrix A(1:N,1:N);
!>               If IPIV(k) = k, no interchange occurred.
!>
!>            b) A pair of consecutive negative entries
!>               IPIV(k) < 0 and IPIV(k-1) < 0 means:
!>               D(k-1:k,k-1:k) is a 2-by-2 diagonal block.
!>               (NOTE: negative entries in IPIV appear ONLY in pairs).
!>               1) If -IPIV(k) != k, rows and columns
!>                  k and -IPIV(k) were interchanged
!>                  in the matrix A(1:N,1:N).
!>                  If -IPIV(k) = k, no interchange occurred.
!>               2) If -IPIV(k-1) != k-1, rows and columns
!>                  k-1 and -IPIV(k-1) were interchanged
!>                  in the matrix A(1:N,1:N).
!>                  If -IPIV(k-1) = k-1, no interchange occurred.
!>
!>            c) In both cases a) and b), always ABS( IPIV(k) ) <= k.
!>
!>            d) NOTE: Any entry IPIV(k) is always NONZERO on output.
!>
!>          If UPLO = 'L',
!>          ( in factorization order, k increases from 1 to N ):
!>            a) A single positive entry IPIV(k) > 0 means:
!>               D(k,k) is a 1-by-1 diagonal block.
!>               If IPIV(k) != k, rows and columns k and IPIV(k) were
!>               interchanged in the matrix A(1:N,1:N).
!>               If IPIV(k) = k, no interchange occurred.
!>
!>            b) A pair of consecutive negative entries
!>               IPIV(k) < 0 and IPIV(k+1) < 0 means:
!>               D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
!>               (NOTE: negative entries in IPIV appear ONLY in pairs).
!>               1) If -IPIV(k) != k, rows and columns
!>                  k and -IPIV(k) were interchanged
!>                  in the matrix A(1:N,1:N).
!>                  If -IPIV(k) = k, no interchange occurred.
!>               2) If -IPIV(k+1) != k+1, rows and columns
!>                  k-1 and -IPIV(k-1) were interchanged
!>                  in the matrix A(1:N,1:N).
!>                  If -IPIV(k+1) = k+1, no interchange occurred.
!>
!>            c) In both cases a) and b), always ABS( IPIV(k) ) >= k.
!>
!>            d) NOTE: Any entry IPIV(k) is always NONZERO on output.
!> 
[out]INFO
!>          INFO is INTEGER
!>          = 0: successful exit
!>
!>          < 0: If INFO = -k, the k-th argument had an illegal value
!>
!>          > 0: If INFO = k, the matrix A is singular, because:
!>                 If UPLO = 'U': column k in the upper
!>                 triangular part of A contains all zeros.
!>                 If UPLO = 'L': column k in the lower
!>                 triangular part of A contains all zeros.
!>
!>               Therefore D(k,k) is exactly zero, and superdiagonal
!>               elements of column k of U (or subdiagonal elements of
!>               column k of L ) are all zeros. 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.
!>
!>               NOTE: INFO only stores the first occurrence of
!>               a singularity, any subsequent occurrence of singularity
!>               is not stored in INFO even though the factorization
!>               always completes.
!> 
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!> TODO: put further details
!> 
Contributors:
!>
!>  December 2016,  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 238 of file csytf2_rk.f.

239*
240* -- LAPACK computational routine --
241* -- LAPACK is a software package provided by Univ. of Tennessee, --
242* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
243*
244* .. Scalar Arguments ..
245 CHARACTER UPLO
246 INTEGER INFO, LDA, N
247* ..
248* .. Array Arguments ..
249 INTEGER IPIV( * )
250 COMPLEX A( LDA, * ), E( * )
251* ..
252*
253* =====================================================================
254*
255* .. Parameters ..
256 REAL ZERO, ONE
257 parameter( zero = 0.0e+0, one = 1.0e+0 )
258 REAL EIGHT, SEVTEN
259 parameter( eight = 8.0e+0, sevten = 17.0e+0 )
260 COMPLEX CONE, CZERO
261 parameter( cone = ( 1.0e+0, 0.0e+0 ),
262 $ czero = ( 0.0e+0, 0.0e+0 ) )
263* ..
264* .. Local Scalars ..
265 LOGICAL UPPER, DONE
266 INTEGER I, IMAX, J, JMAX, ITEMP, K, KK, KP, KSTEP,
267 $ P, II
268 REAL ABSAKK, ALPHA, COLMAX, ROWMAX, STEMP, SFMIN
269 COMPLEX D11, D12, D21, D22, T, WK, WKM1, WKP1, Z
270* ..
271* .. External Functions ..
272 LOGICAL LSAME
273 INTEGER ICAMAX
274 REAL SLAMCH
275 EXTERNAL lsame, icamax, slamch
276* ..
277* .. External Subroutines ..
278 EXTERNAL cscal, cswap, csyr, xerbla
279* ..
280* .. Intrinsic Functions ..
281 INTRINSIC abs, max, sqrt, aimag, real
282* ..
283* .. Statement Functions ..
284 REAL CABS1
285* ..
286* .. Statement Function definitions ..
287 cabs1( z ) = abs( real( z ) ) + abs( aimag( z ) )
288* ..
289* .. Executable Statements ..
290*
291* Test the input parameters.
292*
293 info = 0
294 upper = lsame( uplo, 'U' )
295 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
296 info = -1
297 ELSE IF( n.LT.0 ) THEN
298 info = -2
299 ELSE IF( lda.LT.max( 1, n ) ) THEN
300 info = -4
301 END IF
302 IF( info.NE.0 ) THEN
303 CALL xerbla( 'CSYTF2_RK', -info )
304 RETURN
305 END IF
306*
307* Initialize ALPHA for use in choosing pivot block size.
308*
309 alpha = ( one+sqrt( sevten ) ) / eight
310*
311* Compute machine safe minimum
312*
313 sfmin = slamch( 'S' )
314*
315 IF( upper ) THEN
316*
317* Factorize A as U*D*U**T using the upper triangle of A
318*
319* Initialize the first entry of array E, where superdiagonal
320* elements of D are stored
321*
322 e( 1 ) = czero
323*
324* K is the main loop index, decreasing from N to 1 in steps of
325* 1 or 2
326*
327 k = n
328 10 CONTINUE
329*
330* If K < 1, exit from loop
331*
332 IF( k.LT.1 )
333 $ GO TO 34
334 kstep = 1
335 p = k
336*
337* Determine rows and columns to be interchanged and whether
338* a 1-by-1 or 2-by-2 pivot block will be used
339*
340 absakk = cabs1( a( k, k ) )
341*
342* IMAX is the row-index of the largest off-diagonal element in
343* column K, and COLMAX is its absolute value.
344* Determine both COLMAX and IMAX.
345*
346 IF( k.GT.1 ) THEN
347 imax = icamax( k-1, a( 1, k ), 1 )
348 colmax = cabs1( a( imax, k ) )
349 ELSE
350 colmax = zero
351 END IF
352*
353 IF( (max( absakk, colmax ).EQ.zero) ) THEN
354*
355* Column K is zero or underflow: set INFO and continue
356*
357 IF( info.EQ.0 )
358 $ info = k
359 kp = k
360*
361* Set E( K ) to zero
362*
363 IF( k.GT.1 )
364 $ e( k ) = czero
365*
366 ELSE
367*
368* Test for interchange
369*
370* Equivalent to testing for (used to handle NaN and Inf)
371* ABSAKK.GE.ALPHA*COLMAX
372*
373 IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
374*
375* no interchange,
376* use 1-by-1 pivot block
377*
378 kp = k
379 ELSE
380*
381 done = .false.
382*
383* Loop until pivot found
384*
385 12 CONTINUE
386*
387* Begin pivot search loop body
388*
389* JMAX is the column-index of the largest off-diagonal
390* element in row IMAX, and ROWMAX is its absolute value.
391* Determine both ROWMAX and JMAX.
392*
393 IF( imax.NE.k ) THEN
394 jmax = imax + icamax( k-imax, a( imax, imax+1 ),
395 $ lda )
396 rowmax = cabs1( a( imax, jmax ) )
397 ELSE
398 rowmax = zero
399 END IF
400*
401 IF( imax.GT.1 ) THEN
402 itemp = icamax( imax-1, a( 1, imax ), 1 )
403 stemp = cabs1( a( itemp, imax ) )
404 IF( stemp.GT.rowmax ) THEN
405 rowmax = stemp
406 jmax = itemp
407 END IF
408 END IF
409*
410* Equivalent to testing for (used to handle NaN and Inf)
411* ABS( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX
412*
413 IF( .NOT.( cabs1( a( imax, imax ) ).LT.alpha*rowmax ))
414 $ THEN
415*
416* interchange rows and columns K and IMAX,
417* use 1-by-1 pivot block
418*
419 kp = imax
420 done = .true.
421*
422* Equivalent to testing for ROWMAX .EQ. COLMAX,
423* used to handle NaN and Inf
424*
425 ELSE IF( ( p.EQ.jmax ).OR.( rowmax.LE.colmax ) ) THEN
426*
427* interchange rows and columns K+1 and IMAX,
428* use 2-by-2 pivot block
429*
430 kp = imax
431 kstep = 2
432 done = .true.
433 ELSE
434*
435* Pivot NOT found, set variables and repeat
436*
437 p = imax
438 colmax = rowmax
439 imax = jmax
440 END IF
441*
442* End pivot search loop body
443*
444 IF( .NOT. done ) GOTO 12
445*
446 END IF
447*
448* Swap TWO rows and TWO columns
449*
450* First swap
451*
452 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
453*
454* Interchange rows and column K and P in the leading
455* submatrix A(1:k,1:k) if we have a 2-by-2 pivot
456*
457 IF( p.GT.1 )
458 $ CALL cswap( p-1, a( 1, k ), 1, a( 1, p ), 1 )
459 IF( p.LT.(k-1) )
460 $ CALL cswap( k-p-1, a( p+1, k ), 1, a( p, p+1 ),
461 $ lda )
462 t = a( k, k )
463 a( k, k ) = a( p, p )
464 a( p, p ) = t
465*
466* Convert upper triangle of A into U form by applying
467* the interchanges in columns k+1:N.
468*
469 IF( k.LT.n )
470 $ CALL cswap( n-k, a( k, k+1 ), lda, a( p, k+1 ),
471 $ lda )
472*
473 END IF
474*
475* Second swap
476*
477 kk = k - kstep + 1
478 IF( kp.NE.kk ) THEN
479*
480* Interchange rows and columns KK and KP in the leading
481* submatrix A(1:k,1:k)
482*
483 IF( kp.GT.1 )
484 $ CALL cswap( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
485 IF( ( kk.GT.1 ) .AND. ( kp.LT.(kk-1) ) )
486 $ CALL cswap( kk-kp-1, a( kp+1, kk ), 1, a( kp,
487 $ kp+1 ),
488 $ lda )
489 t = a( kk, kk )
490 a( kk, kk ) = a( kp, kp )
491 a( kp, kp ) = t
492 IF( kstep.EQ.2 ) THEN
493 t = a( k-1, k )
494 a( k-1, k ) = a( kp, k )
495 a( kp, k ) = t
496 END IF
497*
498* Convert upper triangle of A into U form by applying
499* the interchanges in columns k+1:N.
500*
501 IF( k.LT.n )
502 $ CALL cswap( n-k, a( kk, k+1 ), lda, a( kp, k+1 ),
503 $ lda )
504*
505 END IF
506*
507* Update the leading submatrix
508*
509 IF( kstep.EQ.1 ) THEN
510*
511* 1-by-1 pivot block D(k): column k now holds
512*
513* W(k) = U(k)*D(k)
514*
515* where U(k) is the k-th column of U
516*
517 IF( k.GT.1 ) THEN
518*
519* Perform a rank-1 update of A(1:k-1,1:k-1) and
520* store U(k) in column k
521*
522 IF( cabs1( a( k, k ) ).GE.sfmin ) THEN
523*
524* Perform a rank-1 update of A(1:k-1,1:k-1) as
525* A := A - U(k)*D(k)*U(k)**T
526* = A - W(k)*1/D(k)*W(k)**T
527*
528 d11 = cone / a( k, k )
529 CALL csyr( uplo, k-1, -d11, a( 1, k ), 1, a,
530 $ lda )
531*
532* Store U(k) in column k
533*
534 CALL cscal( k-1, d11, a( 1, k ), 1 )
535 ELSE
536*
537* Store L(k) in column K
538*
539 d11 = a( k, k )
540 DO 16 ii = 1, k - 1
541 a( ii, k ) = a( ii, k ) / d11
542 16 CONTINUE
543*
544* Perform a rank-1 update of A(k+1:n,k+1:n) as
545* A := A - U(k)*D(k)*U(k)**T
546* = A - W(k)*(1/D(k))*W(k)**T
547* = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
548*
549 CALL csyr( uplo, k-1, -d11, a( 1, k ), 1, a,
550 $ lda )
551 END IF
552*
553* Store the superdiagonal element of D in array E
554*
555 e( k ) = czero
556*
557 END IF
558*
559 ELSE
560*
561* 2-by-2 pivot block D(k): columns k and k-1 now hold
562*
563* ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
564*
565* where U(k) and U(k-1) are the k-th and (k-1)-th columns
566* of U
567*
568* Perform a rank-2 update of A(1:k-2,1:k-2) as
569*
570* A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T
571* = A - ( ( A(k-1)A(k) )*inv(D(k)) ) * ( A(k-1)A(k) )**T
572*
573* and store L(k) and L(k+1) in columns k and k+1
574*
575 IF( k.GT.2 ) THEN
576*
577 d12 = a( k-1, k )
578 d22 = a( k-1, k-1 ) / d12
579 d11 = a( k, k ) / d12
580 t = cone / ( d11*d22-cone )
581*
582 DO 30 j = k - 2, 1, -1
583*
584 wkm1 = t*( d11*a( j, k-1 )-a( j, k ) )
585 wk = t*( d22*a( j, k )-a( j, k-1 ) )
586*
587 DO 20 i = j, 1, -1
588 a( i, j ) = a( i, j ) - (a( i, k ) / d12 )*wk -
589 $ ( a( i, k-1 ) / d12 )*wkm1
590 20 CONTINUE
591*
592* Store U(k) and U(k-1) in cols k and k-1 for row J
593*
594 a( j, k ) = wk / d12
595 a( j, k-1 ) = wkm1 / d12
596*
597 30 CONTINUE
598*
599 END IF
600*
601* Copy superdiagonal elements of D(K) to E(K) and
602* ZERO out superdiagonal entry of A
603*
604 e( k ) = a( k-1, k )
605 e( k-1 ) = czero
606 a( k-1, k ) = czero
607*
608 END IF
609*
610* End column K is nonsingular
611*
612 END IF
613*
614* Store details of the interchanges in IPIV
615*
616 IF( kstep.EQ.1 ) THEN
617 ipiv( k ) = kp
618 ELSE
619 ipiv( k ) = -p
620 ipiv( k-1 ) = -kp
621 END IF
622*
623* Decrease K and return to the start of the main loop
624*
625 k = k - kstep
626 GO TO 10
627*
628 34 CONTINUE
629*
630 ELSE
631*
632* Factorize A as L*D*L**T using the lower triangle of A
633*
634* Initialize the unused last entry of the subdiagonal array E.
635*
636 e( n ) = czero
637*
638* K is the main loop index, increasing from 1 to N in steps of
639* 1 or 2
640*
641 k = 1
642 40 CONTINUE
643*
644* If K > N, exit from loop
645*
646 IF( k.GT.n )
647 $ GO TO 64
648 kstep = 1
649 p = k
650*
651* Determine rows and columns to be interchanged and whether
652* a 1-by-1 or 2-by-2 pivot block will be used
653*
654 absakk = cabs1( a( k, k ) )
655*
656* IMAX is the row-index of the largest off-diagonal element in
657* column K, and COLMAX is its absolute value.
658* Determine both COLMAX and IMAX.
659*
660 IF( k.LT.n ) THEN
661 imax = k + icamax( n-k, a( k+1, k ), 1 )
662 colmax = cabs1( a( imax, k ) )
663 ELSE
664 colmax = zero
665 END IF
666*
667 IF( ( max( absakk, colmax ).EQ.zero ) ) THEN
668*
669* Column K is zero or underflow: set INFO and continue
670*
671 IF( info.EQ.0 )
672 $ info = k
673 kp = k
674*
675* Set E( K ) to zero
676*
677 IF( k.LT.n )
678 $ e( k ) = czero
679*
680 ELSE
681*
682* Test for interchange
683*
684* Equivalent to testing for (used to handle NaN and Inf)
685* ABSAKK.GE.ALPHA*COLMAX
686*
687 IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
688*
689* no interchange, use 1-by-1 pivot block
690*
691 kp = k
692*
693 ELSE
694*
695 done = .false.
696*
697* Loop until pivot found
698*
699 42 CONTINUE
700*
701* Begin pivot search loop body
702*
703* JMAX is the column-index of the largest off-diagonal
704* element in row IMAX, and ROWMAX is its absolute value.
705* Determine both ROWMAX and JMAX.
706*
707 IF( imax.NE.k ) THEN
708 jmax = k - 1 + icamax( imax-k, a( imax, k ),
709 $ lda )
710 rowmax = cabs1( a( imax, jmax ) )
711 ELSE
712 rowmax = zero
713 END IF
714*
715 IF( imax.LT.n ) THEN
716 itemp = imax + icamax( n-imax, a( imax+1,
717 $ imax ),
718 $ 1 )
719 stemp = cabs1( a( itemp, imax ) )
720 IF( stemp.GT.rowmax ) THEN
721 rowmax = stemp
722 jmax = itemp
723 END IF
724 END IF
725*
726* Equivalent to testing for (used to handle NaN and Inf)
727* ABS( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX
728*
729 IF( .NOT.( cabs1( a( imax, imax ) ).LT.alpha*rowmax ))
730 $ THEN
731*
732* interchange rows and columns K and IMAX,
733* use 1-by-1 pivot block
734*
735 kp = imax
736 done = .true.
737*
738* Equivalent to testing for ROWMAX .EQ. COLMAX,
739* used to handle NaN and Inf
740*
741 ELSE IF( ( p.EQ.jmax ).OR.( rowmax.LE.colmax ) ) THEN
742*
743* interchange rows and columns K+1 and IMAX,
744* use 2-by-2 pivot block
745*
746 kp = imax
747 kstep = 2
748 done = .true.
749 ELSE
750*
751* Pivot NOT found, set variables and repeat
752*
753 p = imax
754 colmax = rowmax
755 imax = jmax
756 END IF
757*
758* End pivot search loop body
759*
760 IF( .NOT. done ) GOTO 42
761*
762 END IF
763*
764* Swap TWO rows and TWO columns
765*
766* First swap
767*
768 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
769*
770* Interchange rows and column K and P in the trailing
771* submatrix A(k:n,k:n) if we have a 2-by-2 pivot
772*
773 IF( p.LT.n )
774 $ CALL cswap( n-p, a( p+1, k ), 1, a( p+1, p ), 1 )
775 IF( p.GT.(k+1) )
776 $ CALL cswap( p-k-1, a( k+1, k ), 1, a( p, k+1 ),
777 $ lda )
778 t = a( k, k )
779 a( k, k ) = a( p, p )
780 a( p, p ) = t
781*
782* Convert lower triangle of A into L form by applying
783* the interchanges in columns 1:k-1.
784*
785 IF ( k.GT.1 )
786 $ CALL cswap( k-1, a( k, 1 ), lda, a( p, 1 ), lda )
787*
788 END IF
789*
790* Second swap
791*
792 kk = k + kstep - 1
793 IF( kp.NE.kk ) THEN
794*
795* Interchange rows and columns KK and KP in the trailing
796* submatrix A(k:n,k:n)
797*
798 IF( kp.LT.n )
799 $ CALL cswap( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ),
800 $ 1 )
801 IF( ( kk.LT.n ) .AND. ( kp.GT.(kk+1) ) )
802 $ CALL cswap( kp-kk-1, a( kk+1, kk ), 1, a( kp,
803 $ kk+1 ),
804 $ lda )
805 t = a( kk, kk )
806 a( kk, kk ) = a( kp, kp )
807 a( kp, kp ) = t
808 IF( kstep.EQ.2 ) THEN
809 t = a( k+1, k )
810 a( k+1, k ) = a( kp, k )
811 a( kp, k ) = t
812 END IF
813*
814* Convert lower triangle of A into L form by applying
815* the interchanges in columns 1:k-1.
816*
817 IF ( k.GT.1 )
818 $ CALL cswap( k-1, a( kk, 1 ), lda, a( kp, 1 ), lda )
819*
820 END IF
821*
822* Update the trailing submatrix
823*
824 IF( kstep.EQ.1 ) THEN
825*
826* 1-by-1 pivot block D(k): column k now holds
827*
828* W(k) = L(k)*D(k)
829*
830* where L(k) is the k-th column of L
831*
832 IF( k.LT.n ) THEN
833*
834* Perform a rank-1 update of A(k+1:n,k+1:n) and
835* store L(k) in column k
836*
837 IF( cabs1( a( k, k ) ).GE.sfmin ) THEN
838*
839* Perform a rank-1 update of A(k+1:n,k+1:n) as
840* A := A - L(k)*D(k)*L(k)**T
841* = A - W(k)*(1/D(k))*W(k)**T
842*
843 d11 = cone / a( k, k )
844 CALL csyr( uplo, n-k, -d11, a( k+1, k ), 1,
845 $ a( k+1, k+1 ), lda )
846*
847* Store L(k) in column k
848*
849 CALL cscal( n-k, d11, a( k+1, k ), 1 )
850 ELSE
851*
852* Store L(k) in column k
853*
854 d11 = a( k, k )
855 DO 46 ii = k + 1, n
856 a( ii, k ) = a( ii, k ) / d11
857 46 CONTINUE
858*
859* Perform a rank-1 update of A(k+1:n,k+1:n) as
860* A := A - L(k)*D(k)*L(k)**T
861* = A - W(k)*(1/D(k))*W(k)**T
862* = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
863*
864 CALL csyr( uplo, n-k, -d11, a( k+1, k ), 1,
865 $ a( k+1, k+1 ), lda )
866 END IF
867*
868* Store the subdiagonal element of D in array E
869*
870 e( k ) = czero
871*
872 END IF
873*
874 ELSE
875*
876* 2-by-2 pivot block D(k): columns k and k+1 now hold
877*
878* ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
879*
880* where L(k) and L(k+1) are the k-th and (k+1)-th columns
881* of L
882*
883*
884* Perform a rank-2 update of A(k+2:n,k+2:n) as
885*
886* A := A - ( L(k) L(k+1) ) * D(k) * ( L(k) L(k+1) )**T
887* = A - ( ( A(k)A(k+1) )*inv(D(k) ) * ( A(k)A(k+1) )**T
888*
889* and store L(k) and L(k+1) in columns k and k+1
890*
891 IF( k.LT.n-1 ) THEN
892*
893 d21 = a( k+1, k )
894 d11 = a( k+1, k+1 ) / d21
895 d22 = a( k, k ) / d21
896 t = cone / ( d11*d22-cone )
897*
898 DO 60 j = k + 2, n
899*
900* Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J
901*
902 wk = t*( d11*a( j, k )-a( j, k+1 ) )
903 wkp1 = t*( d22*a( j, k+1 )-a( j, k ) )
904*
905* Perform a rank-2 update of A(k+2:n,k+2:n)
906*
907 DO 50 i = j, n
908 a( i, j ) = a( i, j ) - ( a( i, k ) / d21 )*wk -
909 $ ( a( i, k+1 ) / d21 )*wkp1
910 50 CONTINUE
911*
912* Store L(k) and L(k+1) in cols k and k+1 for row J
913*
914 a( j, k ) = wk / d21
915 a( j, k+1 ) = wkp1 / d21
916*
917 60 CONTINUE
918*
919 END IF
920*
921* Copy subdiagonal elements of D(K) to E(K) and
922* ZERO out subdiagonal entry of A
923*
924 e( k ) = a( k+1, k )
925 e( k+1 ) = czero
926 a( k+1, k ) = czero
927*
928 END IF
929*
930* End column K is nonsingular
931*
932 END IF
933*
934* Store details of the interchanges in IPIV
935*
936 IF( kstep.EQ.1 ) THEN
937 ipiv( k ) = kp
938 ELSE
939 ipiv( k ) = -p
940 ipiv( k+1 ) = -kp
941 END IF
942*
943* Increase K and return to the start of the main loop
944*
945 k = k + kstep
946 GO TO 40
947*
948 64 CONTINUE
949*
950 END IF
951*
952 RETURN
953*
954* End of CSYTF2_RK
955*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine csyr(uplo, n, alpha, x, incx, a, lda)
CSYR performs the symmetric rank-1 update of a complex symmetric matrix.
Definition csyr.f:133
integer function icamax(n, cx, incx)
ICAMAX
Definition icamax.f:71
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine cscal(n, ca, cx, incx)
CSCAL
Definition cscal.f:78
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81
Here is the call graph for this function:
Here is the caller graph for this function: