1: SUBROUTINE ZSYTF2( UPLO, N, A, LDA, IPIV, INFO )
2: *
3: * -- LAPACK routine (version 3.2) --
4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6: * November 2006
7: *
8: * .. Scalar Arguments ..
9: CHARACTER UPLO
10: INTEGER INFO, LDA, N
11: * ..
12: * .. Array Arguments ..
13: INTEGER IPIV( * )
14: COMPLEX*16 A( LDA, * )
15: * ..
16: *
17: * Purpose
18: * =======
19: *
20: * ZSYTF2 computes the factorization of a complex symmetric matrix A
21: * using the Bunch-Kaufman diagonal pivoting method:
22: *
23: * A = U*D*U' or A = L*D*L'
24: *
25: * where U (or L) is a product of permutation and unit upper (lower)
26: * triangular matrices, U' is the transpose of U, and D is symmetric and
27: * block diagonal with 1-by-1 and 2-by-2 diagonal blocks.
28: *
29: * This is the unblocked version of the algorithm, calling Level 2 BLAS.
30: *
31: * Arguments
32: * =========
33: *
34: * UPLO (input) CHARACTER*1
35: * Specifies whether the upper or lower triangular part of the
36: * symmetric matrix A is stored:
37: * = 'U': Upper triangular
38: * = 'L': Lower triangular
39: *
40: * N (input) INTEGER
41: * The order of the matrix A. N >= 0.
42: *
43: * A (input/output) COMPLEX*16 array, dimension (LDA,N)
44: * On entry, the symmetric matrix A. If UPLO = 'U', the leading
45: * n-by-n upper triangular part of A contains the upper
46: * triangular part of the matrix A, and the strictly lower
47: * triangular part of A is not referenced. If UPLO = 'L', the
48: * leading n-by-n lower triangular part of A contains the lower
49: * triangular part of the matrix A, and the strictly upper
50: * triangular part of A is not referenced.
51: *
52: * On exit, the block diagonal matrix D and the multipliers used
53: * to obtain the factor U or L (see below for further details).
54: *
55: * LDA (input) INTEGER
56: * The leading dimension of the array A. LDA >= max(1,N).
57: *
58: * IPIV (output) INTEGER array, dimension (N)
59: * Details of the interchanges and the block structure of D.
60: * If IPIV(k) > 0, then rows and columns k and IPIV(k) were
61: * interchanged and D(k,k) is a 1-by-1 diagonal block.
62: * If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
63: * columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
64: * is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) =
65: * IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
66: * interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
67: *
68: * INFO (output) INTEGER
69: * = 0: successful exit
70: * < 0: if INFO = -k, the k-th argument had an illegal value
71: * > 0: if INFO = k, D(k,k) is exactly zero. The factorization
72: * has been completed, but the block diagonal matrix D is
73: * exactly singular, and division by zero will occur if it
74: * is used to solve a system of equations.
75: *
76: * Further Details
77: * ===============
78: *
79: * 09-29-06 - patch from
80: * Bobby Cheng, MathWorks
81: *
82: * Replace l.209 and l.377
83: * IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
84: * by
85: * IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
86: *
87: * 1-96 - Based on modifications by J. Lewis, Boeing Computer Services
88: * Company
89: *
90: * If UPLO = 'U', then A = U*D*U', where
91: * U = P(n)*U(n)* ... *P(k)U(k)* ...,
92: * i.e., U is a product of terms P(k)*U(k), where k decreases from n to
93: * 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
94: * and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
95: * defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
96: * that if the diagonal block D(k) is of order s (s = 1 or 2), then
97: *
98: * ( I v 0 ) k-s
99: * U(k) = ( 0 I 0 ) s
100: * ( 0 0 I ) n-k
101: * k-s s n-k
102: *
103: * If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
104: * If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
105: * and A(k,k), and v overwrites A(1:k-2,k-1:k).
106: *
107: * If UPLO = 'L', then A = L*D*L', where
108: * L = P(1)*L(1)* ... *P(k)*L(k)* ...,
109: * i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
110: * n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
111: * and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
112: * defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
113: * that if the diagonal block D(k) is of order s (s = 1 or 2), then
114: *
115: * ( I 0 0 ) k-1
116: * L(k) = ( 0 I 0 ) s
117: * ( 0 v I ) n-k-s+1
118: * k-1 s n-k-s+1
119: *
120: * If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
121: * If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
122: * and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
123: *
124: * =====================================================================
125: *
126: * .. Parameters ..
127: DOUBLE PRECISION ZERO, ONE
128: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
129: DOUBLE PRECISION EIGHT, SEVTEN
130: PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
131: COMPLEX*16 CONE
132: PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) )
133: * ..
134: * .. Local Scalars ..
135: LOGICAL UPPER
136: INTEGER I, IMAX, J, JMAX, K, KK, KP, KSTEP
137: DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, ROWMAX
138: COMPLEX*16 D11, D12, D21, D22, R1, T, WK, WKM1, WKP1, Z
139: * ..
140: * .. External Functions ..
141: LOGICAL DISNAN, LSAME
142: INTEGER IZAMAX
143: EXTERNAL DISNAN, LSAME, IZAMAX
144: * ..
145: * .. External Subroutines ..
146: EXTERNAL XERBLA, ZSCAL, ZSWAP, ZSYR
147: * ..
148: * .. Intrinsic Functions ..
149: INTRINSIC ABS, DBLE, DIMAG, MAX, SQRT
150: * ..
151: * .. Statement Functions ..
152: DOUBLE PRECISION CABS1
153: * ..
154: * .. Statement Function definitions ..
155: CABS1( Z ) = ABS( DBLE( Z ) ) + ABS( DIMAG( Z ) )
156: * ..
157: * .. Executable Statements ..
158: *
159: * Test the input parameters.
160: *
161: INFO = 0
162: UPPER = LSAME( UPLO, 'U' )
163: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
164: INFO = -1
165: ELSE IF( N.LT.0 ) THEN
166: INFO = -2
167: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
168: INFO = -4
169: END IF
170: IF( INFO.NE.0 ) THEN
171: CALL XERBLA( 'ZSYTF2', -INFO )
172: RETURN
173: END IF
174: *
175: * Initialize ALPHA for use in choosing pivot block size.
176: *
177: ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
178: *
179: IF( UPPER ) THEN
180: *
181: * Factorize A as U*D*U' using the upper triangle of A
182: *
183: * K is the main loop index, decreasing from N to 1 in steps of
184: * 1 or 2
185: *
186: K = N
187: 10 CONTINUE
188: *
189: * If K < 1, exit from loop
190: *
191: IF( K.LT.1 )
192: $ GO TO 70
193: KSTEP = 1
194: *
195: * Determine rows and columns to be interchanged and whether
196: * a 1-by-1 or 2-by-2 pivot block will be used
197: *
198: ABSAKK = CABS1( A( K, K ) )
199: *
200: * IMAX is the row-index of the largest off-diagonal element in
201: * column K, and COLMAX is its absolute value
202: *
203: IF( K.GT.1 ) THEN
204: IMAX = IZAMAX( K-1, A( 1, K ), 1 )
205: COLMAX = CABS1( A( IMAX, K ) )
206: ELSE
207: COLMAX = ZERO
208: END IF
209: *
210: IF( MAX( ABSAKK, COLMAX ).EQ.ZERO .OR. DISNAN(ABSAKK) ) THEN
211: *
212: * Column K is zero or contains a NaN: set INFO and continue
213: *
214: IF( INFO.EQ.0 )
215: $ INFO = K
216: KP = K
217: ELSE
218: IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
219: *
220: * no interchange, use 1-by-1 pivot block
221: *
222: KP = K
223: ELSE
224: *
225: * JMAX is the column-index of the largest off-diagonal
226: * element in row IMAX, and ROWMAX is its absolute value
227: *
228: JMAX = IMAX + IZAMAX( K-IMAX, A( IMAX, IMAX+1 ), LDA )
229: ROWMAX = CABS1( A( IMAX, JMAX ) )
230: IF( IMAX.GT.1 ) THEN
231: JMAX = IZAMAX( IMAX-1, A( 1, IMAX ), 1 )
232: ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) )
233: END IF
234: *
235: IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
236: *
237: * no interchange, use 1-by-1 pivot block
238: *
239: KP = K
240: ELSE IF( CABS1( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX ) THEN
241: *
242: * interchange rows and columns K and IMAX, use 1-by-1
243: * pivot block
244: *
245: KP = IMAX
246: ELSE
247: *
248: * interchange rows and columns K-1 and IMAX, use 2-by-2
249: * pivot block
250: *
251: KP = IMAX
252: KSTEP = 2
253: END IF
254: END IF
255: *
256: KK = K - KSTEP + 1
257: IF( KP.NE.KK ) THEN
258: *
259: * Interchange rows and columns KK and KP in the leading
260: * submatrix A(1:k,1:k)
261: *
262: CALL ZSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
263: CALL ZSWAP( KK-KP-1, A( KP+1, KK ), 1, A( KP, KP+1 ),
264: $ LDA )
265: T = A( KK, KK )
266: A( KK, KK ) = A( KP, KP )
267: A( KP, KP ) = T
268: IF( KSTEP.EQ.2 ) THEN
269: T = A( K-1, K )
270: A( K-1, K ) = A( KP, K )
271: A( KP, K ) = T
272: END IF
273: END IF
274: *
275: * Update the leading submatrix
276: *
277: IF( KSTEP.EQ.1 ) THEN
278: *
279: * 1-by-1 pivot block D(k): column k now holds
280: *
281: * W(k) = U(k)*D(k)
282: *
283: * where U(k) is the k-th column of U
284: *
285: * Perform a rank-1 update of A(1:k-1,1:k-1) as
286: *
287: * A := A - U(k)*D(k)*U(k)' = A - W(k)*1/D(k)*W(k)'
288: *
289: R1 = CONE / A( K, K )
290: CALL ZSYR( UPLO, K-1, -R1, A( 1, K ), 1, A, LDA )
291: *
292: * Store U(k) in column k
293: *
294: CALL ZSCAL( K-1, R1, A( 1, K ), 1 )
295: ELSE
296: *
297: * 2-by-2 pivot block D(k): columns k and k-1 now hold
298: *
299: * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
300: *
301: * where U(k) and U(k-1) are the k-th and (k-1)-th columns
302: * of U
303: *
304: * Perform a rank-2 update of A(1:k-2,1:k-2) as
305: *
306: * A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )'
307: * = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )'
308: *
309: IF( K.GT.2 ) THEN
310: *
311: D12 = A( K-1, K )
312: D22 = A( K-1, K-1 ) / D12
313: D11 = A( K, K ) / D12
314: T = CONE / ( D11*D22-CONE )
315: D12 = T / D12
316: *
317: DO 30 J = K - 2, 1, -1
318: WKM1 = D12*( D11*A( J, K-1 )-A( J, K ) )
319: WK = D12*( D22*A( J, K )-A( J, K-1 ) )
320: DO 20 I = J, 1, -1
321: A( I, J ) = A( I, J ) - A( I, K )*WK -
322: $ A( I, K-1 )*WKM1
323: 20 CONTINUE
324: A( J, K ) = WK
325: A( J, K-1 ) = WKM1
326: 30 CONTINUE
327: *
328: END IF
329: *
330: END IF
331: END IF
332: *
333: * Store details of the interchanges in IPIV
334: *
335: IF( KSTEP.EQ.1 ) THEN
336: IPIV( K ) = KP
337: ELSE
338: IPIV( K ) = -KP
339: IPIV( K-1 ) = -KP
340: END IF
341: *
342: * Decrease K and return to the start of the main loop
343: *
344: K = K - KSTEP
345: GO TO 10
346: *
347: ELSE
348: *
349: * Factorize A as L*D*L' using the lower triangle of A
350: *
351: * K is the main loop index, increasing from 1 to N in steps of
352: * 1 or 2
353: *
354: K = 1
355: 40 CONTINUE
356: *
357: * If K > N, exit from loop
358: *
359: IF( K.GT.N )
360: $ GO TO 70
361: KSTEP = 1
362: *
363: * Determine rows and columns to be interchanged and whether
364: * a 1-by-1 or 2-by-2 pivot block will be used
365: *
366: ABSAKK = CABS1( A( K, K ) )
367: *
368: * IMAX is the row-index of the largest off-diagonal element in
369: * column K, and COLMAX is its absolute value
370: *
371: IF( K.LT.N ) THEN
372: IMAX = K + IZAMAX( N-K, A( K+1, K ), 1 )
373: COLMAX = CABS1( A( IMAX, K ) )
374: ELSE
375: COLMAX = ZERO
376: END IF
377: *
378: IF( MAX( ABSAKK, COLMAX ).EQ.ZERO .OR. DISNAN(ABSAKK) ) THEN
379: *
380: * Column K is zero or contains a NaN: set INFO and continue
381: *
382: IF( INFO.EQ.0 )
383: $ INFO = K
384: KP = K
385: ELSE
386: IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
387: *
388: * no interchange, use 1-by-1 pivot block
389: *
390: KP = K
391: ELSE
392: *
393: * JMAX is the column-index of the largest off-diagonal
394: * element in row IMAX, and ROWMAX is its absolute value
395: *
396: JMAX = K - 1 + IZAMAX( IMAX-K, A( IMAX, K ), LDA )
397: ROWMAX = CABS1( A( IMAX, JMAX ) )
398: IF( IMAX.LT.N ) THEN
399: JMAX = IMAX + IZAMAX( N-IMAX, A( IMAX+1, IMAX ), 1 )
400: ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) )
401: END IF
402: *
403: IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
404: *
405: * no interchange, use 1-by-1 pivot block
406: *
407: KP = K
408: ELSE IF( CABS1( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX ) THEN
409: *
410: * interchange rows and columns K and IMAX, use 1-by-1
411: * pivot block
412: *
413: KP = IMAX
414: ELSE
415: *
416: * interchange rows and columns K+1 and IMAX, use 2-by-2
417: * pivot block
418: *
419: KP = IMAX
420: KSTEP = 2
421: END IF
422: END IF
423: *
424: KK = K + KSTEP - 1
425: IF( KP.NE.KK ) THEN
426: *
427: * Interchange rows and columns KK and KP in the trailing
428: * submatrix A(k:n,k:n)
429: *
430: IF( KP.LT.N )
431: $ CALL ZSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
432: CALL ZSWAP( KP-KK-1, A( KK+1, KK ), 1, A( KP, KK+1 ),
433: $ LDA )
434: T = A( KK, KK )
435: A( KK, KK ) = A( KP, KP )
436: A( KP, KP ) = T
437: IF( KSTEP.EQ.2 ) THEN
438: T = A( K+1, K )
439: A( K+1, K ) = A( KP, K )
440: A( KP, K ) = T
441: END IF
442: END IF
443: *
444: * Update the trailing submatrix
445: *
446: IF( KSTEP.EQ.1 ) THEN
447: *
448: * 1-by-1 pivot block D(k): column k now holds
449: *
450: * W(k) = L(k)*D(k)
451: *
452: * where L(k) is the k-th column of L
453: *
454: IF( K.LT.N ) THEN
455: *
456: * Perform a rank-1 update of A(k+1:n,k+1:n) as
457: *
458: * A := A - L(k)*D(k)*L(k)' = A - W(k)*(1/D(k))*W(k)'
459: *
460: R1 = CONE / A( K, K )
461: CALL ZSYR( UPLO, N-K, -R1, A( K+1, K ), 1,
462: $ A( K+1, K+1 ), LDA )
463: *
464: * Store L(k) in column K
465: *
466: CALL ZSCAL( N-K, R1, A( K+1, K ), 1 )
467: END IF
468: ELSE
469: *
470: * 2-by-2 pivot block D(k)
471: *
472: IF( K.LT.N-1 ) THEN
473: *
474: * Perform a rank-2 update of A(k+2:n,k+2:n) as
475: *
476: * A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )'
477: * = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )'
478: *
479: * where L(k) and L(k+1) are the k-th and (k+1)-th
480: * columns of L
481: *
482: D21 = A( K+1, K )
483: D11 = A( K+1, K+1 ) / D21
484: D22 = A( K, K ) / D21
485: T = CONE / ( D11*D22-CONE )
486: D21 = T / D21
487: *
488: DO 60 J = K + 2, N
489: WK = D21*( D11*A( J, K )-A( J, K+1 ) )
490: WKP1 = D21*( D22*A( J, K+1 )-A( J, K ) )
491: DO 50 I = J, N
492: A( I, J ) = A( I, J ) - A( I, K )*WK -
493: $ A( I, K+1 )*WKP1
494: 50 CONTINUE
495: A( J, K ) = WK
496: A( J, K+1 ) = WKP1
497: 60 CONTINUE
498: END IF
499: END IF
500: END IF
501: *
502: * Store details of the interchanges in IPIV
503: *
504: IF( KSTEP.EQ.1 ) THEN
505: IPIV( K ) = KP
506: ELSE
507: IPIV( K ) = -KP
508: IPIV( K+1 ) = -KP
509: END IF
510: *
511: * Increase K and return to the start of the main loop
512: *
513: K = K + KSTEP
514: GO TO 40
515: *
516: END IF
517: *
518: 70 CONTINUE
519: RETURN
520: *
521: * End of ZSYTF2
522: *
523: 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