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