1: *> \brief \b ZSYTF2 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
9: *> Download ZSYTF2 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsytf2.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsytf2.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsytf2.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZSYTF2( 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*16 A( LDA, * )
30: * ..
31: *
32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *> \verbatim
37: *>
38: *> ZSYTF2 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*16 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 complex16SYcomputational
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. DISNAN(ABSAKK) ) THEN
184: *>
185: *> 1-96 - Based on modifications by J. Lewis, Boeing Computer Services
186: *> Company
187: *> \endverbatim
188: *
189: * =====================================================================
190: SUBROUTINE ZSYTF2( 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*16 A( LDA, * )
203: * ..
204: *
205: * =====================================================================
206: *
207: * .. Parameters ..
208: DOUBLE PRECISION ZERO, ONE
209: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
210: DOUBLE PRECISION EIGHT, SEVTEN
211: PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
212: COMPLEX*16 CONE
213: PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) )
214: * ..
215: * .. Local Scalars ..
216: LOGICAL UPPER
217: INTEGER I, IMAX, J, JMAX, K, KK, KP, KSTEP
218: DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, ROWMAX
219: COMPLEX*16 D11, D12, D21, D22, R1, T, WK, WKM1, WKP1, Z
220: * ..
221: * .. External Functions ..
222: LOGICAL DISNAN, LSAME
223: INTEGER IZAMAX
224: EXTERNAL DISNAN, LSAME, IZAMAX
225: * ..
226: * .. External Subroutines ..
227: EXTERNAL XERBLA, ZSCAL, ZSWAP, ZSYR
228: * ..
229: * .. Intrinsic Functions ..
230: INTRINSIC ABS, DBLE, DIMAG, MAX, SQRT
231: * ..
232: * .. Statement Functions ..
233: DOUBLE PRECISION CABS1
234: * ..
235: * .. Statement Function definitions ..
236: CABS1( Z ) = ABS( DBLE( Z ) ) + ABS( DIMAG( 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( 'ZSYTF2', -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 = IZAMAX( 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. DISNAN(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 + IZAMAX( K-IMAX, A( IMAX, IMAX+1 ), LDA )
312: ROWMAX = CABS1( A( IMAX, JMAX ) )
313: IF( IMAX.GT.1 ) THEN
314: JMAX = IZAMAX( 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 ZSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
346: CALL ZSWAP( 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 ZSYR( UPLO, K-1, -R1, A( 1, K ), 1, A, LDA )
374: *
375: * Store U(k) in column k
376: *
377: CALL ZSCAL( 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 + IZAMAX( 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. DISNAN(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 + IZAMAX( IMAX-K, A( IMAX, K ), LDA )
482: ROWMAX = CABS1( A( IMAX, JMAX ) )
483: IF( IMAX.LT.N ) THEN
484: JMAX = IMAX + IZAMAX( 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 ZSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
517: CALL ZSWAP( 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 ZSYR( 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 ZSCAL( 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 ZSYTF2
607: *
608: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zsytf2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack