![[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|
Return to zhptri.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: SUBROUTINE ZHPTRI( UPLO, N, AP, IPIV, WORK, INFO )
2: *
3: * -- LAPACK routine (version 3.2) --
4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6: * November 2006
7: *
8: * .. Scalar Arguments ..
9: CHARACTER UPLO
10: INTEGER INFO, N
11: * ..
12: * .. Array Arguments ..
13: INTEGER IPIV( * )
14: COMPLEX*16 AP( * ), WORK( * )
15: * ..
16: *
17: * Purpose
18: * =======
19: *
20: * ZHPTRI computes the inverse of a complex Hermitian indefinite matrix
21: * A in packed storage using the factorization A = U*D*U**H or
22: * A = L*D*L**H computed by ZHPTRF.
23: *
24: * Arguments
25: * =========
26: *
27: * UPLO (input) CHARACTER*1
28: * Specifies whether the details of the factorization are stored
29: * as an upper or lower triangular matrix.
30: * = 'U': Upper triangular, form is A = U*D*U**H;
31: * = 'L': Lower triangular, form is A = L*D*L**H.
32: *
33: * N (input) INTEGER
34: * The order of the matrix A. N >= 0.
35: *
36: * AP (input/output) COMPLEX*16 array, dimension (N*(N+1)/2)
37: * On entry, the block diagonal matrix D and the multipliers
38: * used to obtain the factor U or L as computed by ZHPTRF,
39: * stored as a packed triangular matrix.
40: *
41: * On exit, if INFO = 0, the (Hermitian) inverse of the original
42: * matrix, stored as a packed triangular matrix. The j-th column
43: * of inv(A) is stored in the array AP as follows:
44: * if UPLO = 'U', AP(i + (j-1)*j/2) = inv(A)(i,j) for 1<=i<=j;
45: * if UPLO = 'L',
46: * AP(i + (j-1)*(2n-j)/2) = inv(A)(i,j) for j<=i<=n.
47: *
48: * IPIV (input) INTEGER array, dimension (N)
49: * Details of the interchanges and the block structure of D
50: * as determined by ZHPTRF.
51: *
52: * WORK (workspace) COMPLEX*16 array, dimension (N)
53: *
54: * INFO (output) INTEGER
55: * = 0: successful exit
56: * < 0: if INFO = -i, the i-th argument had an illegal value
57: * > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
58: * inverse could not be computed.
59: *
60: * =====================================================================
61: *
62: * .. Parameters ..
63: DOUBLE PRECISION ONE
64: COMPLEX*16 CONE, ZERO
65: PARAMETER ( ONE = 1.0D+0, CONE = ( 1.0D+0, 0.0D+0 ),
66: $ ZERO = ( 0.0D+0, 0.0D+0 ) )
67: * ..
68: * .. Local Scalars ..
69: LOGICAL UPPER
70: INTEGER J, K, KC, KCNEXT, KP, KPC, KSTEP, KX, NPP
71: DOUBLE PRECISION AK, AKP1, D, T
72: COMPLEX*16 AKKP1, TEMP
73: * ..
74: * .. External Functions ..
75: LOGICAL LSAME
76: COMPLEX*16 ZDOTC
77: EXTERNAL LSAME, ZDOTC
78: * ..
79: * .. External Subroutines ..
80: EXTERNAL XERBLA, ZCOPY, ZHPMV, ZSWAP
81: * ..
82: * .. Intrinsic Functions ..
83: INTRINSIC ABS, DBLE, DCONJG
84: * ..
85: * .. Executable Statements ..
86: *
87: * Test the input parameters.
88: *
89: INFO = 0
90: UPPER = LSAME( UPLO, 'U' )
91: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
92: INFO = -1
93: ELSE IF( N.LT.0 ) THEN
94: INFO = -2
95: END IF
96: IF( INFO.NE.0 ) THEN
97: CALL XERBLA( 'ZHPTRI', -INFO )
98: RETURN
99: END IF
100: *
101: * Quick return if possible
102: *
103: IF( N.EQ.0 )
104: $ RETURN
105: *
106: * Check that the diagonal matrix D is nonsingular.
107: *
108: IF( UPPER ) THEN
109: *
110: * Upper triangular storage: examine D from bottom to top
111: *
112: KP = N*( N+1 ) / 2
113: DO 10 INFO = N, 1, -1
114: IF( IPIV( INFO ).GT.0 .AND. AP( KP ).EQ.ZERO )
115: $ RETURN
116: KP = KP - INFO
117: 10 CONTINUE
118: ELSE
119: *
120: * Lower triangular storage: examine D from top to bottom.
121: *
122: KP = 1
123: DO 20 INFO = 1, N
124: IF( IPIV( INFO ).GT.0 .AND. AP( KP ).EQ.ZERO )
125: $ RETURN
126: KP = KP + N - INFO + 1
127: 20 CONTINUE
128: END IF
129: INFO = 0
130: *
131: IF( UPPER ) THEN
132: *
133: * Compute inv(A) from the factorization A = U*D*U'.
134: *
135: * K is the main loop index, increasing from 1 to N in steps of
136: * 1 or 2, depending on the size of the diagonal blocks.
137: *
138: K = 1
139: KC = 1
140: 30 CONTINUE
141: *
142: * If K > N, exit from loop.
143: *
144: IF( K.GT.N )
145: $ GO TO 50
146: *
147: KCNEXT = KC + K
148: IF( IPIV( K ).GT.0 ) THEN
149: *
150: * 1 x 1 diagonal block
151: *
152: * Invert the diagonal block.
153: *
154: AP( KC+K-1 ) = ONE / DBLE( AP( KC+K-1 ) )
155: *
156: * Compute column K of the inverse.
157: *
158: IF( K.GT.1 ) THEN
159: CALL ZCOPY( K-1, AP( KC ), 1, WORK, 1 )
160: CALL ZHPMV( UPLO, K-1, -CONE, AP, WORK, 1, ZERO,
161: $ AP( KC ), 1 )
162: AP( KC+K-1 ) = AP( KC+K-1 ) -
163: $ DBLE( ZDOTC( K-1, WORK, 1, AP( KC ), 1 ) )
164: END IF
165: KSTEP = 1
166: ELSE
167: *
168: * 2 x 2 diagonal block
169: *
170: * Invert the diagonal block.
171: *
172: T = ABS( AP( KCNEXT+K-1 ) )
173: AK = DBLE( AP( KC+K-1 ) ) / T
174: AKP1 = DBLE( AP( KCNEXT+K ) ) / T
175: AKKP1 = AP( KCNEXT+K-1 ) / T
176: D = T*( AK*AKP1-ONE )
177: AP( KC+K-1 ) = AKP1 / D
178: AP( KCNEXT+K ) = AK / D
179: AP( KCNEXT+K-1 ) = -AKKP1 / D
180: *
181: * Compute columns K and K+1 of the inverse.
182: *
183: IF( K.GT.1 ) THEN
184: CALL ZCOPY( K-1, AP( KC ), 1, WORK, 1 )
185: CALL ZHPMV( UPLO, K-1, -CONE, AP, WORK, 1, ZERO,
186: $ AP( KC ), 1 )
187: AP( KC+K-1 ) = AP( KC+K-1 ) -
188: $ DBLE( ZDOTC( K-1, WORK, 1, AP( KC ), 1 ) )
189: AP( KCNEXT+K-1 ) = AP( KCNEXT+K-1 ) -
190: $ ZDOTC( K-1, AP( KC ), 1, AP( KCNEXT ),
191: $ 1 )
192: CALL ZCOPY( K-1, AP( KCNEXT ), 1, WORK, 1 )
193: CALL ZHPMV( UPLO, K-1, -CONE, AP, WORK, 1, ZERO,
194: $ AP( KCNEXT ), 1 )
195: AP( KCNEXT+K ) = AP( KCNEXT+K ) -
196: $ DBLE( ZDOTC( K-1, WORK, 1, AP( KCNEXT ),
197: $ 1 ) )
198: END IF
199: KSTEP = 2
200: KCNEXT = KCNEXT + K + 1
201: END IF
202: *
203: KP = ABS( IPIV( K ) )
204: IF( KP.NE.K ) THEN
205: *
206: * Interchange rows and columns K and KP in the leading
207: * submatrix A(1:k+1,1:k+1)
208: *
209: KPC = ( KP-1 )*KP / 2 + 1
210: CALL ZSWAP( KP-1, AP( KC ), 1, AP( KPC ), 1 )
211: KX = KPC + KP - 1
212: DO 40 J = KP + 1, K - 1
213: KX = KX + J - 1
214: TEMP = DCONJG( AP( KC+J-1 ) )
215: AP( KC+J-1 ) = DCONJG( AP( KX ) )
216: AP( KX ) = TEMP
217: 40 CONTINUE
218: AP( KC+KP-1 ) = DCONJG( AP( KC+KP-1 ) )
219: TEMP = AP( KC+K-1 )
220: AP( KC+K-1 ) = AP( KPC+KP-1 )
221: AP( KPC+KP-1 ) = TEMP
222: IF( KSTEP.EQ.2 ) THEN
223: TEMP = AP( KC+K+K-1 )
224: AP( KC+K+K-1 ) = AP( KC+K+KP-1 )
225: AP( KC+K+KP-1 ) = TEMP
226: END IF
227: END IF
228: *
229: K = K + KSTEP
230: KC = KCNEXT
231: GO TO 30
232: 50 CONTINUE
233: *
234: ELSE
235: *
236: * Compute inv(A) from the factorization A = L*D*L'.
237: *
238: * K is the main loop index, increasing from 1 to N in steps of
239: * 1 or 2, depending on the size of the diagonal blocks.
240: *
241: NPP = N*( N+1 ) / 2
242: K = N
243: KC = NPP
244: 60 CONTINUE
245: *
246: * If K < 1, exit from loop.
247: *
248: IF( K.LT.1 )
249: $ GO TO 80
250: *
251: KCNEXT = KC - ( N-K+2 )
252: IF( IPIV( K ).GT.0 ) THEN
253: *
254: * 1 x 1 diagonal block
255: *
256: * Invert the diagonal block.
257: *
258: AP( KC ) = ONE / DBLE( AP( KC ) )
259: *
260: * Compute column K of the inverse.
261: *
262: IF( K.LT.N ) THEN
263: CALL ZCOPY( N-K, AP( KC+1 ), 1, WORK, 1 )
264: CALL ZHPMV( UPLO, N-K, -CONE, AP( KC+N-K+1 ), WORK, 1,
265: $ ZERO, AP( KC+1 ), 1 )
266: AP( KC ) = AP( KC ) - DBLE( ZDOTC( N-K, WORK, 1,
267: $ AP( KC+1 ), 1 ) )
268: END IF
269: KSTEP = 1
270: ELSE
271: *
272: * 2 x 2 diagonal block
273: *
274: * Invert the diagonal block.
275: *
276: T = ABS( AP( KCNEXT+1 ) )
277: AK = DBLE( AP( KCNEXT ) ) / T
278: AKP1 = DBLE( AP( KC ) ) / T
279: AKKP1 = AP( KCNEXT+1 ) / T
280: D = T*( AK*AKP1-ONE )
281: AP( KCNEXT ) = AKP1 / D
282: AP( KC ) = AK / D
283: AP( KCNEXT+1 ) = -AKKP1 / D
284: *
285: * Compute columns K-1 and K of the inverse.
286: *
287: IF( K.LT.N ) THEN
288: CALL ZCOPY( N-K, AP( KC+1 ), 1, WORK, 1 )
289: CALL ZHPMV( UPLO, N-K, -CONE, AP( KC+( N-K+1 ) ), WORK,
290: $ 1, ZERO, AP( KC+1 ), 1 )
291: AP( KC ) = AP( KC ) - DBLE( ZDOTC( N-K, WORK, 1,
292: $ AP( KC+1 ), 1 ) )
293: AP( KCNEXT+1 ) = AP( KCNEXT+1 ) -
294: $ ZDOTC( N-K, AP( KC+1 ), 1,
295: $ AP( KCNEXT+2 ), 1 )
296: CALL ZCOPY( N-K, AP( KCNEXT+2 ), 1, WORK, 1 )
297: CALL ZHPMV( UPLO, N-K, -CONE, AP( KC+( N-K+1 ) ), WORK,
298: $ 1, ZERO, AP( KCNEXT+2 ), 1 )
299: AP( KCNEXT ) = AP( KCNEXT ) -
300: $ DBLE( ZDOTC( N-K, WORK, 1, AP( KCNEXT+2 ),
301: $ 1 ) )
302: END IF
303: KSTEP = 2
304: KCNEXT = KCNEXT - ( N-K+3 )
305: END IF
306: *
307: KP = ABS( IPIV( K ) )
308: IF( KP.NE.K ) THEN
309: *
310: * Interchange rows and columns K and KP in the trailing
311: * submatrix A(k-1:n,k-1:n)
312: *
313: KPC = NPP - ( N-KP+1 )*( N-KP+2 ) / 2 + 1
314: IF( KP.LT.N )
315: $ CALL ZSWAP( N-KP, AP( KC+KP-K+1 ), 1, AP( KPC+1 ), 1 )
316: KX = KC + KP - K
317: DO 70 J = K + 1, KP - 1
318: KX = KX + N - J + 1
319: TEMP = DCONJG( AP( KC+J-K ) )
320: AP( KC+J-K ) = DCONJG( AP( KX ) )
321: AP( KX ) = TEMP
322: 70 CONTINUE
323: AP( KC+KP-K ) = DCONJG( AP( KC+KP-K ) )
324: TEMP = AP( KC )
325: AP( KC ) = AP( KPC )
326: AP( KPC ) = TEMP
327: IF( KSTEP.EQ.2 ) THEN
328: TEMP = AP( KC-N+K-1 )
329: AP( KC-N+K-1 ) = AP( KC-N+KP-1 )
330: AP( KC-N+KP-1 ) = TEMP
331: END IF
332: END IF
333: *
334: K = K - KSTEP
335: KC = KCNEXT
336: GO TO 60
337: 80 CONTINUE
338: END IF
339: *
340: RETURN
341: *
342: * End of ZHPTRI
343: *
344: END