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