LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine csytf2 ( character  UPLO,
integer  N,
complex, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
integer  INFO 
)

CSYTF2 computes the factorization of a real symmetric indefinite matrix, using the diagonal pivoting method (unblocked algorithm).

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

Purpose:
 CSYTF2 computes the factorization of a complex symmetric matrix A
 using the Bunch-Kaufman 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 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) = IPIV(k-1) < 0, then rows and columns
             k-1 and -IPIV(k) were interchanged and 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) = IPIV(k+1) < 0, then rows and columns
             k+1 and -IPIV(k) were interchanged and 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:
  09-29-06 - patch from
    Bobby Cheng, MathWorks

    Replace l.209 and l.377
         IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
    by
         IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. SISNAN(ABSAKK) ) THEN

  1-96 - Based on modifications by J. Lewis, Boeing Computer Services
         Company

Definition at line 193 of file csytf2.f.

193 *
194 * -- LAPACK computational routine (version 3.5.0) --
195 * -- LAPACK is a software package provided by Univ. of Tennessee, --
196 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
197 * November 2013
198 *
199 * .. Scalar Arguments ..
200  CHARACTER uplo
201  INTEGER info, lda, n
202 * ..
203 * .. Array Arguments ..
204  INTEGER ipiv( * )
205  COMPLEX a( lda, * )
206 * ..
207 *
208 * =====================================================================
209 *
210 * .. Parameters ..
211  REAL zero, one
212  parameter ( zero = 0.0e+0, one = 1.0e+0 )
213  REAL eight, sevten
214  parameter ( eight = 8.0e+0, sevten = 17.0e+0 )
215  COMPLEX cone
216  parameter ( cone = ( 1.0e+0, 0.0e+0 ) )
217 * ..
218 * .. Local Scalars ..
219  LOGICAL upper
220  INTEGER i, imax, j, jmax, k, kk, kp, kstep
221  REAL absakk, alpha, colmax, rowmax
222  COMPLEX d11, d12, d21, d22, r1, t, wk, wkm1, wkp1, z
223 * ..
224 * .. External Functions ..
225  LOGICAL lsame, sisnan
226  INTEGER icamax
227  EXTERNAL lsame, icamax, sisnan
228 * ..
229 * .. External Subroutines ..
230  EXTERNAL cscal, cswap, csyr, xerbla
231 * ..
232 * .. Intrinsic Functions ..
233  INTRINSIC abs, aimag, max, REAL, sqrt
234 * ..
235 * .. Statement Functions ..
236  REAL cabs1
237 * ..
238 * .. Statement Function definitions ..
239  cabs1( z ) = abs( REAL( Z ) ) + abs( aimag( z ) )
240 * ..
241 * .. Executable Statements ..
242 *
243 * Test the input parameters.
244 *
245  info = 0
246  upper = lsame( uplo, 'U' )
247  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
248  info = -1
249  ELSE IF( n.LT.0 ) THEN
250  info = -2
251  ELSE IF( lda.LT.max( 1, n ) ) THEN
252  info = -4
253  END IF
254  IF( info.NE.0 ) THEN
255  CALL xerbla( 'CSYTF2', -info )
256  RETURN
257  END IF
258 *
259 * Initialize ALPHA for use in choosing pivot block size.
260 *
261  alpha = ( one+sqrt( sevten ) ) / eight
262 *
263  IF( upper ) THEN
264 *
265 * Factorize A as U*D*U**T using the upper triangle of A
266 *
267 * K is the main loop index, decreasing from N to 1 in steps of
268 * 1 or 2
269 *
270  k = n
271  10 CONTINUE
272 *
273 * If K < 1, exit from loop
274 *
275  IF( k.LT.1 )
276  $ GO TO 70
277  kstep = 1
278 *
279 * Determine rows and columns to be interchanged and whether
280 * a 1-by-1 or 2-by-2 pivot block will be used
281 *
282  absakk = cabs1( a( k, k ) )
283 *
284 * IMAX is the row-index of the largest off-diagonal element in
285 * column K, and COLMAX is its absolute value.
286 * Determine both COLMAX and IMAX.
287 *
288  IF( k.GT.1 ) THEN
289  imax = icamax( k-1, a( 1, k ), 1 )
290  colmax = cabs1( a( imax, k ) )
291  ELSE
292  colmax = zero
293  END IF
294 *
295  IF( max( absakk, colmax ).EQ.zero .OR. sisnan(absakk) ) THEN
296 *
297 * Column K is zero or underflow, or contains a NaN:
298 * set INFO and continue
299 *
300  IF( info.EQ.0 )
301  $ info = k
302  kp = k
303  ELSE
304  IF( absakk.GE.alpha*colmax ) THEN
305 *
306 * no interchange, use 1-by-1 pivot block
307 *
308  kp = k
309  ELSE
310 *
311 * JMAX is the column-index of the largest off-diagonal
312 * element in row IMAX, and ROWMAX is its absolute value
313 *
314  jmax = imax + icamax( k-imax, a( imax, imax+1 ), lda )
315  rowmax = cabs1( a( imax, jmax ) )
316  IF( imax.GT.1 ) THEN
317  jmax = icamax( imax-1, a( 1, imax ), 1 )
318  rowmax = max( rowmax, cabs1( a( jmax, imax ) ) )
319  END IF
320 *
321  IF( absakk.GE.alpha*colmax*( colmax / rowmax ) ) THEN
322 *
323 * no interchange, use 1-by-1 pivot block
324 *
325  kp = k
326  ELSE IF( cabs1( a( imax, imax ) ).GE.alpha*rowmax ) THEN
327 *
328 * interchange rows and columns K and IMAX, use 1-by-1
329 * pivot block
330 *
331  kp = imax
332  ELSE
333 *
334 * interchange rows and columns K-1 and IMAX, use 2-by-2
335 * pivot block
336 *
337  kp = imax
338  kstep = 2
339  END IF
340  END IF
341 *
342  kk = k - kstep + 1
343  IF( kp.NE.kk ) THEN
344 *
345 * Interchange rows and columns KK and KP in the leading
346 * submatrix A(1:k,1:k)
347 *
348  CALL cswap( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
349  CALL cswap( kk-kp-1, a( kp+1, kk ), 1, a( kp, kp+1 ),
350  $ lda )
351  t = a( kk, kk )
352  a( kk, kk ) = a( kp, kp )
353  a( kp, kp ) = t
354  IF( kstep.EQ.2 ) THEN
355  t = a( k-1, k )
356  a( k-1, k ) = a( kp, k )
357  a( kp, k ) = t
358  END IF
359  END IF
360 *
361 * Update the leading submatrix
362 *
363  IF( kstep.EQ.1 ) THEN
364 *
365 * 1-by-1 pivot block D(k): column k now holds
366 *
367 * W(k) = U(k)*D(k)
368 *
369 * where U(k) is the k-th column of U
370 *
371 * Perform a rank-1 update of A(1:k-1,1:k-1) as
372 *
373 * A := A - U(k)*D(k)*U(k)**T = A - W(k)*1/D(k)*W(k)**T
374 *
375  r1 = cone / a( k, k )
376  CALL csyr( uplo, k-1, -r1, a( 1, k ), 1, a, lda )
377 *
378 * Store U(k) in column k
379 *
380  CALL cscal( k-1, r1, a( 1, k ), 1 )
381  ELSE
382 *
383 * 2-by-2 pivot block D(k): columns k and k-1 now hold
384 *
385 * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
386 *
387 * where U(k) and U(k-1) are the k-th and (k-1)-th columns
388 * of U
389 *
390 * Perform a rank-2 update of A(1:k-2,1:k-2) as
391 *
392 * A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T
393 * = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )**T
394 *
395  IF( k.GT.2 ) THEN
396 *
397  d12 = a( k-1, k )
398  d22 = a( k-1, k-1 ) / d12
399  d11 = a( k, k ) / d12
400  t = cone / ( d11*d22-cone )
401  d12 = t / d12
402 *
403  DO 30 j = k - 2, 1, -1
404  wkm1 = d12*( d11*a( j, k-1 )-a( j, k ) )
405  wk = d12*( d22*a( j, k )-a( j, k-1 ) )
406  DO 20 i = j, 1, -1
407  a( i, j ) = a( i, j ) - a( i, k )*wk -
408  $ a( i, k-1 )*wkm1
409  20 CONTINUE
410  a( j, k ) = wk
411  a( j, k-1 ) = wkm1
412  30 CONTINUE
413 *
414  END IF
415 *
416  END IF
417  END IF
418 *
419 * Store details of the interchanges in IPIV
420 *
421  IF( kstep.EQ.1 ) THEN
422  ipiv( k ) = kp
423  ELSE
424  ipiv( k ) = -kp
425  ipiv( k-1 ) = -kp
426  END IF
427 *
428 * Decrease K and return to the start of the main loop
429 *
430  k = k - kstep
431  GO TO 10
432 *
433  ELSE
434 *
435 * Factorize A as L*D*L**T using the lower triangle of A
436 *
437 * K is the main loop index, increasing from 1 to N in steps of
438 * 1 or 2
439 *
440  k = 1
441  40 CONTINUE
442 *
443 * If K > N, exit from loop
444 *
445  IF( k.GT.n )
446  $ GO TO 70
447  kstep = 1
448 *
449 * Determine rows and columns to be interchanged and whether
450 * a 1-by-1 or 2-by-2 pivot block will be used
451 *
452  absakk = cabs1( a( k, k ) )
453 *
454 * IMAX is the row-index of the largest off-diagonal element in
455 * column K, and COLMAX is its absolute value.
456 * Determine both COLMAX and IMAX.
457 *
458  IF( k.LT.n ) THEN
459  imax = k + icamax( n-k, a( k+1, k ), 1 )
460  colmax = cabs1( a( imax, k ) )
461  ELSE
462  colmax = zero
463  END IF
464 *
465  IF( max( absakk, colmax ).EQ.zero .OR. sisnan(absakk) ) THEN
466 *
467 * Column K is zero or underflow, or contains a NaN:
468 * set INFO and continue
469 *
470  IF( info.EQ.0 )
471  $ info = k
472  kp = k
473  ELSE
474  IF( absakk.GE.alpha*colmax ) THEN
475 *
476 * no interchange, use 1-by-1 pivot block
477 *
478  kp = k
479  ELSE
480 *
481 * JMAX is the column-index of the largest off-diagonal
482 * element in row IMAX, and ROWMAX is its absolute value
483 *
484  jmax = k - 1 + icamax( imax-k, a( imax, k ), lda )
485  rowmax = cabs1( a( imax, jmax ) )
486  IF( imax.LT.n ) THEN
487  jmax = imax + icamax( n-imax, a( imax+1, imax ), 1 )
488  rowmax = max( rowmax, cabs1( a( jmax, imax ) ) )
489  END IF
490 *
491  IF( absakk.GE.alpha*colmax*( colmax / rowmax ) ) THEN
492 *
493 * no interchange, use 1-by-1 pivot block
494 *
495  kp = k
496  ELSE IF( cabs1( a( imax, imax ) ).GE.alpha*rowmax ) THEN
497 *
498 * interchange rows and columns K and IMAX, use 1-by-1
499 * pivot block
500 *
501  kp = imax
502  ELSE
503 *
504 * interchange rows and columns K+1 and IMAX, use 2-by-2
505 * pivot block
506 *
507  kp = imax
508  kstep = 2
509  END IF
510  END IF
511 *
512  kk = k + kstep - 1
513  IF( kp.NE.kk ) THEN
514 *
515 * Interchange rows and columns KK and KP in the trailing
516 * submatrix A(k:n,k:n)
517 *
518  IF( kp.LT.n )
519  $ CALL cswap( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
520  CALL cswap( kp-kk-1, a( kk+1, kk ), 1, a( kp, kk+1 ),
521  $ lda )
522  t = a( kk, kk )
523  a( kk, kk ) = a( kp, kp )
524  a( kp, kp ) = t
525  IF( kstep.EQ.2 ) THEN
526  t = a( k+1, k )
527  a( k+1, k ) = a( kp, k )
528  a( kp, k ) = t
529  END IF
530  END IF
531 *
532 * Update the trailing submatrix
533 *
534  IF( kstep.EQ.1 ) THEN
535 *
536 * 1-by-1 pivot block D(k): column k now holds
537 *
538 * W(k) = L(k)*D(k)
539 *
540 * where L(k) is the k-th column of L
541 *
542  IF( k.LT.n ) THEN
543 *
544 * Perform a rank-1 update of A(k+1:n,k+1:n) as
545 *
546 * A := A - L(k)*D(k)*L(k)**T = A - W(k)*(1/D(k))*W(k)**T
547 *
548  r1 = cone / a( k, k )
549  CALL csyr( uplo, n-k, -r1, a( k+1, k ), 1,
550  $ a( k+1, k+1 ), lda )
551 *
552 * Store L(k) in column K
553 *
554  CALL cscal( n-k, r1, a( k+1, k ), 1 )
555  END IF
556  ELSE
557 *
558 * 2-by-2 pivot block D(k)
559 *
560  IF( k.LT.n-1 ) THEN
561 *
562 * Perform a rank-2 update of A(k+2:n,k+2:n) as
563 *
564 * A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )**T
565 * = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )**T
566 *
567 * where L(k) and L(k+1) are the k-th and (k+1)-th
568 * columns of L
569 *
570  d21 = a( k+1, k )
571  d11 = a( k+1, k+1 ) / d21
572  d22 = a( k, k ) / d21
573  t = cone / ( d11*d22-cone )
574  d21 = t / d21
575 *
576  DO 60 j = k + 2, n
577  wk = d21*( d11*a( j, k )-a( j, k+1 ) )
578  wkp1 = d21*( d22*a( j, k+1 )-a( j, k ) )
579  DO 50 i = j, n
580  a( i, j ) = a( i, j ) - a( i, k )*wk -
581  $ a( i, k+1 )*wkp1
582  50 CONTINUE
583  a( j, k ) = wk
584  a( j, k+1 ) = wkp1
585  60 CONTINUE
586  END IF
587  END IF
588  END IF
589 *
590 * Store details of the interchanges in IPIV
591 *
592  IF( kstep.EQ.1 ) THEN
593  ipiv( k ) = kp
594  ELSE
595  ipiv( k ) = -kp
596  ipiv( k+1 ) = -kp
597  END IF
598 *
599 * Increase K and return to the start of the main loop
600 *
601  k = k + kstep
602  GO TO 40
603 *
604  END IF
605 *
606  70 CONTINUE
607  RETURN
608 *
609 * End of CSYTF2
610 *
logical function sisnan(SIN)
SISNAN tests input for NaN.
Definition: sisnan.f:61
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 xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:54
integer function icamax(N, CX, INCX)
ICAMAX
Definition: icamax.f:53
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:52
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function: