LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ csytf2()

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.
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 190 of file csytf2.f.

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