Annotation of rpl/lapack/lapack/zhptri.f, revision 1.18
1.9 bertrand 1: *> \brief \b ZHPTRI
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 ZHPTRI + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhptri.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhptri.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhptri.f">
1.9 bertrand 15: *> [TXT]</a>
1.15 bertrand 16: *> \endhtmlonly
1.9 bertrand 17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZHPTRI( UPLO, N, AP, IPIV, WORK, 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( * ), WORK( * )
30: * ..
1.15 bertrand 31: *
1.9 bertrand 32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *> \verbatim
37: *>
38: *> ZHPTRI computes the inverse of a complex Hermitian indefinite matrix
39: *> A in packed storage using the factorization A = U*D*U**H or
40: *> A = L*D*L**H computed by ZHPTRF.
41: *> \endverbatim
42: *
43: * Arguments:
44: * ==========
45: *
46: *> \param[in] UPLO
47: *> \verbatim
48: *> UPLO is CHARACTER*1
49: *> Specifies whether the details of the factorization are stored
50: *> as an upper or lower triangular matrix.
51: *> = 'U': Upper triangular, form is A = U*D*U**H;
52: *> = 'L': Lower triangular, form is A = L*D*L**H.
53: *> \endverbatim
54: *>
55: *> \param[in] N
56: *> \verbatim
57: *> N is INTEGER
58: *> The order of the matrix A. N >= 0.
59: *> \endverbatim
60: *>
61: *> \param[in,out] AP
62: *> \verbatim
63: *> AP is COMPLEX*16 array, dimension (N*(N+1)/2)
64: *> On entry, the block diagonal matrix D and the multipliers
65: *> used to obtain the factor U or L as computed by ZHPTRF,
66: *> stored as a packed triangular matrix.
67: *>
68: *> On exit, if INFO = 0, the (Hermitian) inverse of the original
69: *> matrix, stored as a packed triangular matrix. The j-th column
70: *> of inv(A) is stored in the array AP as follows:
71: *> if UPLO = 'U', AP(i + (j-1)*j/2) = inv(A)(i,j) for 1<=i<=j;
72: *> if UPLO = 'L',
73: *> AP(i + (j-1)*(2n-j)/2) = inv(A)(i,j) for j<=i<=n.
74: *> \endverbatim
75: *>
76: *> \param[in] IPIV
77: *> \verbatim
78: *> IPIV is INTEGER array, dimension (N)
79: *> Details of the interchanges and the block structure of D
80: *> as determined by ZHPTRF.
81: *> \endverbatim
82: *>
83: *> \param[out] WORK
84: *> \verbatim
85: *> WORK is COMPLEX*16 array, dimension (N)
86: *> \endverbatim
87: *>
88: *> \param[out] INFO
89: *> \verbatim
90: *> INFO is INTEGER
91: *> = 0: successful exit
92: *> < 0: if INFO = -i, the i-th argument had an illegal value
93: *> > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
94: *> inverse could not be computed.
95: *> \endverbatim
96: *
97: * Authors:
98: * ========
99: *
1.15 bertrand 100: *> \author Univ. of Tennessee
101: *> \author Univ. of California Berkeley
102: *> \author Univ. of Colorado Denver
103: *> \author NAG Ltd.
1.9 bertrand 104: *
105: *> \ingroup complex16OTHERcomputational
106: *
107: * =====================================================================
1.1 bertrand 108: SUBROUTINE ZHPTRI( UPLO, N, AP, IPIV, WORK, INFO )
109: *
1.18 ! bertrand 110: * -- LAPACK computational routine --
1.1 bertrand 111: * -- LAPACK is a software package provided by Univ. of Tennessee, --
112: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
113: *
114: * .. Scalar Arguments ..
115: CHARACTER UPLO
116: INTEGER INFO, N
117: * ..
118: * .. Array Arguments ..
119: INTEGER IPIV( * )
120: COMPLEX*16 AP( * ), WORK( * )
121: * ..
122: *
123: * =====================================================================
124: *
125: * .. Parameters ..
126: DOUBLE PRECISION ONE
127: COMPLEX*16 CONE, ZERO
128: PARAMETER ( ONE = 1.0D+0, CONE = ( 1.0D+0, 0.0D+0 ),
129: $ ZERO = ( 0.0D+0, 0.0D+0 ) )
130: * ..
131: * .. Local Scalars ..
132: LOGICAL UPPER
133: INTEGER J, K, KC, KCNEXT, KP, KPC, KSTEP, KX, NPP
134: DOUBLE PRECISION AK, AKP1, D, T
135: COMPLEX*16 AKKP1, TEMP
136: * ..
137: * .. External Functions ..
138: LOGICAL LSAME
139: COMPLEX*16 ZDOTC
140: EXTERNAL LSAME, ZDOTC
141: * ..
142: * .. External Subroutines ..
143: EXTERNAL XERBLA, ZCOPY, ZHPMV, ZSWAP
144: * ..
145: * .. Intrinsic Functions ..
146: INTRINSIC ABS, DBLE, DCONJG
147: * ..
148: * .. Executable Statements ..
149: *
150: * Test the input parameters.
151: *
152: INFO = 0
153: UPPER = LSAME( UPLO, 'U' )
154: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
155: INFO = -1
156: ELSE IF( N.LT.0 ) THEN
157: INFO = -2
158: END IF
159: IF( INFO.NE.0 ) THEN
160: CALL XERBLA( 'ZHPTRI', -INFO )
161: RETURN
162: END IF
163: *
164: * Quick return if possible
165: *
166: IF( N.EQ.0 )
167: $ RETURN
168: *
169: * Check that the diagonal matrix D is nonsingular.
170: *
171: IF( UPPER ) THEN
172: *
173: * Upper triangular storage: examine D from bottom to top
174: *
175: KP = N*( N+1 ) / 2
176: DO 10 INFO = N, 1, -1
177: IF( IPIV( INFO ).GT.0 .AND. AP( KP ).EQ.ZERO )
178: $ RETURN
179: KP = KP - INFO
180: 10 CONTINUE
181: ELSE
182: *
183: * Lower triangular storage: examine D from top to bottom.
184: *
185: KP = 1
186: DO 20 INFO = 1, N
187: IF( IPIV( INFO ).GT.0 .AND. AP( KP ).EQ.ZERO )
188: $ RETURN
189: KP = KP + N - INFO + 1
190: 20 CONTINUE
191: END IF
192: INFO = 0
193: *
194: IF( UPPER ) THEN
195: *
1.8 bertrand 196: * Compute inv(A) from the factorization A = U*D*U**H.
1.1 bertrand 197: *
198: * K is the main loop index, increasing from 1 to N in steps of
199: * 1 or 2, depending on the size of the diagonal blocks.
200: *
201: K = 1
202: KC = 1
203: 30 CONTINUE
204: *
205: * If K > N, exit from loop.
206: *
207: IF( K.GT.N )
208: $ GO TO 50
209: *
210: KCNEXT = KC + K
211: IF( IPIV( K ).GT.0 ) THEN
212: *
213: * 1 x 1 diagonal block
214: *
215: * Invert the diagonal block.
216: *
217: AP( KC+K-1 ) = ONE / DBLE( AP( KC+K-1 ) )
218: *
219: * Compute column K of the inverse.
220: *
221: IF( K.GT.1 ) THEN
222: CALL ZCOPY( K-1, AP( KC ), 1, WORK, 1 )
223: CALL ZHPMV( UPLO, K-1, -CONE, AP, WORK, 1, ZERO,
224: $ AP( KC ), 1 )
225: AP( KC+K-1 ) = AP( KC+K-1 ) -
226: $ DBLE( ZDOTC( K-1, WORK, 1, AP( KC ), 1 ) )
227: END IF
228: KSTEP = 1
229: ELSE
230: *
231: * 2 x 2 diagonal block
232: *
233: * Invert the diagonal block.
234: *
235: T = ABS( AP( KCNEXT+K-1 ) )
236: AK = DBLE( AP( KC+K-1 ) ) / T
237: AKP1 = DBLE( AP( KCNEXT+K ) ) / T
238: AKKP1 = AP( KCNEXT+K-1 ) / T
239: D = T*( AK*AKP1-ONE )
240: AP( KC+K-1 ) = AKP1 / D
241: AP( KCNEXT+K ) = AK / D
242: AP( KCNEXT+K-1 ) = -AKKP1 / D
243: *
244: * Compute columns K and K+1 of the inverse.
245: *
246: IF( K.GT.1 ) THEN
247: CALL ZCOPY( K-1, AP( KC ), 1, WORK, 1 )
248: CALL ZHPMV( UPLO, K-1, -CONE, AP, WORK, 1, ZERO,
249: $ AP( KC ), 1 )
250: AP( KC+K-1 ) = AP( KC+K-1 ) -
251: $ DBLE( ZDOTC( K-1, WORK, 1, AP( KC ), 1 ) )
252: AP( KCNEXT+K-1 ) = AP( KCNEXT+K-1 ) -
253: $ ZDOTC( K-1, AP( KC ), 1, AP( KCNEXT ),
254: $ 1 )
255: CALL ZCOPY( K-1, AP( KCNEXT ), 1, WORK, 1 )
256: CALL ZHPMV( UPLO, K-1, -CONE, AP, WORK, 1, ZERO,
257: $ AP( KCNEXT ), 1 )
258: AP( KCNEXT+K ) = AP( KCNEXT+K ) -
259: $ DBLE( ZDOTC( K-1, WORK, 1, AP( KCNEXT ),
260: $ 1 ) )
261: END IF
262: KSTEP = 2
263: KCNEXT = KCNEXT + K + 1
264: END IF
265: *
266: KP = ABS( IPIV( K ) )
267: IF( KP.NE.K ) THEN
268: *
269: * Interchange rows and columns K and KP in the leading
270: * submatrix A(1:k+1,1:k+1)
271: *
272: KPC = ( KP-1 )*KP / 2 + 1
273: CALL ZSWAP( KP-1, AP( KC ), 1, AP( KPC ), 1 )
274: KX = KPC + KP - 1
275: DO 40 J = KP + 1, K - 1
276: KX = KX + J - 1
277: TEMP = DCONJG( AP( KC+J-1 ) )
278: AP( KC+J-1 ) = DCONJG( AP( KX ) )
279: AP( KX ) = TEMP
280: 40 CONTINUE
281: AP( KC+KP-1 ) = DCONJG( AP( KC+KP-1 ) )
282: TEMP = AP( KC+K-1 )
283: AP( KC+K-1 ) = AP( KPC+KP-1 )
284: AP( KPC+KP-1 ) = TEMP
285: IF( KSTEP.EQ.2 ) THEN
286: TEMP = AP( KC+K+K-1 )
287: AP( KC+K+K-1 ) = AP( KC+K+KP-1 )
288: AP( KC+K+KP-1 ) = TEMP
289: END IF
290: END IF
291: *
292: K = K + KSTEP
293: KC = KCNEXT
294: GO TO 30
295: 50 CONTINUE
296: *
297: ELSE
298: *
1.8 bertrand 299: * Compute inv(A) from the factorization A = L*D*L**H.
1.1 bertrand 300: *
301: * K is the main loop index, increasing from 1 to N in steps of
302: * 1 or 2, depending on the size of the diagonal blocks.
303: *
304: NPP = N*( N+1 ) / 2
305: K = N
306: KC = NPP
307: 60 CONTINUE
308: *
309: * If K < 1, exit from loop.
310: *
311: IF( K.LT.1 )
312: $ GO TO 80
313: *
314: KCNEXT = KC - ( N-K+2 )
315: IF( IPIV( K ).GT.0 ) THEN
316: *
317: * 1 x 1 diagonal block
318: *
319: * Invert the diagonal block.
320: *
321: AP( KC ) = ONE / DBLE( AP( KC ) )
322: *
323: * Compute column K of the inverse.
324: *
325: IF( K.LT.N ) THEN
326: CALL ZCOPY( N-K, AP( KC+1 ), 1, WORK, 1 )
327: CALL ZHPMV( UPLO, N-K, -CONE, AP( KC+N-K+1 ), WORK, 1,
328: $ ZERO, AP( KC+1 ), 1 )
329: AP( KC ) = AP( KC ) - DBLE( ZDOTC( N-K, WORK, 1,
330: $ AP( KC+1 ), 1 ) )
331: END IF
332: KSTEP = 1
333: ELSE
334: *
335: * 2 x 2 diagonal block
336: *
337: * Invert the diagonal block.
338: *
339: T = ABS( AP( KCNEXT+1 ) )
340: AK = DBLE( AP( KCNEXT ) ) / T
341: AKP1 = DBLE( AP( KC ) ) / T
342: AKKP1 = AP( KCNEXT+1 ) / T
343: D = T*( AK*AKP1-ONE )
344: AP( KCNEXT ) = AKP1 / D
345: AP( KC ) = AK / D
346: AP( KCNEXT+1 ) = -AKKP1 / D
347: *
348: * Compute columns K-1 and K of the inverse.
349: *
350: IF( K.LT.N ) THEN
351: CALL ZCOPY( N-K, AP( KC+1 ), 1, WORK, 1 )
352: CALL ZHPMV( UPLO, N-K, -CONE, AP( KC+( N-K+1 ) ), WORK,
353: $ 1, ZERO, AP( KC+1 ), 1 )
354: AP( KC ) = AP( KC ) - DBLE( ZDOTC( N-K, WORK, 1,
355: $ AP( KC+1 ), 1 ) )
356: AP( KCNEXT+1 ) = AP( KCNEXT+1 ) -
357: $ ZDOTC( N-K, AP( KC+1 ), 1,
358: $ AP( KCNEXT+2 ), 1 )
359: CALL ZCOPY( N-K, AP( KCNEXT+2 ), 1, WORK, 1 )
360: CALL ZHPMV( UPLO, N-K, -CONE, AP( KC+( N-K+1 ) ), WORK,
361: $ 1, ZERO, AP( KCNEXT+2 ), 1 )
362: AP( KCNEXT ) = AP( KCNEXT ) -
363: $ DBLE( ZDOTC( N-K, WORK, 1, AP( KCNEXT+2 ),
364: $ 1 ) )
365: END IF
366: KSTEP = 2
367: KCNEXT = KCNEXT - ( N-K+3 )
368: END IF
369: *
370: KP = ABS( IPIV( K ) )
371: IF( KP.NE.K ) THEN
372: *
373: * Interchange rows and columns K and KP in the trailing
374: * submatrix A(k-1:n,k-1:n)
375: *
376: KPC = NPP - ( N-KP+1 )*( N-KP+2 ) / 2 + 1
377: IF( KP.LT.N )
378: $ CALL ZSWAP( N-KP, AP( KC+KP-K+1 ), 1, AP( KPC+1 ), 1 )
379: KX = KC + KP - K
380: DO 70 J = K + 1, KP - 1
381: KX = KX + N - J + 1
382: TEMP = DCONJG( AP( KC+J-K ) )
383: AP( KC+J-K ) = DCONJG( AP( KX ) )
384: AP( KX ) = TEMP
385: 70 CONTINUE
386: AP( KC+KP-K ) = DCONJG( AP( KC+KP-K ) )
387: TEMP = AP( KC )
388: AP( KC ) = AP( KPC )
389: AP( KPC ) = TEMP
390: IF( KSTEP.EQ.2 ) THEN
391: TEMP = AP( KC-N+K-1 )
392: AP( KC-N+K-1 ) = AP( KC-N+KP-1 )
393: AP( KC-N+KP-1 ) = TEMP
394: END IF
395: END IF
396: *
397: K = K - KSTEP
398: KC = KCNEXT
399: GO TO 60
400: 80 CONTINUE
401: END IF
402: *
403: RETURN
404: *
405: * End of ZHPTRI
406: *
407: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zhptri.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack