1: SUBROUTINE ZHETF2( 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: * ZHETF2 computes the factorization of a complex Hermitian 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 conjugate transpose of U, and D is
27: * Hermitian and 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: * Hermitian 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 Hermitian 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.210 and l.393
83: * IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
84: * by
85: * IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
86: *
87: * 01-01-96 - Based on modifications by
88: * J. Lewis, Boeing Computer Services Company
89: * A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
90: *
91: * If UPLO = 'U', then A = U*D*U', where
92: * U = P(n)*U(n)* ... *P(k)U(k)* ...,
93: * i.e., U is a product of terms P(k)*U(k), where k decreases from n to
94: * 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
95: * and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
96: * defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
97: * that if the diagonal block D(k) is of order s (s = 1 or 2), then
98: *
99: * ( I v 0 ) k-s
100: * U(k) = ( 0 I 0 ) s
101: * ( 0 0 I ) n-k
102: * k-s s n-k
103: *
104: * If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
105: * If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
106: * and A(k,k), and v overwrites A(1:k-2,k-1:k).
107: *
108: * If UPLO = 'L', then A = L*D*L', where
109: * L = P(1)*L(1)* ... *P(k)*L(k)* ...,
110: * i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
111: * n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
112: * and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
113: * defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
114: * that if the diagonal block D(k) is of order s (s = 1 or 2), then
115: *
116: * ( I 0 0 ) k-1
117: * L(k) = ( 0 I 0 ) s
118: * ( 0 v I ) n-k-s+1
119: * k-1 s n-k-s+1
120: *
121: * If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
122: * If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
123: * and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
124: *
125: * =====================================================================
126: *
127: * .. Parameters ..
128: DOUBLE PRECISION ZERO, ONE
129: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
130: DOUBLE PRECISION EIGHT, SEVTEN
131: PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
132: * ..
133: * .. Local Scalars ..
134: LOGICAL UPPER
135: INTEGER I, IMAX, J, JMAX, K, KK, KP, KSTEP
136: DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D, D11, D22, R1, ROWMAX,
137: $ TT
138: COMPLEX*16 D12, D21, T, WK, WKM1, WKP1, ZDUM
139: * ..
140: * .. External Functions ..
141: LOGICAL LSAME, DISNAN
142: INTEGER IZAMAX
143: DOUBLE PRECISION DLAPY2
144: EXTERNAL LSAME, IZAMAX, DLAPY2, DISNAN
145: * ..
146: * .. External Subroutines ..
147: EXTERNAL XERBLA, ZDSCAL, ZHER, ZSWAP
148: * ..
149: * .. Intrinsic Functions ..
150: INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, SQRT
151: * ..
152: * .. Statement Functions ..
153: DOUBLE PRECISION CABS1
154: * ..
155: * .. Statement Function definitions ..
156: CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
157: * ..
158: * .. Executable Statements ..
159: *
160: * Test the input parameters.
161: *
162: INFO = 0
163: UPPER = LSAME( UPLO, 'U' )
164: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
165: INFO = -1
166: ELSE IF( N.LT.0 ) THEN
167: INFO = -2
168: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
169: INFO = -4
170: END IF
171: IF( INFO.NE.0 ) THEN
172: CALL XERBLA( 'ZHETF2', -INFO )
173: RETURN
174: END IF
175: *
176: * Initialize ALPHA for use in choosing pivot block size.
177: *
178: ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
179: *
180: IF( UPPER ) THEN
181: *
182: * Factorize A as U*D*U' using the upper triangle of A
183: *
184: * K is the main loop index, decreasing from N to 1 in steps of
185: * 1 or 2
186: *
187: K = N
188: 10 CONTINUE
189: *
190: * If K < 1, exit from loop
191: *
192: IF( K.LT.1 )
193: $ GO TO 90
194: KSTEP = 1
195: *
196: * Determine rows and columns to be interchanged and whether
197: * a 1-by-1 or 2-by-2 pivot block will be used
198: *
199: ABSAKK = ABS( DBLE( A( K, K ) ) )
200: *
201: * IMAX is the row-index of the largest off-diagonal element in
202: * column K, and COLMAX is its absolute value
203: *
204: IF( K.GT.1 ) THEN
205: IMAX = IZAMAX( K-1, A( 1, K ), 1 )
206: COLMAX = CABS1( A( IMAX, K ) )
207: ELSE
208: COLMAX = ZERO
209: END IF
210: *
211: IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
212: *
213: * Column K is zero or contains a NaN: set INFO and continue
214: *
215: IF( INFO.EQ.0 )
216: $ INFO = K
217: KP = K
218: A( K, K ) = DBLE( A( K, K ) )
219: ELSE
220: IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
221: *
222: * no interchange, use 1-by-1 pivot block
223: *
224: KP = K
225: ELSE
226: *
227: * JMAX is the column-index of the largest off-diagonal
228: * element in row IMAX, and ROWMAX is its absolute value
229: *
230: JMAX = IMAX + IZAMAX( K-IMAX, A( IMAX, IMAX+1 ), LDA )
231: ROWMAX = CABS1( A( IMAX, JMAX ) )
232: IF( IMAX.GT.1 ) THEN
233: JMAX = IZAMAX( IMAX-1, A( 1, IMAX ), 1 )
234: ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) )
235: END IF
236: *
237: IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
238: *
239: * no interchange, use 1-by-1 pivot block
240: *
241: KP = K
242: ELSE IF( ABS( DBLE( A( IMAX, IMAX ) ) ).GE.ALPHA*ROWMAX )
243: $ THEN
244: *
245: * interchange rows and columns K and IMAX, use 1-by-1
246: * pivot block
247: *
248: KP = IMAX
249: ELSE
250: *
251: * interchange rows and columns K-1 and IMAX, use 2-by-2
252: * pivot block
253: *
254: KP = IMAX
255: KSTEP = 2
256: END IF
257: END IF
258: *
259: KK = K - KSTEP + 1
260: IF( KP.NE.KK ) THEN
261: *
262: * Interchange rows and columns KK and KP in the leading
263: * submatrix A(1:k,1:k)
264: *
265: CALL ZSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
266: DO 20 J = KP + 1, KK - 1
267: T = DCONJG( A( J, KK ) )
268: A( J, KK ) = DCONJG( A( KP, J ) )
269: A( KP, J ) = T
270: 20 CONTINUE
271: A( KP, KK ) = DCONJG( A( KP, KK ) )
272: R1 = DBLE( A( KK, KK ) )
273: A( KK, KK ) = DBLE( A( KP, KP ) )
274: A( KP, KP ) = R1
275: IF( KSTEP.EQ.2 ) THEN
276: A( K, K ) = DBLE( A( K, K ) )
277: T = A( K-1, K )
278: A( K-1, K ) = A( KP, K )
279: A( KP, K ) = T
280: END IF
281: ELSE
282: A( K, K ) = DBLE( A( K, K ) )
283: IF( KSTEP.EQ.2 )
284: $ A( K-1, K-1 ) = DBLE( A( K-1, K-1 ) )
285: END IF
286: *
287: * Update the leading submatrix
288: *
289: IF( KSTEP.EQ.1 ) THEN
290: *
291: * 1-by-1 pivot block D(k): column k now holds
292: *
293: * W(k) = U(k)*D(k)
294: *
295: * where U(k) is the k-th column of U
296: *
297: * Perform a rank-1 update of A(1:k-1,1:k-1) as
298: *
299: * A := A - U(k)*D(k)*U(k)' = A - W(k)*1/D(k)*W(k)'
300: *
301: R1 = ONE / DBLE( A( K, K ) )
302: CALL ZHER( UPLO, K-1, -R1, A( 1, K ), 1, A, LDA )
303: *
304: * Store U(k) in column k
305: *
306: CALL ZDSCAL( K-1, R1, A( 1, K ), 1 )
307: ELSE
308: *
309: * 2-by-2 pivot block D(k): columns k and k-1 now hold
310: *
311: * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
312: *
313: * where U(k) and U(k-1) are the k-th and (k-1)-th columns
314: * of U
315: *
316: * Perform a rank-2 update of A(1:k-2,1:k-2) as
317: *
318: * A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )'
319: * = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )'
320: *
321: IF( K.GT.2 ) THEN
322: *
323: D = DLAPY2( DBLE( A( K-1, K ) ),
324: $ DIMAG( A( K-1, K ) ) )
325: D22 = DBLE( A( K-1, K-1 ) ) / D
326: D11 = DBLE( A( K, K ) ) / D
327: TT = ONE / ( D11*D22-ONE )
328: D12 = A( K-1, K ) / D
329: D = TT / D
330: *
331: DO 40 J = K - 2, 1, -1
332: WKM1 = D*( D11*A( J, K-1 )-DCONJG( D12 )*
333: $ A( J, K ) )
334: WK = D*( D22*A( J, K )-D12*A( J, K-1 ) )
335: DO 30 I = J, 1, -1
336: A( I, J ) = A( I, J ) - A( I, K )*DCONJG( WK ) -
337: $ A( I, K-1 )*DCONJG( WKM1 )
338: 30 CONTINUE
339: A( J, K ) = WK
340: A( J, K-1 ) = WKM1
341: A( J, J ) = DCMPLX( DBLE( A( J, J ) ), 0.0D+0 )
342: 40 CONTINUE
343: *
344: END IF
345: *
346: END IF
347: END IF
348: *
349: * Store details of the interchanges in IPIV
350: *
351: IF( KSTEP.EQ.1 ) THEN
352: IPIV( K ) = KP
353: ELSE
354: IPIV( K ) = -KP
355: IPIV( K-1 ) = -KP
356: END IF
357: *
358: * Decrease K and return to the start of the main loop
359: *
360: K = K - KSTEP
361: GO TO 10
362: *
363: ELSE
364: *
365: * Factorize A as L*D*L' using the lower triangle of A
366: *
367: * K is the main loop index, increasing from 1 to N in steps of
368: * 1 or 2
369: *
370: K = 1
371: 50 CONTINUE
372: *
373: * If K > N, exit from loop
374: *
375: IF( K.GT.N )
376: $ GO TO 90
377: KSTEP = 1
378: *
379: * Determine rows and columns to be interchanged and whether
380: * a 1-by-1 or 2-by-2 pivot block will be used
381: *
382: ABSAKK = ABS( DBLE( A( K, K ) ) )
383: *
384: * IMAX is the row-index of the largest off-diagonal element in
385: * column K, and COLMAX is its absolute value
386: *
387: IF( K.LT.N ) THEN
388: IMAX = K + IZAMAX( N-K, A( K+1, K ), 1 )
389: COLMAX = CABS1( A( IMAX, K ) )
390: ELSE
391: COLMAX = ZERO
392: END IF
393: *
394: IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
395: *
396: * Column K is zero or contains a NaN: set INFO and continue
397: *
398: IF( INFO.EQ.0 )
399: $ INFO = K
400: KP = K
401: A( K, K ) = DBLE( A( K, K ) )
402: ELSE
403: IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
404: *
405: * no interchange, use 1-by-1 pivot block
406: *
407: KP = K
408: ELSE
409: *
410: * JMAX is the column-index of the largest off-diagonal
411: * element in row IMAX, and ROWMAX is its absolute value
412: *
413: JMAX = K - 1 + IZAMAX( IMAX-K, A( IMAX, K ), LDA )
414: ROWMAX = CABS1( A( IMAX, JMAX ) )
415: IF( IMAX.LT.N ) THEN
416: JMAX = IMAX + IZAMAX( N-IMAX, A( IMAX+1, IMAX ), 1 )
417: ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) )
418: END IF
419: *
420: IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
421: *
422: * no interchange, use 1-by-1 pivot block
423: *
424: KP = K
425: ELSE IF( ABS( DBLE( A( IMAX, IMAX ) ) ).GE.ALPHA*ROWMAX )
426: $ THEN
427: *
428: * interchange rows and columns K and IMAX, use 1-by-1
429: * pivot block
430: *
431: KP = IMAX
432: ELSE
433: *
434: * interchange rows and columns K+1 and IMAX, use 2-by-2
435: * pivot block
436: *
437: KP = IMAX
438: KSTEP = 2
439: END IF
440: END IF
441: *
442: KK = K + KSTEP - 1
443: IF( KP.NE.KK ) THEN
444: *
445: * Interchange rows and columns KK and KP in the trailing
446: * submatrix A(k:n,k:n)
447: *
448: IF( KP.LT.N )
449: $ CALL ZSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
450: DO 60 J = KK + 1, KP - 1
451: T = DCONJG( A( J, KK ) )
452: A( J, KK ) = DCONJG( A( KP, J ) )
453: A( KP, J ) = T
454: 60 CONTINUE
455: A( KP, KK ) = DCONJG( A( KP, KK ) )
456: R1 = DBLE( A( KK, KK ) )
457: A( KK, KK ) = DBLE( A( KP, KP ) )
458: A( KP, KP ) = R1
459: IF( KSTEP.EQ.2 ) THEN
460: A( K, K ) = DBLE( A( K, K ) )
461: T = A( K+1, K )
462: A( K+1, K ) = A( KP, K )
463: A( KP, K ) = T
464: END IF
465: ELSE
466: A( K, K ) = DBLE( A( K, K ) )
467: IF( KSTEP.EQ.2 )
468: $ A( K+1, K+1 ) = DBLE( A( K+1, K+1 ) )
469: END IF
470: *
471: * Update the trailing submatrix
472: *
473: IF( KSTEP.EQ.1 ) THEN
474: *
475: * 1-by-1 pivot block D(k): column k now holds
476: *
477: * W(k) = L(k)*D(k)
478: *
479: * where L(k) is the k-th column of L
480: *
481: IF( K.LT.N ) THEN
482: *
483: * Perform a rank-1 update of A(k+1:n,k+1:n) as
484: *
485: * A := A - L(k)*D(k)*L(k)' = A - W(k)*(1/D(k))*W(k)'
486: *
487: R1 = ONE / DBLE( A( K, K ) )
488: CALL ZHER( UPLO, N-K, -R1, A( K+1, K ), 1,
489: $ A( K+1, K+1 ), LDA )
490: *
491: * Store L(k) in column K
492: *
493: CALL ZDSCAL( N-K, R1, A( K+1, K ), 1 )
494: END IF
495: ELSE
496: *
497: * 2-by-2 pivot block D(k)
498: *
499: IF( K.LT.N-1 ) THEN
500: *
501: * Perform a rank-2 update of A(k+2:n,k+2:n) as
502: *
503: * A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )'
504: * = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )'
505: *
506: * where L(k) and L(k+1) are the k-th and (k+1)-th
507: * columns of L
508: *
509: D = DLAPY2( DBLE( A( K+1, K ) ),
510: $ DIMAG( A( K+1, K ) ) )
511: D11 = DBLE( A( K+1, K+1 ) ) / D
512: D22 = DBLE( A( K, K ) ) / D
513: TT = ONE / ( D11*D22-ONE )
514: D21 = A( K+1, K ) / D
515: D = TT / D
516: *
517: DO 80 J = K + 2, N
518: WK = D*( D11*A( J, K )-D21*A( J, K+1 ) )
519: WKP1 = D*( D22*A( J, K+1 )-DCONJG( D21 )*
520: $ A( J, K ) )
521: DO 70 I = J, N
522: A( I, J ) = A( I, J ) - A( I, K )*DCONJG( WK ) -
523: $ A( I, K+1 )*DCONJG( WKP1 )
524: 70 CONTINUE
525: A( J, K ) = WK
526: A( J, K+1 ) = WKP1
527: A( J, J ) = DCMPLX( DBLE( A( J, J ) ), 0.0D+0 )
528: 80 CONTINUE
529: END IF
530: END IF
531: END IF
532: *
533: * Store details of the interchanges in IPIV
534: *
535: IF( KSTEP.EQ.1 ) THEN
536: IPIV( K ) = KP
537: ELSE
538: IPIV( K ) = -KP
539: IPIV( K+1 ) = -KP
540: END IF
541: *
542: * Increase K and return to the start of the main loop
543: *
544: K = K + KSTEP
545: GO TO 50
546: *
547: END IF
548: *
549: 90 CONTINUE
550: RETURN
551: *
552: * End of ZHETF2
553: *
554: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zhetf2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack