LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dsytri.f
Go to the documentation of this file.
1*> \brief \b DSYTRI
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DSYTRI + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsytri.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsytri.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsytri.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DSYTRI( UPLO, N, A, LDA, IPIV, WORK, INFO )
22*
23* .. Scalar Arguments ..
24* CHARACTER UPLO
25* INTEGER INFO, LDA, N
26* ..
27* .. Array Arguments ..
28* INTEGER IPIV( * )
29* DOUBLE PRECISION A( LDA, * ), WORK( * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> DSYTRI computes the inverse of a real symmetric indefinite matrix
39*> A using the factorization A = U*D*U**T or A = L*D*L**T computed by
40*> DSYTRF.
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] A
62*> \verbatim
63*> A is DOUBLE PRECISION array, dimension (LDA,N)
64*> On entry, the block diagonal matrix D and the multipliers
65*> used to obtain the factor U or L as computed by DSYTRF.
66*>
67*> On exit, if INFO = 0, the (symmetric) inverse of the original
68*> matrix. If UPLO = 'U', the upper triangular part of the
69*> inverse is formed and the part of A below the diagonal is not
70*> referenced; if UPLO = 'L' the lower triangular part of the
71*> inverse is formed and the part of A above the diagonal is
72*> not referenced.
73*> \endverbatim
74*>
75*> \param[in] LDA
76*> \verbatim
77*> LDA is INTEGER
78*> The leading dimension of the array A. LDA >= max(1,N).
79*> \endverbatim
80*>
81*> \param[in] IPIV
82*> \verbatim
83*> IPIV is INTEGER array, dimension (N)
84*> Details of the interchanges and the block structure of D
85*> as determined by DSYTRF.
86*> \endverbatim
87*>
88*> \param[out] WORK
89*> \verbatim
90*> WORK is DOUBLE PRECISION array, dimension (N)
91*> \endverbatim
92*>
93*> \param[out] INFO
94*> \verbatim
95*> INFO is INTEGER
96*> = 0: successful exit
97*> < 0: if INFO = -i, the i-th argument had an illegal value
98*> > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
99*> inverse could not be computed.
100*> \endverbatim
101*
102* Authors:
103* ========
104*
105*> \author Univ. of Tennessee
106*> \author Univ. of California Berkeley
107*> \author Univ. of Colorado Denver
108*> \author NAG Ltd.
109*
110*> \ingroup hetri
111*
112* =====================================================================
113 SUBROUTINE dsytri( UPLO, N, A, LDA, IPIV, WORK, INFO )
114*
115* -- LAPACK computational routine --
116* -- LAPACK is a software package provided by Univ. of Tennessee, --
117* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
118*
119* .. Scalar Arguments ..
120 CHARACTER UPLO
121 INTEGER INFO, LDA, N
122* ..
123* .. Array Arguments ..
124 INTEGER IPIV( * )
125 DOUBLE PRECISION A( LDA, * ), WORK( * )
126* ..
127*
128* =====================================================================
129*
130* .. Parameters ..
131 DOUBLE PRECISION ONE, ZERO
132 parameter( one = 1.0d+0, zero = 0.0d+0 )
133* ..
134* .. Local Scalars ..
135 LOGICAL UPPER
136 INTEGER K, KP, KSTEP
137 DOUBLE PRECISION AK, AKKP1, AKP1, D, T, TEMP
138* ..
139* .. External Functions ..
140 LOGICAL LSAME
141 DOUBLE PRECISION DDOT
142 EXTERNAL lsame, ddot
143* ..
144* .. External Subroutines ..
145 EXTERNAL dcopy, dswap, dsymv, xerbla
146* ..
147* .. Intrinsic Functions ..
148 INTRINSIC abs, max
149* ..
150* .. Executable Statements ..
151*
152* Test the input parameters.
153*
154 info = 0
155 upper = lsame( uplo, 'U' )
156 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
157 info = -1
158 ELSE IF( n.LT.0 ) THEN
159 info = -2
160 ELSE IF( lda.LT.max( 1, n ) ) THEN
161 info = -4
162 END IF
163 IF( info.NE.0 ) THEN
164 CALL xerbla( 'DSYTRI', -info )
165 RETURN
166 END IF
167*
168* Quick return if possible
169*
170 IF( n.EQ.0 )
171 $ RETURN
172*
173* Check that the diagonal matrix D is nonsingular.
174*
175 IF( upper ) THEN
176*
177* Upper triangular storage: examine D from bottom to top
178*
179 DO 10 info = n, 1, -1
180 IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.zero )
181 $ RETURN
182 10 CONTINUE
183 ELSE
184*
185* Lower triangular storage: examine D from top to bottom.
186*
187 DO 20 info = 1, n
188 IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.zero )
189 $ RETURN
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 30 CONTINUE
203*
204* If K > N, exit from loop.
205*
206 IF( k.GT.n )
207 $ GO TO 40
208*
209 IF( ipiv( k ).GT.0 ) THEN
210*
211* 1 x 1 diagonal block
212*
213* Invert the diagonal block.
214*
215 a( k, k ) = one / a( k, k )
216*
217* Compute column K of the inverse.
218*
219 IF( k.GT.1 ) THEN
220 CALL dcopy( k-1, a( 1, k ), 1, work, 1 )
221 CALL dsymv( uplo, k-1, -one, a, lda, work, 1, zero,
222 $ a( 1, k ), 1 )
223 a( k, k ) = a( k, k ) - ddot( k-1, work, 1, a( 1, k ),
224 $ 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 = abs( a( k, k+1 ) )
234 ak = a( k, k ) / t
235 akp1 = a( k+1, k+1 ) / t
236 akkp1 = a( k, k+1 ) / t
237 d = t*( ak*akp1-one )
238 a( k, k ) = akp1 / d
239 a( k+1, k+1 ) = ak / d
240 a( k, k+1 ) = -akkp1 / d
241*
242* Compute columns K and K+1 of the inverse.
243*
244 IF( k.GT.1 ) THEN
245 CALL dcopy( k-1, a( 1, k ), 1, work, 1 )
246 CALL dsymv( uplo, k-1, -one, a, lda, work, 1, zero,
247 $ a( 1, k ), 1 )
248 a( k, k ) = a( k, k ) - ddot( k-1, work, 1, a( 1, k ),
249 $ 1 )
250 a( k, k+1 ) = a( k, k+1 ) -
251 $ ddot( k-1, a( 1, k ), 1, a( 1, k+1 ), 1 )
252 CALL dcopy( k-1, a( 1, k+1 ), 1, work, 1 )
253 CALL dsymv( uplo, k-1, -one, a, lda, work, 1, zero,
254 $ a( 1, k+1 ), 1 )
255 a( k+1, k+1 ) = a( k+1, k+1 ) -
256 $ ddot( k-1, work, 1, a( 1, k+1 ), 1 )
257 END IF
258 kstep = 2
259 END IF
260*
261 kp = abs( ipiv( k ) )
262 IF( kp.NE.k ) THEN
263*
264* Interchange rows and columns K and KP in the leading
265* submatrix A(1:k+1,1:k+1)
266*
267 CALL dswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
268 CALL dswap( k-kp-1, a( kp+1, k ), 1, a( kp, kp+1 ), lda )
269 temp = a( k, k )
270 a( k, k ) = a( kp, kp )
271 a( kp, kp ) = temp
272 IF( kstep.EQ.2 ) THEN
273 temp = a( k, k+1 )
274 a( k, k+1 ) = a( kp, k+1 )
275 a( kp, k+1 ) = temp
276 END IF
277 END IF
278*
279 k = k + kstep
280 GO TO 30
281 40 CONTINUE
282*
283 ELSE
284*
285* Compute inv(A) from the factorization A = L*D*L**T.
286*
287* K is the main loop index, increasing from 1 to N in steps of
288* 1 or 2, depending on the size of the diagonal blocks.
289*
290 k = n
291 50 CONTINUE
292*
293* If K < 1, exit from loop.
294*
295 IF( k.LT.1 )
296 $ GO TO 60
297*
298 IF( ipiv( k ).GT.0 ) THEN
299*
300* 1 x 1 diagonal block
301*
302* Invert the diagonal block.
303*
304 a( k, k ) = one / a( k, k )
305*
306* Compute column K of the inverse.
307*
308 IF( k.LT.n ) THEN
309 CALL dcopy( n-k, a( k+1, k ), 1, work, 1 )
310 CALL dsymv( uplo, n-k, -one, a( k+1, k+1 ), lda, work, 1,
311 $ zero, a( k+1, k ), 1 )
312 a( k, k ) = a( k, k ) - ddot( n-k, work, 1, a( k+1, k ),
313 $ 1 )
314 END IF
315 kstep = 1
316 ELSE
317*
318* 2 x 2 diagonal block
319*
320* Invert the diagonal block.
321*
322 t = abs( a( k, k-1 ) )
323 ak = a( k-1, k-1 ) / t
324 akp1 = a( k, k ) / t
325 akkp1 = a( k, k-1 ) / t
326 d = t*( ak*akp1-one )
327 a( k-1, k-1 ) = akp1 / d
328 a( k, k ) = ak / d
329 a( k, k-1 ) = -akkp1 / d
330*
331* Compute columns K-1 and K of the inverse.
332*
333 IF( k.LT.n ) THEN
334 CALL dcopy( n-k, a( k+1, k ), 1, work, 1 )
335 CALL dsymv( uplo, n-k, -one, a( k+1, k+1 ), lda, work, 1,
336 $ zero, a( k+1, k ), 1 )
337 a( k, k ) = a( k, k ) - ddot( n-k, work, 1, a( k+1, k ),
338 $ 1 )
339 a( k, k-1 ) = a( k, k-1 ) -
340 $ ddot( n-k, a( k+1, k ), 1, a( k+1, k-1 ),
341 $ 1 )
342 CALL dcopy( n-k, a( k+1, k-1 ), 1, work, 1 )
343 CALL dsymv( uplo, n-k, -one, a( k+1, k+1 ), lda, work, 1,
344 $ zero, a( k+1, k-1 ), 1 )
345 a( k-1, k-1 ) = a( k-1, k-1 ) -
346 $ ddot( n-k, work, 1, a( k+1, k-1 ), 1 )
347 END IF
348 kstep = 2
349 END IF
350*
351 kp = abs( ipiv( k ) )
352 IF( kp.NE.k ) THEN
353*
354* Interchange rows and columns K and KP in the trailing
355* submatrix A(k-1:n,k-1:n)
356*
357 IF( kp.LT.n )
358 $ CALL dswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
359 CALL dswap( kp-k-1, a( k+1, k ), 1, a( kp, k+1 ), lda )
360 temp = a( k, k )
361 a( k, k ) = a( kp, kp )
362 a( kp, kp ) = temp
363 IF( kstep.EQ.2 ) THEN
364 temp = a( k, k-1 )
365 a( k, k-1 ) = a( kp, k-1 )
366 a( kp, k-1 ) = temp
367 END IF
368 END IF
369*
370 k = k - kstep
371 GO TO 50
372 60 CONTINUE
373 END IF
374*
375 RETURN
376*
377* End of DSYTRI
378*
379 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
subroutine dsymv(uplo, n, alpha, a, lda, x, incx, beta, y, incy)
DSYMV
Definition dsymv.f:152
subroutine dsytri(uplo, n, a, lda, ipiv, work, info)
DSYTRI
Definition dsytri.f:114
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82