1: *> \brief \b DSPTRI
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
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">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DSPTRI( 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: * DOUBLE PRECISION AP( * ), WORK( * )
30: * ..
31: *
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: *
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 doubleOTHERcomputational
108: *
109: * =====================================================================
110: SUBROUTINE DSPTRI( 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: DOUBLE PRECISION AP( * ), WORK( * )
124: * ..
125: *
126: * =====================================================================
127: *
128: * .. Parameters ..
129: DOUBLE PRECISION ONE, ZERO
130: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
131: * ..
132: * .. Local Scalars ..
133: LOGICAL UPPER
134: INTEGER J, K, KC, KCNEXT, KP, KPC, KSTEP, KX, NPP
135: DOUBLE PRECISION AK, AKKP1, AKP1, D, T, TEMP
136: * ..
137: * .. External Functions ..
138: LOGICAL LSAME
139: DOUBLE PRECISION DDOT
140: EXTERNAL LSAME, DDOT
141: * ..
142: * .. External Subroutines ..
143: EXTERNAL DCOPY, DSPMV, DSWAP, XERBLA
144: * ..
145: * .. Intrinsic Functions ..
146: INTRINSIC ABS
147: * ..
148: * .. Executable Statements ..
149: *
150: * Test the input parameters.
151: *
152: INFO = 0
153: UPPER = LSAME( UPLO, 'U' )
154: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
155: INFO = -1
156: ELSE IF( N.LT.0 ) THEN
157: INFO = -2
158: END IF
159: IF( INFO.NE.0 ) THEN
160: CALL XERBLA( 'DSPTRI', -INFO )
161: RETURN
162: END IF
163: *
164: * Quick return if possible
165: *
166: IF( N.EQ.0 )
167: $ RETURN
168: *
169: * Check that the diagonal matrix D is nonsingular.
170: *
171: IF( UPPER ) THEN
172: *
173: * Upper triangular storage: examine D from bottom to top
174: *
175: KP = N*( N+1 ) / 2
176: DO 10 INFO = N, 1, -1
177: IF( IPIV( INFO ).GT.0 .AND. AP( KP ).EQ.ZERO )
178: $ RETURN
179: KP = KP - INFO
180: 10 CONTINUE
181: ELSE
182: *
183: * Lower triangular storage: examine D from top to bottom.
184: *
185: KP = 1
186: DO 20 INFO = 1, N
187: IF( IPIV( INFO ).GT.0 .AND. AP( KP ).EQ.ZERO )
188: $ RETURN
189: KP = KP + N - INFO + 1
190: 20 CONTINUE
191: END IF
192: INFO = 0
193: *
194: IF( UPPER ) THEN
195: *
196: * Compute inv(A) from the factorization A = U*D*U**T.
197: *
198: * K is the main loop index, increasing from 1 to N in steps of
199: * 1 or 2, depending on the size of the diagonal blocks.
200: *
201: K = 1
202: KC = 1
203: 30 CONTINUE
204: *
205: * If K > N, exit from loop.
206: *
207: IF( K.GT.N )
208: $ GO TO 50
209: *
210: KCNEXT = KC + K
211: IF( IPIV( K ).GT.0 ) THEN
212: *
213: * 1 x 1 diagonal block
214: *
215: * Invert the diagonal block.
216: *
217: AP( KC+K-1 ) = ONE / AP( KC+K-1 )
218: *
219: * Compute column K of the inverse.
220: *
221: IF( K.GT.1 ) THEN
222: CALL DCOPY( K-1, AP( KC ), 1, WORK, 1 )
223: CALL DSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO, AP( KC ),
224: $ 1 )
225: AP( KC+K-1 ) = AP( KC+K-1 ) -
226: $ DDOT( K-1, WORK, 1, AP( KC ), 1 )
227: END IF
228: KSTEP = 1
229: ELSE
230: *
231: * 2 x 2 diagonal block
232: *
233: * Invert the diagonal block.
234: *
235: T = ABS( AP( KCNEXT+K-1 ) )
236: AK = AP( KC+K-1 ) / T
237: AKP1 = AP( KCNEXT+K ) / T
238: AKKP1 = AP( KCNEXT+K-1 ) / T
239: D = T*( AK*AKP1-ONE )
240: AP( KC+K-1 ) = AKP1 / D
241: AP( KCNEXT+K ) = AK / D
242: AP( KCNEXT+K-1 ) = -AKKP1 / D
243: *
244: * Compute columns K and K+1 of the inverse.
245: *
246: IF( K.GT.1 ) THEN
247: CALL DCOPY( K-1, AP( KC ), 1, WORK, 1 )
248: CALL DSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO, AP( KC ),
249: $ 1 )
250: AP( KC+K-1 ) = AP( KC+K-1 ) -
251: $ DDOT( K-1, WORK, 1, AP( KC ), 1 )
252: AP( KCNEXT+K-1 ) = AP( KCNEXT+K-1 ) -
253: $ DDOT( K-1, AP( KC ), 1, AP( KCNEXT ),
254: $ 1 )
255: CALL DCOPY( K-1, AP( KCNEXT ), 1, WORK, 1 )
256: CALL DSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO,
257: $ AP( KCNEXT ), 1 )
258: AP( KCNEXT+K ) = AP( KCNEXT+K ) -
259: $ DDOT( K-1, WORK, 1, AP( KCNEXT ), 1 )
260: END IF
261: KSTEP = 2
262: KCNEXT = KCNEXT + K + 1
263: END IF
264: *
265: KP = ABS( IPIV( K ) )
266: IF( KP.NE.K ) THEN
267: *
268: * Interchange rows and columns K and KP in the leading
269: * submatrix A(1:k+1,1:k+1)
270: *
271: KPC = ( KP-1 )*KP / 2 + 1
272: CALL DSWAP( KP-1, AP( KC ), 1, AP( KPC ), 1 )
273: KX = KPC + KP - 1
274: DO 40 J = KP + 1, K - 1
275: KX = KX + J - 1
276: TEMP = AP( KC+J-1 )
277: AP( KC+J-1 ) = AP( KX )
278: AP( KX ) = TEMP
279: 40 CONTINUE
280: TEMP = AP( KC+K-1 )
281: AP( KC+K-1 ) = AP( KPC+KP-1 )
282: AP( KPC+KP-1 ) = TEMP
283: IF( KSTEP.EQ.2 ) THEN
284: TEMP = AP( KC+K+K-1 )
285: AP( KC+K+K-1 ) = AP( KC+K+KP-1 )
286: AP( KC+K+KP-1 ) = TEMP
287: END IF
288: END IF
289: *
290: K = K + KSTEP
291: KC = KCNEXT
292: GO TO 30
293: 50 CONTINUE
294: *
295: ELSE
296: *
297: * Compute inv(A) from the factorization A = L*D*L**T.
298: *
299: * K is the main loop index, increasing from 1 to N in steps of
300: * 1 or 2, depending on the size of the diagonal blocks.
301: *
302: NPP = N*( N+1 ) / 2
303: K = N
304: KC = NPP
305: 60 CONTINUE
306: *
307: * If K < 1, exit from loop.
308: *
309: IF( K.LT.1 )
310: $ GO TO 80
311: *
312: KCNEXT = KC - ( N-K+2 )
313: IF( IPIV( K ).GT.0 ) THEN
314: *
315: * 1 x 1 diagonal block
316: *
317: * Invert the diagonal block.
318: *
319: AP( KC ) = ONE / AP( KC )
320: *
321: * Compute column K of the inverse.
322: *
323: IF( K.LT.N ) THEN
324: CALL DCOPY( N-K, AP( KC+1 ), 1, WORK, 1 )
325: CALL DSPMV( UPLO, N-K, -ONE, AP( KC+N-K+1 ), WORK, 1,
326: $ ZERO, AP( KC+1 ), 1 )
327: AP( KC ) = AP( KC ) - DDOT( N-K, WORK, 1, AP( KC+1 ), 1 )
328: END IF
329: KSTEP = 1
330: ELSE
331: *
332: * 2 x 2 diagonal block
333: *
334: * Invert the diagonal block.
335: *
336: T = ABS( AP( KCNEXT+1 ) )
337: AK = AP( KCNEXT ) / T
338: AKP1 = AP( KC ) / T
339: AKKP1 = AP( KCNEXT+1 ) / T
340: D = T*( AK*AKP1-ONE )
341: AP( KCNEXT ) = AKP1 / D
342: AP( KC ) = AK / D
343: AP( KCNEXT+1 ) = -AKKP1 / D
344: *
345: * Compute columns K-1 and K of the inverse.
346: *
347: IF( K.LT.N ) THEN
348: CALL DCOPY( N-K, AP( KC+1 ), 1, WORK, 1 )
349: CALL DSPMV( UPLO, N-K, -ONE, AP( KC+( N-K+1 ) ), WORK, 1,
350: $ ZERO, AP( KC+1 ), 1 )
351: AP( KC ) = AP( KC ) - DDOT( N-K, WORK, 1, AP( KC+1 ), 1 )
352: AP( KCNEXT+1 ) = AP( KCNEXT+1 ) -
353: $ DDOT( N-K, AP( KC+1 ), 1,
354: $ AP( KCNEXT+2 ), 1 )
355: CALL DCOPY( N-K, AP( KCNEXT+2 ), 1, WORK, 1 )
356: CALL DSPMV( UPLO, N-K, -ONE, AP( KC+( N-K+1 ) ), WORK, 1,
357: $ ZERO, AP( KCNEXT+2 ), 1 )
358: AP( KCNEXT ) = AP( KCNEXT ) -
359: $ DDOT( 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 DSWAP( 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 DSPTRI
400: *
401: 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