LAPACK 3.11.0
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 240 of file csytf2_rk.f.

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