Annotation of rpl/lapack/lapack/zhetf2.f, revision 1.15
1.14 bertrand 1: *> \brief \b ZHETF2 computes the factorization of a complex Hermitian matrix, using the diagonal pivoting method (unblocked algorithm, calling Level 2 BLAS).
1.9 bertrand 2: *
3: * =========== DOCUMENTATION ===========
4: *
1.14 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.9 bertrand 7: *
8: *> \htmlonly
1.14 bertrand 9: *> Download ZHETF2 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetf2.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetf2.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetf2.f">
1.9 bertrand 15: *> [TXT]</a>
1.14 bertrand 16: *> \endhtmlonly
1.9 bertrand 17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZHETF2( UPLO, N, A, LDA, IPIV, INFO )
1.14 bertrand 22: *
1.9 bertrand 23: * .. Scalar Arguments ..
24: * CHARACTER UPLO
25: * INTEGER INFO, LDA, N
26: * ..
27: * .. Array Arguments ..
28: * INTEGER IPIV( * )
29: * COMPLEX*16 A( LDA, * )
30: * ..
1.14 bertrand 31: *
1.9 bertrand 32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *> \verbatim
37: *>
38: *> ZHETF2 computes the factorization of a complex Hermitian matrix A
39: *> using the Bunch-Kaufman diagonal pivoting method:
40: *>
41: *> A = U*D*U**H or A = L*D*L**H
42: *>
43: *> where U (or L) is a product of permutation and unit upper (lower)
44: *> triangular matrices, U**H is the conjugate transpose of U, and D is
45: *> Hermitian and 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: *> Hermitian 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 Hermitian 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.
1.14 bertrand 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.
1.9 bertrand 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: *
1.14 bertrand 125: *> \author Univ. of Tennessee
126: *> \author Univ. of California Berkeley
127: *> \author Univ. of Colorado Denver
128: *> \author NAG Ltd.
1.9 bertrand 129: *
1.14 bertrand 130: *> \date November 2013
1.9 bertrand 131: *
132: *> \ingroup complex16HEcomputational
133: *
134: *> \par Further Details:
135: * =====================
136: *>
137: *> \verbatim
138: *>
139: *> If UPLO = 'U', then A = U*D*U**H, where
140: *> U = P(n)*U(n)* ... *P(k)U(k)* ...,
141: *> i.e., U is a product of terms P(k)*U(k), where k decreases from n to
142: *> 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
143: *> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
144: *> defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
145: *> that if the diagonal block D(k) is of order s (s = 1 or 2), then
146: *>
147: *> ( I v 0 ) k-s
148: *> U(k) = ( 0 I 0 ) s
149: *> ( 0 0 I ) n-k
150: *> k-s s n-k
151: *>
152: *> If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
153: *> If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
154: *> and A(k,k), and v overwrites A(1:k-2,k-1:k).
155: *>
156: *> If UPLO = 'L', then A = L*D*L**H, where
157: *> L = P(1)*L(1)* ... *P(k)*L(k)* ...,
158: *> i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
159: *> n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
160: *> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
161: *> defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
162: *> that if the diagonal block D(k) is of order s (s = 1 or 2), then
163: *>
164: *> ( I 0 0 ) k-1
165: *> L(k) = ( 0 I 0 ) s
166: *> ( 0 v I ) n-k-s+1
167: *> k-1 s n-k-s+1
168: *>
169: *> If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
170: *> If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
171: *> and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
172: *> \endverbatim
173: *
174: *> \par Contributors:
175: * ==================
176: *>
177: *> \verbatim
178: *> 09-29-06 - patch from
179: *> Bobby Cheng, MathWorks
180: *>
181: *> Replace l.210 and l.393
182: *> IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
183: *> by
184: *> IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
185: *>
186: *> 01-01-96 - Based on modifications by
187: *> J. Lewis, Boeing Computer Services Company
188: *> A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
189: *> \endverbatim
190: *
191: * =====================================================================
1.1 bertrand 192: SUBROUTINE ZHETF2( UPLO, N, A, LDA, IPIV, INFO )
193: *
1.14 bertrand 194: * -- LAPACK computational routine (version 3.5.0) --
1.1 bertrand 195: * -- LAPACK is a software package provided by Univ. of Tennessee, --
196: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.14 bertrand 197: * November 2013
1.1 bertrand 198: *
199: * .. Scalar Arguments ..
200: CHARACTER UPLO
201: INTEGER INFO, LDA, N
202: * ..
203: * .. Array Arguments ..
204: INTEGER IPIV( * )
205: COMPLEX*16 A( LDA, * )
206: * ..
207: *
208: * =====================================================================
209: *
210: * .. Parameters ..
211: DOUBLE PRECISION ZERO, ONE
212: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
213: DOUBLE PRECISION EIGHT, SEVTEN
214: PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
215: * ..
216: * .. Local Scalars ..
217: LOGICAL UPPER
218: INTEGER I, IMAX, J, JMAX, K, KK, KP, KSTEP
219: DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D, D11, D22, R1, ROWMAX,
220: $ TT
221: COMPLEX*16 D12, D21, T, WK, WKM1, WKP1, ZDUM
222: * ..
223: * .. External Functions ..
224: LOGICAL LSAME, DISNAN
225: INTEGER IZAMAX
226: DOUBLE PRECISION DLAPY2
227: EXTERNAL LSAME, IZAMAX, DLAPY2, DISNAN
228: * ..
229: * .. External Subroutines ..
230: EXTERNAL XERBLA, ZDSCAL, ZHER, ZSWAP
231: * ..
232: * .. Intrinsic Functions ..
233: INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, SQRT
234: * ..
235: * .. Statement Functions ..
236: DOUBLE PRECISION CABS1
237: * ..
238: * .. Statement Function definitions ..
239: CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
240: * ..
241: * .. Executable Statements ..
242: *
243: * Test the input parameters.
244: *
245: INFO = 0
246: UPPER = LSAME( UPLO, 'U' )
247: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
248: INFO = -1
249: ELSE IF( N.LT.0 ) THEN
250: INFO = -2
251: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
252: INFO = -4
253: END IF
254: IF( INFO.NE.0 ) THEN
255: CALL XERBLA( 'ZHETF2', -INFO )
256: RETURN
257: END IF
258: *
259: * Initialize ALPHA for use in choosing pivot block size.
260: *
261: ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
262: *
263: IF( UPPER ) THEN
264: *
1.8 bertrand 265: * Factorize A as U*D*U**H using the upper triangle of A
1.1 bertrand 266: *
267: * K is the main loop index, decreasing from N to 1 in steps of
268: * 1 or 2
269: *
270: K = N
271: 10 CONTINUE
272: *
273: * If K < 1, exit from loop
274: *
275: IF( K.LT.1 )
276: $ GO TO 90
277: KSTEP = 1
278: *
279: * Determine rows and columns to be interchanged and whether
280: * a 1-by-1 or 2-by-2 pivot block will be used
281: *
282: ABSAKK = ABS( DBLE( A( K, K ) ) )
283: *
284: * IMAX is the row-index of the largest off-diagonal element in
1.14 bertrand 285: * column K, and COLMAX is its absolute value.
286: * Determine both COLMAX and IMAX.
1.1 bertrand 287: *
288: IF( K.GT.1 ) THEN
289: IMAX = IZAMAX( K-1, A( 1, K ), 1 )
290: COLMAX = CABS1( A( IMAX, K ) )
291: ELSE
292: COLMAX = ZERO
293: END IF
294: *
295: IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
296: *
1.14 bertrand 297: * Column K is zero or underflow, or contains a NaN:
298: * set INFO and continue
1.1 bertrand 299: *
300: IF( INFO.EQ.0 )
301: $ INFO = K
302: KP = K
303: A( K, K ) = DBLE( A( K, K ) )
304: ELSE
1.14 bertrand 305: *
306: * ============================================================
307: *
308: * Test for interchange
309: *
1.1 bertrand 310: IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
311: *
312: * no interchange, use 1-by-1 pivot block
313: *
314: KP = K
315: ELSE
316: *
317: * JMAX is the column-index of the largest off-diagonal
1.14 bertrand 318: * element in row IMAX, and ROWMAX is its absolute value.
319: * Determine only ROWMAX.
1.1 bertrand 320: *
321: JMAX = IMAX + IZAMAX( K-IMAX, A( IMAX, IMAX+1 ), LDA )
322: ROWMAX = CABS1( A( IMAX, JMAX ) )
323: IF( IMAX.GT.1 ) THEN
324: JMAX = IZAMAX( IMAX-1, A( 1, IMAX ), 1 )
325: ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) )
326: END IF
327: *
328: IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
329: *
330: * no interchange, use 1-by-1 pivot block
331: *
332: KP = K
1.14 bertrand 333: *
1.1 bertrand 334: ELSE IF( ABS( DBLE( A( IMAX, IMAX ) ) ).GE.ALPHA*ROWMAX )
335: $ THEN
336: *
337: * interchange rows and columns K and IMAX, use 1-by-1
338: * pivot block
339: *
340: KP = IMAX
341: ELSE
342: *
343: * interchange rows and columns K-1 and IMAX, use 2-by-2
344: * pivot block
345: *
346: KP = IMAX
347: KSTEP = 2
348: END IF
1.14 bertrand 349: *
1.1 bertrand 350: END IF
351: *
1.14 bertrand 352: * ============================================================
353: *
1.1 bertrand 354: KK = K - KSTEP + 1
355: IF( KP.NE.KK ) THEN
356: *
357: * Interchange rows and columns KK and KP in the leading
358: * submatrix A(1:k,1:k)
359: *
360: CALL ZSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
361: DO 20 J = KP + 1, KK - 1
362: T = DCONJG( A( J, KK ) )
363: A( J, KK ) = DCONJG( A( KP, J ) )
364: A( KP, J ) = T
365: 20 CONTINUE
366: A( KP, KK ) = DCONJG( A( KP, KK ) )
367: R1 = DBLE( A( KK, KK ) )
368: A( KK, KK ) = DBLE( A( KP, KP ) )
369: A( KP, KP ) = R1
370: IF( KSTEP.EQ.2 ) THEN
371: A( K, K ) = DBLE( A( K, K ) )
372: T = A( K-1, K )
373: A( K-1, K ) = A( KP, K )
374: A( KP, K ) = T
375: END IF
376: ELSE
377: A( K, K ) = DBLE( A( K, K ) )
378: IF( KSTEP.EQ.2 )
379: $ A( K-1, K-1 ) = DBLE( A( K-1, K-1 ) )
380: END IF
381: *
382: * Update the leading submatrix
383: *
384: IF( KSTEP.EQ.1 ) THEN
385: *
386: * 1-by-1 pivot block D(k): column k now holds
387: *
388: * W(k) = U(k)*D(k)
389: *
390: * where U(k) is the k-th column of U
391: *
392: * Perform a rank-1 update of A(1:k-1,1:k-1) as
393: *
1.8 bertrand 394: * A := A - U(k)*D(k)*U(k)**H = A - W(k)*1/D(k)*W(k)**H
1.1 bertrand 395: *
396: R1 = ONE / DBLE( A( K, K ) )
397: CALL ZHER( UPLO, K-1, -R1, A( 1, K ), 1, A, LDA )
398: *
399: * Store U(k) in column k
400: *
401: CALL ZDSCAL( K-1, R1, A( 1, K ), 1 )
402: ELSE
403: *
404: * 2-by-2 pivot block D(k): columns k and k-1 now hold
405: *
406: * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
407: *
408: * where U(k) and U(k-1) are the k-th and (k-1)-th columns
409: * of U
410: *
411: * Perform a rank-2 update of A(1:k-2,1:k-2) as
412: *
1.8 bertrand 413: * A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**H
414: * = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )**H
1.1 bertrand 415: *
416: IF( K.GT.2 ) THEN
417: *
418: D = DLAPY2( DBLE( A( K-1, K ) ),
419: $ DIMAG( A( K-1, K ) ) )
420: D22 = DBLE( A( K-1, K-1 ) ) / D
421: D11 = DBLE( A( K, K ) ) / D
422: TT = ONE / ( D11*D22-ONE )
423: D12 = A( K-1, K ) / D
424: D = TT / D
425: *
426: DO 40 J = K - 2, 1, -1
427: WKM1 = D*( D11*A( J, K-1 )-DCONJG( D12 )*
428: $ A( J, K ) )
429: WK = D*( D22*A( J, K )-D12*A( J, K-1 ) )
430: DO 30 I = J, 1, -1
431: A( I, J ) = A( I, J ) - A( I, K )*DCONJG( WK ) -
432: $ A( I, K-1 )*DCONJG( WKM1 )
433: 30 CONTINUE
434: A( J, K ) = WK
435: A( J, K-1 ) = WKM1
436: A( J, J ) = DCMPLX( DBLE( A( J, J ) ), 0.0D+0 )
437: 40 CONTINUE
438: *
439: END IF
440: *
441: END IF
442: END IF
443: *
444: * Store details of the interchanges in IPIV
445: *
446: IF( KSTEP.EQ.1 ) THEN
447: IPIV( K ) = KP
448: ELSE
449: IPIV( K ) = -KP
450: IPIV( K-1 ) = -KP
451: END IF
452: *
453: * Decrease K and return to the start of the main loop
454: *
455: K = K - KSTEP
456: GO TO 10
457: *
458: ELSE
459: *
1.8 bertrand 460: * Factorize A as L*D*L**H using the lower triangle of A
1.1 bertrand 461: *
462: * K is the main loop index, increasing from 1 to N in steps of
463: * 1 or 2
464: *
465: K = 1
466: 50 CONTINUE
467: *
468: * If K > N, exit from loop
469: *
470: IF( K.GT.N )
471: $ GO TO 90
472: KSTEP = 1
473: *
474: * Determine rows and columns to be interchanged and whether
475: * a 1-by-1 or 2-by-2 pivot block will be used
476: *
477: ABSAKK = ABS( DBLE( A( K, K ) ) )
478: *
479: * IMAX is the row-index of the largest off-diagonal element in
1.14 bertrand 480: * column K, and COLMAX is its absolute value.
481: * Determine both COLMAX and IMAX.
1.1 bertrand 482: *
483: IF( K.LT.N ) THEN
484: IMAX = K + IZAMAX( N-K, A( K+1, K ), 1 )
485: COLMAX = CABS1( A( IMAX, K ) )
486: ELSE
487: COLMAX = ZERO
488: END IF
489: *
490: IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
491: *
1.14 bertrand 492: * Column K is zero or underflow, or contains a NaN:
493: * set INFO and continue
1.1 bertrand 494: *
495: IF( INFO.EQ.0 )
496: $ INFO = K
497: KP = K
498: A( K, K ) = DBLE( A( K, K ) )
499: ELSE
1.14 bertrand 500: *
501: * ============================================================
502: *
503: * Test for interchange
504: *
1.1 bertrand 505: IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
506: *
507: * no interchange, use 1-by-1 pivot block
508: *
509: KP = K
510: ELSE
511: *
512: * JMAX is the column-index of the largest off-diagonal
1.14 bertrand 513: * element in row IMAX, and ROWMAX is its absolute value.
514: * Determine only ROWMAX.
1.1 bertrand 515: *
516: JMAX = K - 1 + IZAMAX( IMAX-K, A( IMAX, K ), LDA )
517: ROWMAX = CABS1( A( IMAX, JMAX ) )
518: IF( IMAX.LT.N ) THEN
519: JMAX = IMAX + IZAMAX( N-IMAX, A( IMAX+1, IMAX ), 1 )
520: ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) )
521: END IF
522: *
523: IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
524: *
525: * no interchange, use 1-by-1 pivot block
526: *
527: KP = K
1.14 bertrand 528: *
1.1 bertrand 529: ELSE IF( ABS( DBLE( A( IMAX, IMAX ) ) ).GE.ALPHA*ROWMAX )
530: $ THEN
531: *
532: * interchange rows and columns K and IMAX, use 1-by-1
533: * pivot block
534: *
535: KP = IMAX
536: ELSE
537: *
538: * interchange rows and columns K+1 and IMAX, use 2-by-2
539: * pivot block
540: *
541: KP = IMAX
542: KSTEP = 2
543: END IF
1.14 bertrand 544: *
1.1 bertrand 545: END IF
546: *
1.14 bertrand 547: * ============================================================
548: *
1.1 bertrand 549: KK = K + KSTEP - 1
550: IF( KP.NE.KK ) THEN
551: *
552: * Interchange rows and columns KK and KP in the trailing
553: * submatrix A(k:n,k:n)
554: *
555: IF( KP.LT.N )
556: $ CALL ZSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
557: DO 60 J = KK + 1, KP - 1
558: T = DCONJG( A( J, KK ) )
559: A( J, KK ) = DCONJG( A( KP, J ) )
560: A( KP, J ) = T
561: 60 CONTINUE
562: A( KP, KK ) = DCONJG( A( KP, KK ) )
563: R1 = DBLE( A( KK, KK ) )
564: A( KK, KK ) = DBLE( A( KP, KP ) )
565: A( KP, KP ) = R1
566: IF( KSTEP.EQ.2 ) THEN
567: A( K, K ) = DBLE( A( K, K ) )
568: T = A( K+1, K )
569: A( K+1, K ) = A( KP, K )
570: A( KP, K ) = T
571: END IF
572: ELSE
573: A( K, K ) = DBLE( A( K, K ) )
574: IF( KSTEP.EQ.2 )
575: $ A( K+1, K+1 ) = DBLE( A( K+1, K+1 ) )
576: END IF
577: *
578: * Update the trailing submatrix
579: *
580: IF( KSTEP.EQ.1 ) THEN
581: *
582: * 1-by-1 pivot block D(k): column k now holds
583: *
584: * W(k) = L(k)*D(k)
585: *
586: * where L(k) is the k-th column of L
587: *
588: IF( K.LT.N ) THEN
589: *
590: * Perform a rank-1 update of A(k+1:n,k+1:n) as
591: *
1.8 bertrand 592: * A := A - L(k)*D(k)*L(k)**H = A - W(k)*(1/D(k))*W(k)**H
1.1 bertrand 593: *
594: R1 = ONE / DBLE( A( K, K ) )
595: CALL ZHER( UPLO, N-K, -R1, A( K+1, K ), 1,
596: $ A( K+1, K+1 ), LDA )
597: *
598: * Store L(k) in column K
599: *
600: CALL ZDSCAL( N-K, R1, A( K+1, K ), 1 )
601: END IF
602: ELSE
603: *
604: * 2-by-2 pivot block D(k)
605: *
606: IF( K.LT.N-1 ) THEN
607: *
608: * Perform a rank-2 update of A(k+2:n,k+2:n) as
609: *
1.8 bertrand 610: * A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )**H
611: * = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )**H
1.1 bertrand 612: *
613: * where L(k) and L(k+1) are the k-th and (k+1)-th
614: * columns of L
615: *
616: D = DLAPY2( DBLE( A( K+1, K ) ),
617: $ DIMAG( A( K+1, K ) ) )
618: D11 = DBLE( A( K+1, K+1 ) ) / D
619: D22 = DBLE( A( K, K ) ) / D
620: TT = ONE / ( D11*D22-ONE )
621: D21 = A( K+1, K ) / D
622: D = TT / D
623: *
624: DO 80 J = K + 2, N
625: WK = D*( D11*A( J, K )-D21*A( J, K+1 ) )
626: WKP1 = D*( D22*A( J, K+1 )-DCONJG( D21 )*
627: $ A( J, K ) )
628: DO 70 I = J, N
629: A( I, J ) = A( I, J ) - A( I, K )*DCONJG( WK ) -
630: $ A( I, K+1 )*DCONJG( WKP1 )
631: 70 CONTINUE
632: A( J, K ) = WK
633: A( J, K+1 ) = WKP1
634: A( J, J ) = DCMPLX( DBLE( A( J, J ) ), 0.0D+0 )
635: 80 CONTINUE
636: END IF
637: END IF
638: END IF
639: *
640: * Store details of the interchanges in IPIV
641: *
642: IF( KSTEP.EQ.1 ) THEN
643: IPIV( K ) = KP
644: ELSE
645: IPIV( K ) = -KP
646: IPIV( K+1 ) = -KP
647: END IF
648: *
649: * Increase K and return to the start of the main loop
650: *
651: K = K + KSTEP
652: GO TO 50
653: *
654: END IF
655: *
656: 90 CONTINUE
657: RETURN
658: *
659: * End of ZHETF2
660: *
661: 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