Annotation of rpl/lapack/lapack/zsptri.f, revision 1.18
1.9 bertrand 1: *> \brief \b ZSPTRI
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 ZSPTRI + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsptri.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsptri.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsptri.f">
1.9 bertrand 15: *> [TXT]</a>
1.15 bertrand 16: *> \endhtmlonly
1.9 bertrand 17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZSPTRI( 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: *> ZSPTRI computes the inverse of a complex symmetric indefinite matrix
39: *> A in packed storage using the factorization A = U*D*U**T or
40: *> A = L*D*L**T computed by ZSPTRF.
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**T;
52: *> = 'L': Lower triangular, form is A = L*D*L**T.
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 ZSPTRF,
66: *> stored as a packed triangular matrix.
67: *>
68: *> On exit, if INFO = 0, the (symmetric) 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 ZSPTRF.
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 ZSPTRI( 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: COMPLEX*16 ONE, ZERO
127: PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ),
128: $ ZERO = ( 0.0D+0, 0.0D+0 ) )
129: * ..
130: * .. Local Scalars ..
131: LOGICAL UPPER
132: INTEGER J, K, KC, KCNEXT, KP, KPC, KSTEP, KX, NPP
133: COMPLEX*16 AK, AKKP1, AKP1, D, T, TEMP
134: * ..
135: * .. External Functions ..
136: LOGICAL LSAME
137: COMPLEX*16 ZDOTU
138: EXTERNAL LSAME, ZDOTU
139: * ..
140: * .. External Subroutines ..
141: EXTERNAL XERBLA, ZCOPY, ZSPMV, ZSWAP
142: * ..
143: * .. Intrinsic Functions ..
144: INTRINSIC ABS
145: * ..
146: * .. Executable Statements ..
147: *
148: * Test the input parameters.
149: *
150: INFO = 0
151: UPPER = LSAME( UPLO, 'U' )
152: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
153: INFO = -1
154: ELSE IF( N.LT.0 ) THEN
155: INFO = -2
156: END IF
157: IF( INFO.NE.0 ) THEN
158: CALL XERBLA( 'ZSPTRI', -INFO )
159: RETURN
160: END IF
161: *
162: * Quick return if possible
163: *
164: IF( N.EQ.0 )
165: $ RETURN
166: *
167: * Check that the diagonal matrix D is nonsingular.
168: *
169: IF( UPPER ) THEN
170: *
171: * Upper triangular storage: examine D from bottom to top
172: *
173: KP = N*( N+1 ) / 2
174: DO 10 INFO = N, 1, -1
175: IF( IPIV( INFO ).GT.0 .AND. AP( KP ).EQ.ZERO )
176: $ RETURN
177: KP = KP - INFO
178: 10 CONTINUE
179: ELSE
180: *
181: * Lower triangular storage: examine D from top to bottom.
182: *
183: KP = 1
184: DO 20 INFO = 1, N
185: IF( IPIV( INFO ).GT.0 .AND. AP( KP ).EQ.ZERO )
186: $ RETURN
187: KP = KP + N - INFO + 1
188: 20 CONTINUE
189: END IF
190: INFO = 0
191: *
192: IF( UPPER ) THEN
193: *
1.8 bertrand 194: * Compute inv(A) from the factorization A = U*D*U**T.
1.1 bertrand 195: *
196: * K is the main loop index, increasing from 1 to N in steps of
197: * 1 or 2, depending on the size of the diagonal blocks.
198: *
199: K = 1
200: KC = 1
201: 30 CONTINUE
202: *
203: * If K > N, exit from loop.
204: *
205: IF( K.GT.N )
206: $ GO TO 50
207: *
208: KCNEXT = KC + K
209: IF( IPIV( K ).GT.0 ) THEN
210: *
211: * 1 x 1 diagonal block
212: *
213: * Invert the diagonal block.
214: *
215: AP( KC+K-1 ) = ONE / AP( KC+K-1 )
216: *
217: * Compute column K of the inverse.
218: *
219: IF( K.GT.1 ) THEN
220: CALL ZCOPY( K-1, AP( KC ), 1, WORK, 1 )
221: CALL ZSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO, AP( KC ),
222: $ 1 )
223: AP( KC+K-1 ) = AP( KC+K-1 ) -
224: $ ZDOTU( K-1, WORK, 1, AP( KC ), 1 )
225: END IF
226: KSTEP = 1
227: ELSE
228: *
229: * 2 x 2 diagonal block
230: *
231: * Invert the diagonal block.
232: *
233: T = AP( KCNEXT+K-1 )
234: AK = AP( KC+K-1 ) / T
235: AKP1 = AP( KCNEXT+K ) / T
236: AKKP1 = AP( KCNEXT+K-1 ) / T
237: D = T*( AK*AKP1-ONE )
238: AP( KC+K-1 ) = AKP1 / D
239: AP( KCNEXT+K ) = AK / D
240: AP( KCNEXT+K-1 ) = -AKKP1 / D
241: *
242: * Compute columns K and K+1 of the inverse.
243: *
244: IF( K.GT.1 ) THEN
245: CALL ZCOPY( K-1, AP( KC ), 1, WORK, 1 )
246: CALL ZSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO, AP( KC ),
247: $ 1 )
248: AP( KC+K-1 ) = AP( KC+K-1 ) -
249: $ ZDOTU( K-1, WORK, 1, AP( KC ), 1 )
250: AP( KCNEXT+K-1 ) = AP( KCNEXT+K-1 ) -
251: $ ZDOTU( K-1, AP( KC ), 1, AP( KCNEXT ),
252: $ 1 )
253: CALL ZCOPY( K-1, AP( KCNEXT ), 1, WORK, 1 )
254: CALL ZSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO,
255: $ AP( KCNEXT ), 1 )
256: AP( KCNEXT+K ) = AP( KCNEXT+K ) -
257: $ ZDOTU( K-1, WORK, 1, AP( KCNEXT ), 1 )
258: END IF
259: KSTEP = 2
260: KCNEXT = KCNEXT + K + 1
261: END IF
262: *
263: KP = ABS( IPIV( K ) )
264: IF( KP.NE.K ) THEN
265: *
266: * Interchange rows and columns K and KP in the leading
267: * submatrix A(1:k+1,1:k+1)
268: *
269: KPC = ( KP-1 )*KP / 2 + 1
270: CALL ZSWAP( KP-1, AP( KC ), 1, AP( KPC ), 1 )
271: KX = KPC + KP - 1
272: DO 40 J = KP + 1, K - 1
273: KX = KX + J - 1
274: TEMP = AP( KC+J-1 )
275: AP( KC+J-1 ) = AP( KX )
276: AP( KX ) = TEMP
277: 40 CONTINUE
278: TEMP = AP( KC+K-1 )
279: AP( KC+K-1 ) = AP( KPC+KP-1 )
280: AP( KPC+KP-1 ) = TEMP
281: IF( KSTEP.EQ.2 ) THEN
282: TEMP = AP( KC+K+K-1 )
283: AP( KC+K+K-1 ) = AP( KC+K+KP-1 )
284: AP( KC+K+KP-1 ) = TEMP
285: END IF
286: END IF
287: *
288: K = K + KSTEP
289: KC = KCNEXT
290: GO TO 30
291: 50 CONTINUE
292: *
293: ELSE
294: *
1.8 bertrand 295: * Compute inv(A) from the factorization A = L*D*L**T.
1.1 bertrand 296: *
297: * K is the main loop index, increasing from 1 to N in steps of
298: * 1 or 2, depending on the size of the diagonal blocks.
299: *
300: NPP = N*( N+1 ) / 2
301: K = N
302: KC = NPP
303: 60 CONTINUE
304: *
305: * If K < 1, exit from loop.
306: *
307: IF( K.LT.1 )
308: $ GO TO 80
309: *
310: KCNEXT = KC - ( N-K+2 )
311: IF( IPIV( K ).GT.0 ) THEN
312: *
313: * 1 x 1 diagonal block
314: *
315: * Invert the diagonal block.
316: *
317: AP( KC ) = ONE / AP( KC )
318: *
319: * Compute column K of the inverse.
320: *
321: IF( K.LT.N ) THEN
322: CALL ZCOPY( N-K, AP( KC+1 ), 1, WORK, 1 )
323: CALL ZSPMV( UPLO, N-K, -ONE, AP( KC+N-K+1 ), WORK, 1,
324: $ ZERO, AP( KC+1 ), 1 )
325: AP( KC ) = AP( KC ) - ZDOTU( N-K, WORK, 1, AP( KC+1 ),
326: $ 1 )
327: END IF
328: KSTEP = 1
329: ELSE
330: *
331: * 2 x 2 diagonal block
332: *
333: * Invert the diagonal block.
334: *
335: T = AP( KCNEXT+1 )
336: AK = AP( KCNEXT ) / T
337: AKP1 = AP( KC ) / T
338: AKKP1 = AP( KCNEXT+1 ) / T
339: D = T*( AK*AKP1-ONE )
340: AP( KCNEXT ) = AKP1 / D
341: AP( KC ) = AK / D
342: AP( KCNEXT+1 ) = -AKKP1 / D
343: *
344: * Compute columns K-1 and K of the inverse.
345: *
346: IF( K.LT.N ) THEN
347: CALL ZCOPY( N-K, AP( KC+1 ), 1, WORK, 1 )
348: CALL ZSPMV( UPLO, N-K, -ONE, AP( KC+( N-K+1 ) ), WORK, 1,
349: $ ZERO, AP( KC+1 ), 1 )
350: AP( KC ) = AP( KC ) - ZDOTU( N-K, WORK, 1, AP( KC+1 ),
351: $ 1 )
352: AP( KCNEXT+1 ) = AP( KCNEXT+1 ) -
353: $ ZDOTU( N-K, AP( KC+1 ), 1,
354: $ AP( KCNEXT+2 ), 1 )
355: CALL ZCOPY( N-K, AP( KCNEXT+2 ), 1, WORK, 1 )
356: CALL ZSPMV( UPLO, N-K, -ONE, AP( KC+( N-K+1 ) ), WORK, 1,
357: $ ZERO, AP( KCNEXT+2 ), 1 )
358: AP( KCNEXT ) = AP( KCNEXT ) -
359: $ ZDOTU( N-K, WORK, 1, AP( KCNEXT+2 ), 1 )
360: END IF
361: KSTEP = 2
362: KCNEXT = KCNEXT - ( N-K+3 )
363: END IF
364: *
365: KP = ABS( IPIV( K ) )
366: IF( KP.NE.K ) THEN
367: *
368: * Interchange rows and columns K and KP in the trailing
369: * submatrix A(k-1:n,k-1:n)
370: *
371: KPC = NPP - ( N-KP+1 )*( N-KP+2 ) / 2 + 1
372: IF( KP.LT.N )
373: $ CALL ZSWAP( N-KP, AP( KC+KP-K+1 ), 1, AP( KPC+1 ), 1 )
374: KX = KC + KP - K
375: DO 70 J = K + 1, KP - 1
376: KX = KX + N - J + 1
377: TEMP = AP( KC+J-K )
378: AP( KC+J-K ) = AP( KX )
379: AP( KX ) = TEMP
380: 70 CONTINUE
381: TEMP = AP( KC )
382: AP( KC ) = AP( KPC )
383: AP( KPC ) = TEMP
384: IF( KSTEP.EQ.2 ) THEN
385: TEMP = AP( KC-N+K-1 )
386: AP( KC-N+K-1 ) = AP( KC-N+KP-1 )
387: AP( KC-N+KP-1 ) = TEMP
388: END IF
389: END IF
390: *
391: K = K - KSTEP
392: KC = KCNEXT
393: GO TO 60
394: 80 CONTINUE
395: END IF
396: *
397: RETURN
398: *
399: * End of ZSPTRI
400: *
401: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zsptri.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack