LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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 doubleSYcomputational
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)
XERBLA
Definition: xerbla.f:60
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.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
Definition: dsytri_rook.f:129