LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ zsytf2_rk()

subroutine zsytf2_rk ( character  UPLO,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( * )  E,
integer, dimension( * )  IPIV,
integer  INFO 
)

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

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

Purpose:
 ZSYTF2_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*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, 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*16 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 zsytf2_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*16 A( LDA, * ), E( * )
253 * ..
254 *
255 * =====================================================================
256 *
257 * .. Parameters ..
258  DOUBLE PRECISION ZERO, ONE
259  parameter( zero = 0.0d+0, one = 1.0d+0 )
260  DOUBLE PRECISION EIGHT, SEVTEN
261  parameter( eight = 8.0d+0, sevten = 17.0d+0 )
262  COMPLEX*16 CONE, CZERO
263  parameter( cone = ( 1.0d+0, 0.0d+0 ),
264  $ czero = ( 0.0d+0, 0.0d+0 ) )
265 * ..
266 * .. Local Scalars ..
267  LOGICAL UPPER, DONE
268  INTEGER I, IMAX, J, JMAX, ITEMP, K, KK, KP, KSTEP,
269  $ P, II
270  DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, ROWMAX, DTEMP, SFMIN
271  COMPLEX*16 D11, D12, D21, D22, T, WK, WKM1, WKP1, Z
272 * ..
273 * .. External Functions ..
274  LOGICAL LSAME
275  INTEGER IZAMAX
276  DOUBLE PRECISION DLAMCH
277  EXTERNAL lsame, izamax, dlamch
278 * ..
279 * .. External Subroutines ..
280  EXTERNAL zscal, zswap, zsyr, xerbla
281 * ..
282 * .. Intrinsic Functions ..
283  INTRINSIC abs, max, sqrt, dimag, dble
284 * ..
285 * .. Statement Functions ..
286  DOUBLE PRECISION CABS1
287 * ..
288 * .. Statement Function definitions ..
289  cabs1( z ) = abs( dble( z ) ) + abs( dimag( 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( 'ZSYTF2_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 = dlamch( '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 = izamax( 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 + izamax( 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 = izamax( imax-1, a( 1, imax ), 1 )
405  dtemp = cabs1( a( itemp, imax ) )
406  IF( dtemp.GT.rowmax ) THEN
407  rowmax = dtemp
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 zswap( p-1, a( 1, k ), 1, a( 1, p ), 1 )
461  IF( p.LT.(k-1) )
462  $ CALL zswap( 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 zswap( 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 zswap( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
486  IF( ( kk.GT.1 ) .AND. ( kp.LT.(kk-1) ) )
487  $ CALL zswap( 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 zswap( 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 zsyr( uplo, k-1, -d11, a( 1, k ), 1, a, lda )
530 *
531 * Store U(k) in column k
532 *
533  CALL zscal( 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 zsyr( 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 + izamax( 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 + izamax( 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 + izamax( n-imax, a( imax+1, imax ),
714  $ 1 )
715  dtemp = cabs1( a( itemp, imax ) )
716  IF( dtemp.GT.rowmax ) THEN
717  rowmax = dtemp
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 zswap( n-p, a( p+1, k ), 1, a( p+1, p ), 1 )
771  IF( p.GT.(k+1) )
772  $ CALL zswap( 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 zswap( 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 zswap( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
795  IF( ( kk.LT.n ) .AND. ( kp.GT.(kk+1) ) )
796  $ CALL zswap( 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 zswap( 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 zsyr( 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 zscal( 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 zsyr( 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 ZSYTF2_RK
948 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:71
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:81
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:78
subroutine zsyr(UPLO, N, ALPHA, X, INCX, A, LDA)
ZSYR performs the symmetric rank-1 update of a complex symmetric matrix.
Definition: zsyr.f:135
Here is the call graph for this function:
Here is the caller graph for this function: