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