LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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: