LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ 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.
Date
December 2016
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 243 of file csytf2_rk.f.

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