1: *> \brief \b ZHPTRI
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
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">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZHPTRI( UPLO, N, AP, IPIV, WORK, INFO )
22: *
23: * .. Scalar Arguments ..
24: * CHARACTER UPLO
25: * INTEGER INFO, N
26: * ..
27: * .. Array Arguments ..
28: * INTEGER IPIV( * )
29: * COMPLEX*16 AP( * ), WORK( * )
30: * ..
31: *
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: *
100: *> \author Univ. of Tennessee
101: *> \author Univ. of California Berkeley
102: *> \author Univ. of Colorado Denver
103: *> \author NAG Ltd.
104: *
105: *> \date December 2016
106: *
107: *> \ingroup complex16OTHERcomputational
108: *
109: * =====================================================================
110: SUBROUTINE ZHPTRI( UPLO, N, AP, IPIV, WORK, INFO )
111: *
112: * -- LAPACK computational routine (version 3.7.0) --
113: * -- LAPACK is a software package provided by Univ. of Tennessee, --
114: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
115: * December 2016
116: *
117: * .. Scalar Arguments ..
118: CHARACTER UPLO
119: INTEGER INFO, N
120: * ..
121: * .. Array Arguments ..
122: INTEGER IPIV( * )
123: COMPLEX*16 AP( * ), WORK( * )
124: * ..
125: *
126: * =====================================================================
127: *
128: * .. Parameters ..
129: DOUBLE PRECISION ONE
130: COMPLEX*16 CONE, ZERO
131: PARAMETER ( ONE = 1.0D+0, CONE = ( 1.0D+0, 0.0D+0 ),
132: $ ZERO = ( 0.0D+0, 0.0D+0 ) )
133: * ..
134: * .. Local Scalars ..
135: LOGICAL UPPER
136: INTEGER J, K, KC, KCNEXT, KP, KPC, KSTEP, KX, NPP
137: DOUBLE PRECISION AK, AKP1, D, T
138: COMPLEX*16 AKKP1, TEMP
139: * ..
140: * .. External Functions ..
141: LOGICAL LSAME
142: COMPLEX*16 ZDOTC
143: EXTERNAL LSAME, ZDOTC
144: * ..
145: * .. External Subroutines ..
146: EXTERNAL XERBLA, ZCOPY, ZHPMV, ZSWAP
147: * ..
148: * .. Intrinsic Functions ..
149: INTRINSIC ABS, DBLE, DCONJG
150: * ..
151: * .. Executable Statements ..
152: *
153: * Test the input parameters.
154: *
155: INFO = 0
156: UPPER = LSAME( UPLO, 'U' )
157: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
158: INFO = -1
159: ELSE IF( N.LT.0 ) THEN
160: INFO = -2
161: END IF
162: IF( INFO.NE.0 ) THEN
163: CALL XERBLA( 'ZHPTRI', -INFO )
164: RETURN
165: END IF
166: *
167: * Quick return if possible
168: *
169: IF( N.EQ.0 )
170: $ RETURN
171: *
172: * Check that the diagonal matrix D is nonsingular.
173: *
174: IF( UPPER ) THEN
175: *
176: * Upper triangular storage: examine D from bottom to top
177: *
178: KP = N*( N+1 ) / 2
179: DO 10 INFO = N, 1, -1
180: IF( IPIV( INFO ).GT.0 .AND. AP( KP ).EQ.ZERO )
181: $ RETURN
182: KP = KP - INFO
183: 10 CONTINUE
184: ELSE
185: *
186: * Lower triangular storage: examine D from top to bottom.
187: *
188: KP = 1
189: DO 20 INFO = 1, N
190: IF( IPIV( INFO ).GT.0 .AND. AP( KP ).EQ.ZERO )
191: $ RETURN
192: KP = KP + N - INFO + 1
193: 20 CONTINUE
194: END IF
195: INFO = 0
196: *
197: IF( UPPER ) THEN
198: *
199: * Compute inv(A) from the factorization A = U*D*U**H.
200: *
201: * K is the main loop index, increasing from 1 to N in steps of
202: * 1 or 2, depending on the size of the diagonal blocks.
203: *
204: K = 1
205: KC = 1
206: 30 CONTINUE
207: *
208: * If K > N, exit from loop.
209: *
210: IF( K.GT.N )
211: $ GO TO 50
212: *
213: KCNEXT = KC + K
214: IF( IPIV( K ).GT.0 ) THEN
215: *
216: * 1 x 1 diagonal block
217: *
218: * Invert the diagonal block.
219: *
220: AP( KC+K-1 ) = ONE / DBLE( AP( KC+K-1 ) )
221: *
222: * Compute column K of the inverse.
223: *
224: IF( K.GT.1 ) THEN
225: CALL ZCOPY( K-1, AP( KC ), 1, WORK, 1 )
226: CALL ZHPMV( UPLO, K-1, -CONE, AP, WORK, 1, ZERO,
227: $ AP( KC ), 1 )
228: AP( KC+K-1 ) = AP( KC+K-1 ) -
229: $ DBLE( ZDOTC( K-1, WORK, 1, AP( KC ), 1 ) )
230: END IF
231: KSTEP = 1
232: ELSE
233: *
234: * 2 x 2 diagonal block
235: *
236: * Invert the diagonal block.
237: *
238: T = ABS( AP( KCNEXT+K-1 ) )
239: AK = DBLE( AP( KC+K-1 ) ) / T
240: AKP1 = DBLE( AP( KCNEXT+K ) ) / T
241: AKKP1 = AP( KCNEXT+K-1 ) / T
242: D = T*( AK*AKP1-ONE )
243: AP( KC+K-1 ) = AKP1 / D
244: AP( KCNEXT+K ) = AK / D
245: AP( KCNEXT+K-1 ) = -AKKP1 / D
246: *
247: * Compute columns K and K+1 of the inverse.
248: *
249: IF( K.GT.1 ) THEN
250: CALL ZCOPY( K-1, AP( KC ), 1, WORK, 1 )
251: CALL ZHPMV( UPLO, K-1, -CONE, AP, WORK, 1, ZERO,
252: $ AP( KC ), 1 )
253: AP( KC+K-1 ) = AP( KC+K-1 ) -
254: $ DBLE( ZDOTC( K-1, WORK, 1, AP( KC ), 1 ) )
255: AP( KCNEXT+K-1 ) = AP( KCNEXT+K-1 ) -
256: $ ZDOTC( K-1, AP( KC ), 1, AP( KCNEXT ),
257: $ 1 )
258: CALL ZCOPY( K-1, AP( KCNEXT ), 1, WORK, 1 )
259: CALL ZHPMV( UPLO, K-1, -CONE, AP, WORK, 1, ZERO,
260: $ AP( KCNEXT ), 1 )
261: AP( KCNEXT+K ) = AP( KCNEXT+K ) -
262: $ DBLE( ZDOTC( K-1, WORK, 1, AP( KCNEXT ),
263: $ 1 ) )
264: END IF
265: KSTEP = 2
266: KCNEXT = KCNEXT + K + 1
267: END IF
268: *
269: KP = ABS( IPIV( K ) )
270: IF( KP.NE.K ) THEN
271: *
272: * Interchange rows and columns K and KP in the leading
273: * submatrix A(1:k+1,1:k+1)
274: *
275: KPC = ( KP-1 )*KP / 2 + 1
276: CALL ZSWAP( KP-1, AP( KC ), 1, AP( KPC ), 1 )
277: KX = KPC + KP - 1
278: DO 40 J = KP + 1, K - 1
279: KX = KX + J - 1
280: TEMP = DCONJG( AP( KC+J-1 ) )
281: AP( KC+J-1 ) = DCONJG( AP( KX ) )
282: AP( KX ) = TEMP
283: 40 CONTINUE
284: AP( KC+KP-1 ) = DCONJG( AP( KC+KP-1 ) )
285: TEMP = AP( KC+K-1 )
286: AP( KC+K-1 ) = AP( KPC+KP-1 )
287: AP( KPC+KP-1 ) = TEMP
288: IF( KSTEP.EQ.2 ) THEN
289: TEMP = AP( KC+K+K-1 )
290: AP( KC+K+K-1 ) = AP( KC+K+KP-1 )
291: AP( KC+K+KP-1 ) = TEMP
292: END IF
293: END IF
294: *
295: K = K + KSTEP
296: KC = KCNEXT
297: GO TO 30
298: 50 CONTINUE
299: *
300: ELSE
301: *
302: * Compute inv(A) from the factorization A = L*D*L**H.
303: *
304: * K is the main loop index, increasing from 1 to N in steps of
305: * 1 or 2, depending on the size of the diagonal blocks.
306: *
307: NPP = N*( N+1 ) / 2
308: K = N
309: KC = NPP
310: 60 CONTINUE
311: *
312: * If K < 1, exit from loop.
313: *
314: IF( K.LT.1 )
315: $ GO TO 80
316: *
317: KCNEXT = KC - ( N-K+2 )
318: IF( IPIV( K ).GT.0 ) THEN
319: *
320: * 1 x 1 diagonal block
321: *
322: * Invert the diagonal block.
323: *
324: AP( KC ) = ONE / DBLE( AP( KC ) )
325: *
326: * Compute column K of the inverse.
327: *
328: IF( K.LT.N ) THEN
329: CALL ZCOPY( N-K, AP( KC+1 ), 1, WORK, 1 )
330: CALL ZHPMV( UPLO, N-K, -CONE, AP( KC+N-K+1 ), WORK, 1,
331: $ ZERO, AP( KC+1 ), 1 )
332: AP( KC ) = AP( KC ) - DBLE( ZDOTC( N-K, WORK, 1,
333: $ AP( KC+1 ), 1 ) )
334: END IF
335: KSTEP = 1
336: ELSE
337: *
338: * 2 x 2 diagonal block
339: *
340: * Invert the diagonal block.
341: *
342: T = ABS( AP( KCNEXT+1 ) )
343: AK = DBLE( AP( KCNEXT ) ) / T
344: AKP1 = DBLE( AP( KC ) ) / T
345: AKKP1 = AP( KCNEXT+1 ) / T
346: D = T*( AK*AKP1-ONE )
347: AP( KCNEXT ) = AKP1 / D
348: AP( KC ) = AK / D
349: AP( KCNEXT+1 ) = -AKKP1 / D
350: *
351: * Compute columns K-1 and K of the inverse.
352: *
353: IF( K.LT.N ) THEN
354: CALL ZCOPY( N-K, AP( KC+1 ), 1, WORK, 1 )
355: CALL ZHPMV( UPLO, N-K, -CONE, AP( KC+( N-K+1 ) ), WORK,
356: $ 1, ZERO, AP( KC+1 ), 1 )
357: AP( KC ) = AP( KC ) - DBLE( ZDOTC( N-K, WORK, 1,
358: $ AP( KC+1 ), 1 ) )
359: AP( KCNEXT+1 ) = AP( KCNEXT+1 ) -
360: $ ZDOTC( N-K, AP( KC+1 ), 1,
361: $ AP( KCNEXT+2 ), 1 )
362: CALL ZCOPY( N-K, AP( KCNEXT+2 ), 1, WORK, 1 )
363: CALL ZHPMV( UPLO, N-K, -CONE, AP( KC+( N-K+1 ) ), WORK,
364: $ 1, ZERO, AP( KCNEXT+2 ), 1 )
365: AP( KCNEXT ) = AP( KCNEXT ) -
366: $ DBLE( ZDOTC( N-K, WORK, 1, AP( KCNEXT+2 ),
367: $ 1 ) )
368: END IF
369: KSTEP = 2
370: KCNEXT = KCNEXT - ( N-K+3 )
371: END IF
372: *
373: KP = ABS( IPIV( K ) )
374: IF( KP.NE.K ) THEN
375: *
376: * Interchange rows and columns K and KP in the trailing
377: * submatrix A(k-1:n,k-1:n)
378: *
379: KPC = NPP - ( N-KP+1 )*( N-KP+2 ) / 2 + 1
380: IF( KP.LT.N )
381: $ CALL ZSWAP( N-KP, AP( KC+KP-K+1 ), 1, AP( KPC+1 ), 1 )
382: KX = KC + KP - K
383: DO 70 J = K + 1, KP - 1
384: KX = KX + N - J + 1
385: TEMP = DCONJG( AP( KC+J-K ) )
386: AP( KC+J-K ) = DCONJG( AP( KX ) )
387: AP( KX ) = TEMP
388: 70 CONTINUE
389: AP( KC+KP-K ) = DCONJG( AP( KC+KP-K ) )
390: TEMP = AP( KC )
391: AP( KC ) = AP( KPC )
392: AP( KPC ) = TEMP
393: IF( KSTEP.EQ.2 ) THEN
394: TEMP = AP( KC-N+K-1 )
395: AP( KC-N+K-1 ) = AP( KC-N+KP-1 )
396: AP( KC-N+KP-1 ) = TEMP
397: END IF
398: END IF
399: *
400: K = K - KSTEP
401: KC = KCNEXT
402: GO TO 60
403: 80 CONTINUE
404: END IF
405: *
406: RETURN
407: *
408: * End of ZHPTRI
409: *
410: 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