LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dsptri.f
Go to the documentation of this file.
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*> \ingroup hptri
106*
107* =====================================================================
108 SUBROUTINE dsptri( UPLO, N, AP, IPIV, WORK, INFO )
109*
110* -- LAPACK computational routine --
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*
193* Compute inv(A) from the factorization A = U*D*U**T.
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*
294* Compute inv(A) from the factorization A = L*D*L**T.
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
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
subroutine dspmv(uplo, n, alpha, ap, x, incx, beta, y, incy)
DSPMV
Definition dspmv.f:147
subroutine dsptri(uplo, n, ap, ipiv, work, info)
DSPTRI
Definition dsptri.f:109
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82