LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ zsytf2_rook()

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

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

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

Purpose:
 ZSYTF2_ROOK computes the factorization of a complex symmetric matrix A
 using the bounded Bunch-Kaufman ("rook") diagonal pivoting method:

    A = U*D*U**T  or  A = L*D*L**T

 where U (or L) is a product of permutation and unit upper (lower)
 triangular matrices, U**T is the transpose of U, and D is symmetric and
 block diagonal with 1-by-1 and 2-by-2 diagonal blocks.

 This is the unblocked version of the algorithm, calling Level 2 BLAS.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          Specifies whether the upper or lower triangular part of the
          symmetric matrix A is stored:
          = 'U':  Upper triangular
          = 'L':  Lower triangular
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
          n-by-n upper triangular part of A contains the upper
          triangular part of the matrix A, and the strictly lower
          triangular part of A is not referenced.  If UPLO = 'L', the
          leading n-by-n lower triangular part of A contains the lower
          triangular part of the matrix A, and the strictly upper
          triangular part of A is not referenced.

          On exit, the block diagonal matrix D and the multipliers used
          to obtain the factor U or L (see below for further details).
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[out]IPIV
          IPIV is INTEGER array, dimension (N)
          Details of the interchanges and the block structure of D.

          If UPLO = 'U':
             If IPIV(k) > 0, then rows and columns k and IPIV(k)
             were interchanged and D(k,k) is a 1-by-1 diagonal block.

             If IPIV(k) < 0 and IPIV(k-1) < 0, then rows and
             columns k and -IPIV(k) were interchanged and rows and
             columns k-1 and -IPIV(k-1) were inerchaged,
             D(k-1:k,k-1:k) is a 2-by-2 diagonal block.

          If UPLO = 'L':
             If IPIV(k) > 0, then rows and columns k and IPIV(k)
             were interchanged and D(k,k) is a 1-by-1 diagonal block.

             If IPIV(k) < 0 and IPIV(k+1) < 0, then rows and
             columns k and -IPIV(k) were interchanged and rows and
             columns k+1 and -IPIV(k+1) were inerchaged,
             D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -k, the k-th argument had an illegal value
          > 0: if INFO = k, D(k,k) is exactly zero.  The factorization
               has been completed, but the block diagonal matrix D is
               exactly singular, and division by zero will occur if it
               is used to solve a system of equations.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2013
Further Details:
  If UPLO = 'U', then A = U*D*U**T, where
     U = P(n)*U(n)* ... *P(k)U(k)* ...,
  i.e., U is a product of terms P(k)*U(k), where k decreases from n to
  1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
  that if the diagonal block D(k) is of order s (s = 1 or 2), then

             (   I    v    0   )   k-s
     U(k) =  (   0    I    0   )   s
             (   0    0    I   )   n-k
                k-s   s   n-k

  If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
  If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
  and A(k,k), and v overwrites A(1:k-2,k-1:k).

  If UPLO = 'L', then A = L*D*L**T, where
     L = P(1)*L(1)* ... *P(k)*L(k)* ...,
  i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
  n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
  defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
  that if the diagonal block D(k) is of order s (s = 1 or 2), then

             (   I    0     0   )  k-1
     L(k) =  (   0    I     0   )  s
             (   0    v     I   )  n-k-s+1
                k-1   s  n-k-s+1

  If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
  If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
  and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
Contributors:
  November 2013,     Igor Kozachenko,
                  Computer Science Division,
                  University of California, Berkeley

  September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
                  School of Mathematics,
                  University of Manchester

  01-01-96 - Based on modifications by
    J. Lewis, Boeing Computer Services Company
    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville abd , USA

Definition at line 196 of file zsytf2_rook.f.

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