LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dsytri_rook.f
Go to the documentation of this file.
1*> \brief \b DSYTRI_ROOK
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DSYTRI_ROOK + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsytri_rook.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsytri_rook.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsytri_rook.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DSYTRI_ROOK( 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_ROOK computes the inverse of a real symmetric
39*> matrix A using the factorization A = U*D*U**T or A = L*D*L**T
40*> computed by DSYTRF_ROOK.
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_ROOK.
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_ROOK.
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_rook
111*
112*> \par Contributors:
113* ==================
114*>
115*> \verbatim
116*>
117*> April 2012, Igor Kozachenko,
118*> Computer Science Division,
119*> University of California, Berkeley
120*>
121*> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
122*> School of Mathematics,
123*> University of Manchester
124*>
125*> \endverbatim
126*
127* =====================================================================
128 SUBROUTINE dsytri_rook( UPLO, N, A, LDA, IPIV, WORK, INFO )
129*
130* -- LAPACK computational routine --
131* -- LAPACK is a software package provided by Univ. of Tennessee, --
132* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
133*
134* .. Scalar Arguments ..
135 CHARACTER UPLO
136 INTEGER INFO, LDA, N
137* ..
138* .. Array Arguments ..
139 INTEGER IPIV( * )
140 DOUBLE PRECISION A( LDA, * ), WORK( * )
141* ..
142*
143* =====================================================================
144*
145* .. Parameters ..
146 DOUBLE PRECISION ONE, ZERO
147 parameter( one = 1.0d+0, zero = 0.0d+0 )
148* ..
149* .. Local Scalars ..
150 LOGICAL UPPER
151 INTEGER K, KP, KSTEP
152 DOUBLE PRECISION AK, AKKP1, AKP1, D, T, TEMP
153* ..
154* .. External Functions ..
155 LOGICAL LSAME
156 DOUBLE PRECISION DDOT
157 EXTERNAL lsame, ddot
158* ..
159* .. External Subroutines ..
160 EXTERNAL dcopy, dswap, dsymv, xerbla
161* ..
162* .. Intrinsic Functions ..
163 INTRINSIC abs, max
164* ..
165* .. Executable Statements ..
166*
167* Test the input parameters.
168*
169 info = 0
170 upper = lsame( uplo, 'U' )
171 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
172 info = -1
173 ELSE IF( n.LT.0 ) THEN
174 info = -2
175 ELSE IF( lda.LT.max( 1, n ) ) THEN
176 info = -4
177 END IF
178 IF( info.NE.0 ) THEN
179 CALL xerbla( 'DSYTRI_ROOK', -info )
180 RETURN
181 END IF
182*
183* Quick return if possible
184*
185 IF( n.EQ.0 )
186 $ RETURN
187*
188* Check that the diagonal matrix D is nonsingular.
189*
190 IF( upper ) THEN
191*
192* Upper triangular storage: examine D from bottom to top
193*
194 DO 10 info = n, 1, -1
195 IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.zero )
196 $ RETURN
197 10 CONTINUE
198 ELSE
199*
200* Lower triangular storage: examine D from top to bottom.
201*
202 DO 20 info = 1, n
203 IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.zero )
204 $ RETURN
205 20 CONTINUE
206 END IF
207 info = 0
208*
209 IF( upper ) THEN
210*
211* Compute inv(A) from the factorization A = U*D*U**T.
212*
213* K is the main loop index, increasing from 1 to N in steps of
214* 1 or 2, depending on the size of the diagonal blocks.
215*
216 k = 1
217 30 CONTINUE
218*
219* If K > N, exit from loop.
220*
221 IF( k.GT.n )
222 $ GO TO 40
223*
224 IF( ipiv( k ).GT.0 ) THEN
225*
226* 1 x 1 diagonal block
227*
228* Invert the diagonal block.
229*
230 a( k, k ) = one / a( k, k )
231*
232* Compute column K of the inverse.
233*
234 IF( k.GT.1 ) THEN
235 CALL dcopy( k-1, a( 1, k ), 1, work, 1 )
236 CALL dsymv( uplo, k-1, -one, a, lda, work, 1, zero,
237 $ a( 1, k ), 1 )
238 a( k, k ) = a( k, k ) - ddot( k-1, work, 1, a( 1, k ),
239 $ 1 )
240 END IF
241 kstep = 1
242 ELSE
243*
244* 2 x 2 diagonal block
245*
246* Invert the diagonal block.
247*
248 t = abs( a( k, k+1 ) )
249 ak = a( k, k ) / t
250 akp1 = a( k+1, k+1 ) / t
251 akkp1 = a( k, k+1 ) / t
252 d = t*( ak*akp1-one )
253 a( k, k ) = akp1 / d
254 a( k+1, k+1 ) = ak / d
255 a( k, k+1 ) = -akkp1 / d
256*
257* Compute columns K and K+1 of the inverse.
258*
259 IF( k.GT.1 ) THEN
260 CALL dcopy( k-1, a( 1, k ), 1, work, 1 )
261 CALL dsymv( uplo, k-1, -one, a, lda, work, 1, zero,
262 $ a( 1, k ), 1 )
263 a( k, k ) = a( k, k ) - ddot( k-1, work, 1, a( 1, k ),
264 $ 1 )
265 a( k, k+1 ) = a( k, k+1 ) -
266 $ ddot( k-1, a( 1, k ), 1, a( 1, k+1 ), 1 )
267 CALL dcopy( k-1, a( 1, k+1 ), 1, work, 1 )
268 CALL dsymv( uplo, k-1, -one, a, lda, work, 1, zero,
269 $ a( 1, k+1 ), 1 )
270 a( k+1, k+1 ) = a( k+1, k+1 ) -
271 $ ddot( k-1, work, 1, a( 1, k+1 ), 1 )
272 END IF
273 kstep = 2
274 END IF
275*
276 IF( kstep.EQ.1 ) THEN
277*
278* Interchange rows and columns K and IPIV(K) in the leading
279* submatrix A(1:k+1,1:k+1)
280*
281 kp = ipiv( k )
282 IF( kp.NE.k ) THEN
283 IF( kp.GT.1 )
284 $ CALL dswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
285 CALL dswap( k-kp-1, a( kp+1, k ), 1, a( kp, kp+1 ), lda )
286 temp = a( k, k )
287 a( k, k ) = a( kp, kp )
288 a( kp, kp ) = temp
289 END IF
290 ELSE
291*
292* Interchange rows and columns K and K+1 with -IPIV(K) and
293* -IPIV(K+1)in the leading submatrix A(1:k+1,1:k+1)
294*
295 kp = -ipiv( k )
296 IF( kp.NE.k ) THEN
297 IF( kp.GT.1 )
298 $ CALL dswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
299 CALL dswap( k-kp-1, a( kp+1, k ), 1, a( kp, kp+1 ), lda )
300*
301 temp = a( k, k )
302 a( k, k ) = a( kp, kp )
303 a( kp, kp ) = temp
304 temp = a( k, k+1 )
305 a( k, k+1 ) = a( kp, k+1 )
306 a( kp, k+1 ) = temp
307 END IF
308*
309 k = k + 1
310 kp = -ipiv( k )
311 IF( kp.NE.k ) THEN
312 IF( kp.GT.1 )
313 $ CALL dswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
314 CALL dswap( k-kp-1, a( kp+1, k ), 1, a( kp, kp+1 ), lda )
315 temp = a( k, k )
316 a( k, k ) = a( kp, kp )
317 a( kp, kp ) = temp
318 END IF
319 END IF
320*
321 k = k + 1
322 GO TO 30
323 40 CONTINUE
324*
325 ELSE
326*
327* Compute inv(A) from the factorization A = L*D*L**T.
328*
329* K is the main loop index, increasing from 1 to N in steps of
330* 1 or 2, depending on the size of the diagonal blocks.
331*
332 k = n
333 50 CONTINUE
334*
335* If K < 1, exit from loop.
336*
337 IF( k.LT.1 )
338 $ GO TO 60
339*
340 IF( ipiv( k ).GT.0 ) THEN
341*
342* 1 x 1 diagonal block
343*
344* Invert the diagonal block.
345*
346 a( k, k ) = one / a( k, k )
347*
348* Compute column K of the inverse.
349*
350 IF( k.LT.n ) THEN
351 CALL dcopy( n-k, a( k+1, k ), 1, work, 1 )
352 CALL dsymv( uplo, n-k, -one, a( k+1, k+1 ), lda, work, 1,
353 $ zero, a( k+1, k ), 1 )
354 a( k, k ) = a( k, k ) - ddot( n-k, work, 1, a( k+1, k ),
355 $ 1 )
356 END IF
357 kstep = 1
358 ELSE
359*
360* 2 x 2 diagonal block
361*
362* Invert the diagonal block.
363*
364 t = abs( a( k, k-1 ) )
365 ak = a( k-1, k-1 ) / t
366 akp1 = a( k, k ) / t
367 akkp1 = a( k, k-1 ) / t
368 d = t*( ak*akp1-one )
369 a( k-1, k-1 ) = akp1 / d
370 a( k, k ) = ak / d
371 a( k, k-1 ) = -akkp1 / d
372*
373* Compute columns K-1 and K of the inverse.
374*
375 IF( k.LT.n ) THEN
376 CALL dcopy( n-k, a( k+1, k ), 1, work, 1 )
377 CALL dsymv( uplo, n-k, -one, a( k+1, k+1 ), lda, work, 1,
378 $ zero, a( k+1, k ), 1 )
379 a( k, k ) = a( k, k ) - ddot( n-k, work, 1, a( k+1, k ),
380 $ 1 )
381 a( k, k-1 ) = a( k, k-1 ) -
382 $ ddot( n-k, a( k+1, k ), 1, a( k+1, k-1 ),
383 $ 1 )
384 CALL dcopy( n-k, a( k+1, k-1 ), 1, work, 1 )
385 CALL dsymv( uplo, n-k, -one, a( k+1, k+1 ), lda, work, 1,
386 $ zero, a( k+1, k-1 ), 1 )
387 a( k-1, k-1 ) = a( k-1, k-1 ) -
388 $ ddot( n-k, work, 1, a( k+1, k-1 ), 1 )
389 END IF
390 kstep = 2
391 END IF
392*
393 IF( kstep.EQ.1 ) THEN
394*
395* Interchange rows and columns K and IPIV(K) in the trailing
396* submatrix A(k-1:n,k-1:n)
397*
398 kp = ipiv( k )
399 IF( kp.NE.k ) THEN
400 IF( kp.LT.n )
401 $ CALL dswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
402 CALL dswap( kp-k-1, a( k+1, k ), 1, a( kp, k+1 ), lda )
403 temp = a( k, k )
404 a( k, k ) = a( kp, kp )
405 a( kp, kp ) = temp
406 END IF
407 ELSE
408*
409* Interchange rows and columns K and K-1 with -IPIV(K) and
410* -IPIV(K-1) in the trailing submatrix A(k-1:n,k-1:n)
411*
412 kp = -ipiv( k )
413 IF( kp.NE.k ) THEN
414 IF( kp.LT.n )
415 $ CALL dswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
416 CALL dswap( kp-k-1, a( k+1, k ), 1, a( kp, k+1 ), lda )
417*
418 temp = a( k, k )
419 a( k, k ) = a( kp, kp )
420 a( kp, kp ) = temp
421 temp = a( k, k-1 )
422 a( k, k-1 ) = a( kp, k-1 )
423 a( kp, k-1 ) = temp
424 END IF
425*
426 k = k - 1
427 kp = -ipiv( k )
428 IF( kp.NE.k ) THEN
429 IF( kp.LT.n )
430 $ CALL dswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
431 CALL dswap( kp-k-1, a( k+1, k ), 1, a( kp, k+1 ), lda )
432 temp = a( k, k )
433 a( k, k ) = a( kp, kp )
434 a( kp, kp ) = temp
435 END IF
436 END IF
437*
438 k = k - 1
439 GO TO 50
440 60 CONTINUE
441 END IF
442*
443 RETURN
444*
445* End of DSYTRI_ROOK
446*
447 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_rook(uplo, n, a, lda, ipiv, work, info)
DSYTRI_ROOK
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82