LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
csytf2.f
Go to the documentation of this file.
1*> \brief \b CSYTF2 computes the factorization of a real symmetric indefinite matrix, using the diagonal pivoting method (unblocked algorithm).
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csytf2.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csytf2.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csytf2.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CSYTF2( UPLO, N, A, LDA, IPIV, INFO )
22*
23* .. Scalar Arguments ..
24* CHARACTER UPLO
25* INTEGER INFO, LDA, N
26* ..
27* .. Array Arguments ..
28* INTEGER IPIV( * )
29* COMPLEX A( LDA, * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> CSYTF2 computes the factorization of a complex symmetric matrix A
39*> using the Bunch-Kaufman diagonal pivoting method:
40*>
41*> A = U*D*U**T or A = L*D*L**T
42*>
43*> where U (or L) is a product of permutation and unit upper (lower)
44*> triangular matrices, U**T is the transpose of U, and D is symmetric and
45*> block diagonal with 1-by-1 and 2-by-2 diagonal blocks.
46*>
47*> This is the unblocked version of the algorithm, calling Level 2 BLAS.
48*> \endverbatim
49*
50* Arguments:
51* ==========
52*
53*> \param[in] UPLO
54*> \verbatim
55*> UPLO is CHARACTER*1
56*> Specifies whether the upper or lower triangular part of the
57*> symmetric matrix A is stored:
58*> = 'U': Upper triangular
59*> = 'L': Lower triangular
60*> \endverbatim
61*>
62*> \param[in] N
63*> \verbatim
64*> N is INTEGER
65*> The order of the matrix A. N >= 0.
66*> \endverbatim
67*>
68*> \param[in,out] A
69*> \verbatim
70*> A is COMPLEX array, dimension (LDA,N)
71*> On entry, the symmetric matrix A. If UPLO = 'U', the leading
72*> n-by-n upper triangular part of A contains the upper
73*> triangular part of the matrix A, and the strictly lower
74*> triangular part of A is not referenced. If UPLO = 'L', the
75*> leading n-by-n lower triangular part of A contains the lower
76*> triangular part of the matrix A, and the strictly upper
77*> triangular part of A is not referenced.
78*>
79*> On exit, the block diagonal matrix D and the multipliers used
80*> to obtain the factor U or L (see below for further details).
81*> \endverbatim
82*>
83*> \param[in] LDA
84*> \verbatim
85*> LDA is INTEGER
86*> The leading dimension of the array A. LDA >= max(1,N).
87*> \endverbatim
88*>
89*> \param[out] IPIV
90*> \verbatim
91*> IPIV is INTEGER array, dimension (N)
92*> Details of the interchanges and the block structure of D.
93*>
94*> If UPLO = 'U':
95*> If IPIV(k) > 0, then rows and columns k and IPIV(k) were
96*> interchanged and D(k,k) is a 1-by-1 diagonal block.
97*>
98*> If IPIV(k) = IPIV(k-1) < 0, then rows and columns
99*> k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
100*> is a 2-by-2 diagonal block.
101*>
102*> If UPLO = 'L':
103*> If IPIV(k) > 0, then rows and columns k and IPIV(k) were
104*> interchanged and D(k,k) is a 1-by-1 diagonal block.
105*>
106*> If IPIV(k) = IPIV(k+1) < 0, then rows and columns
107*> k+1 and -IPIV(k) were interchanged and D(k:k+1,k:k+1)
108*> is a 2-by-2 diagonal block.
109*> \endverbatim
110*>
111*> \param[out] INFO
112*> \verbatim
113*> INFO is INTEGER
114*> = 0: successful exit
115*> < 0: if INFO = -k, the k-th argument had an illegal value
116*> > 0: if INFO = k, D(k,k) is exactly zero. The factorization
117*> has been completed, but the block diagonal matrix D is
118*> exactly singular, and division by zero will occur if it
119*> is used to solve a system of equations.
120*> \endverbatim
121*
122* Authors:
123* ========
124*
125*> \author Univ. of Tennessee
126*> \author Univ. of California Berkeley
127*> \author Univ. of Colorado Denver
128*> \author NAG Ltd.
129*
130*> \ingroup complexSYcomputational
131*
132*> \par Further Details:
133* =====================
134*>
135*> \verbatim
136*>
137*> If UPLO = 'U', then A = U*D*U**T, where
138*> U = P(n)*U(n)* ... *P(k)U(k)* ...,
139*> i.e., U is a product of terms P(k)*U(k), where k decreases from n to
140*> 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
141*> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
142*> defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
143*> that if the diagonal block D(k) is of order s (s = 1 or 2), then
144*>
145*> ( I v 0 ) k-s
146*> U(k) = ( 0 I 0 ) s
147*> ( 0 0 I ) n-k
148*> k-s s n-k
149*>
150*> If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
151*> If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
152*> and A(k,k), and v overwrites A(1:k-2,k-1:k).
153*>
154*> If UPLO = 'L', then A = L*D*L**T, where
155*> L = P(1)*L(1)* ... *P(k)*L(k)* ...,
156*> i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
157*> n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
158*> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
159*> defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
160*> that if the diagonal block D(k) is of order s (s = 1 or 2), then
161*>
162*> ( I 0 0 ) k-1
163*> L(k) = ( 0 I 0 ) s
164*> ( 0 v I ) n-k-s+1
165*> k-1 s n-k-s+1
166*>
167*> If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
168*> If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
169*> and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
170*> \endverbatim
171*
172*> \par Contributors:
173* ==================
174*>
175*> \verbatim
176*>
177*> 09-29-06 - patch from
178*> Bobby Cheng, MathWorks
179*>
180*> Replace l.209 and l.377
181*> IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
182*> by
183*> IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. SISNAN(ABSAKK) ) THEN
184*>
185*> 1-96 - Based on modifications by J. Lewis, Boeing Computer Services
186*> Company
187*> \endverbatim
188*
189* =====================================================================
190 SUBROUTINE csytf2( UPLO, N, A, LDA, IPIV, INFO )
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*
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*
608 END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
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
subroutine csytf2(UPLO, N, A, LDA, IPIV, INFO)
CSYTF2 computes the factorization of a real symmetric indefinite matrix, using the diagonal pivoting ...
Definition: csytf2.f:191