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