LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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 *> \date November 2011
111 *
112 *> \ingroup doubleSYcomputational
113 *
114 * =====================================================================
115  SUBROUTINE dsytri( UPLO, N, A, LDA, IPIV, WORK, INFO )
116 *
117 * -- LAPACK computational routine (version 3.4.0) --
118 * -- LAPACK is a software package provided by Univ. of Tennessee, --
119 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
120 * November 2011
121 *
122 * .. Scalar Arguments ..
123  CHARACTER uplo
124  INTEGER info, lda, n
125 * ..
126 * .. Array Arguments ..
127  INTEGER ipiv( * )
128  DOUBLE PRECISION a( lda, * ), work( * )
129 * ..
130 *
131 * =====================================================================
132 *
133 * .. Parameters ..
134  DOUBLE PRECISION one, zero
135  parameter( one = 1.0d+0, zero = 0.0d+0 )
136 * ..
137 * .. Local Scalars ..
138  LOGICAL upper
139  INTEGER k, kp, kstep
140  DOUBLE PRECISION ak, akkp1, akp1, d, t, temp
141 * ..
142 * .. External Functions ..
143  LOGICAL lsame
144  DOUBLE PRECISION ddot
145  EXTERNAL lsame, ddot
146 * ..
147 * .. External Subroutines ..
148  EXTERNAL dcopy, dswap, dsymv, xerbla
149 * ..
150 * .. Intrinsic Functions ..
151  INTRINSIC abs, max
152 * ..
153 * .. Executable Statements ..
154 *
155 * Test the input parameters.
156 *
157  info = 0
158  upper = lsame( uplo, 'U' )
159  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
160  info = -1
161  ELSE IF( n.LT.0 ) THEN
162  info = -2
163  ELSE IF( lda.LT.max( 1, n ) ) THEN
164  info = -4
165  END IF
166  IF( info.NE.0 ) THEN
167  CALL xerbla( 'DSYTRI', -info )
168  return
169  END IF
170 *
171 * Quick return if possible
172 *
173  IF( n.EQ.0 )
174  $ return
175 *
176 * Check that the diagonal matrix D is nonsingular.
177 *
178  IF( upper ) THEN
179 *
180 * Upper triangular storage: examine D from bottom to top
181 *
182  DO 10 info = n, 1, -1
183  IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.zero )
184  $ return
185  10 continue
186  ELSE
187 *
188 * Lower triangular storage: examine D from top to bottom.
189 *
190  DO 20 info = 1, n
191  IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.zero )
192  $ return
193  20 continue
194  END IF
195  info = 0
196 *
197  IF( upper ) THEN
198 *
199 * Compute inv(A) from the factorization A = U*D*U**T.
200 *
201 * K is the main loop index, increasing from 1 to N in steps of
202 * 1 or 2, depending on the size of the diagonal blocks.
203 *
204  k = 1
205  30 continue
206 *
207 * If K > N, exit from loop.
208 *
209  IF( k.GT.n )
210  $ go to 40
211 *
212  IF( ipiv( k ).GT.0 ) THEN
213 *
214 * 1 x 1 diagonal block
215 *
216 * Invert the diagonal block.
217 *
218  a( k, k ) = one / a( k, k )
219 *
220 * Compute column K of the inverse.
221 *
222  IF( k.GT.1 ) THEN
223  CALL dcopy( k-1, a( 1, k ), 1, work, 1 )
224  CALL dsymv( uplo, k-1, -one, a, lda, work, 1, zero,
225  $ a( 1, k ), 1 )
226  a( k, k ) = a( k, k ) - ddot( k-1, work, 1, a( 1, k ),
227  $ 1 )
228  END IF
229  kstep = 1
230  ELSE
231 *
232 * 2 x 2 diagonal block
233 *
234 * Invert the diagonal block.
235 *
236  t = abs( a( k, k+1 ) )
237  ak = a( k, k ) / t
238  akp1 = a( k+1, k+1 ) / t
239  akkp1 = a( k, k+1 ) / t
240  d = t*( ak*akp1-one )
241  a( k, k ) = akp1 / d
242  a( k+1, k+1 ) = ak / d
243  a( k, k+1 ) = -akkp1 / d
244 *
245 * Compute columns K and K+1 of the inverse.
246 *
247  IF( k.GT.1 ) THEN
248  CALL dcopy( k-1, a( 1, k ), 1, work, 1 )
249  CALL dsymv( uplo, k-1, -one, a, lda, work, 1, zero,
250  $ a( 1, k ), 1 )
251  a( k, k ) = a( k, k ) - ddot( k-1, work, 1, a( 1, k ),
252  $ 1 )
253  a( k, k+1 ) = a( k, k+1 ) -
254  $ ddot( k-1, a( 1, k ), 1, a( 1, k+1 ), 1 )
255  CALL dcopy( k-1, a( 1, k+1 ), 1, work, 1 )
256  CALL dsymv( uplo, k-1, -one, a, lda, work, 1, zero,
257  $ a( 1, k+1 ), 1 )
258  a( k+1, k+1 ) = a( k+1, k+1 ) -
259  $ ddot( k-1, work, 1, a( 1, k+1 ), 1 )
260  END IF
261  kstep = 2
262  END IF
263 *
264  kp = abs( ipiv( k ) )
265  IF( kp.NE.k ) THEN
266 *
267 * Interchange rows and columns K and KP in the leading
268 * submatrix A(1:k+1,1:k+1)
269 *
270  CALL dswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
271  CALL dswap( k-kp-1, a( kp+1, k ), 1, a( kp, kp+1 ), lda )
272  temp = a( k, k )
273  a( k, k ) = a( kp, kp )
274  a( kp, kp ) = temp
275  IF( kstep.EQ.2 ) THEN
276  temp = a( k, k+1 )
277  a( k, k+1 ) = a( kp, k+1 )
278  a( kp, k+1 ) = temp
279  END IF
280  END IF
281 *
282  k = k + kstep
283  go to 30
284  40 continue
285 *
286  ELSE
287 *
288 * Compute inv(A) from the factorization A = L*D*L**T.
289 *
290 * K is the main loop index, increasing from 1 to N in steps of
291 * 1 or 2, depending on the size of the diagonal blocks.
292 *
293  k = n
294  50 continue
295 *
296 * If K < 1, exit from loop.
297 *
298  IF( k.LT.1 )
299  $ go to 60
300 *
301  IF( ipiv( k ).GT.0 ) THEN
302 *
303 * 1 x 1 diagonal block
304 *
305 * Invert the diagonal block.
306 *
307  a( k, k ) = one / a( k, k )
308 *
309 * Compute column K of the inverse.
310 *
311  IF( k.LT.n ) THEN
312  CALL dcopy( n-k, a( k+1, k ), 1, work, 1 )
313  CALL dsymv( uplo, n-k, -one, a( k+1, k+1 ), lda, work, 1,
314  $ zero, a( k+1, k ), 1 )
315  a( k, k ) = a( k, k ) - ddot( n-k, work, 1, a( k+1, k ),
316  $ 1 )
317  END IF
318  kstep = 1
319  ELSE
320 *
321 * 2 x 2 diagonal block
322 *
323 * Invert the diagonal block.
324 *
325  t = abs( a( k, k-1 ) )
326  ak = a( k-1, k-1 ) / t
327  akp1 = a( k, k ) / t
328  akkp1 = a( k, k-1 ) / t
329  d = t*( ak*akp1-one )
330  a( k-1, k-1 ) = akp1 / d
331  a( k, k ) = ak / d
332  a( k, k-1 ) = -akkp1 / d
333 *
334 * Compute columns K-1 and K of the inverse.
335 *
336  IF( k.LT.n ) THEN
337  CALL dcopy( n-k, a( k+1, k ), 1, work, 1 )
338  CALL dsymv( uplo, n-k, -one, a( k+1, k+1 ), lda, work, 1,
339  $ zero, a( k+1, k ), 1 )
340  a( k, k ) = a( k, k ) - ddot( n-k, work, 1, a( k+1, k ),
341  $ 1 )
342  a( k, k-1 ) = a( k, k-1 ) -
343  $ ddot( n-k, a( k+1, k ), 1, a( k+1, k-1 ),
344  $ 1 )
345  CALL dcopy( n-k, a( k+1, k-1 ), 1, work, 1 )
346  CALL dsymv( uplo, n-k, -one, a( k+1, k+1 ), lda, work, 1,
347  $ zero, a( k+1, k-1 ), 1 )
348  a( k-1, k-1 ) = a( k-1, k-1 ) -
349  $ ddot( n-k, work, 1, a( k+1, k-1 ), 1 )
350  END IF
351  kstep = 2
352  END IF
353 *
354  kp = abs( ipiv( k ) )
355  IF( kp.NE.k ) THEN
356 *
357 * Interchange rows and columns K and KP in the trailing
358 * submatrix A(k-1:n,k-1:n)
359 *
360  IF( kp.LT.n )
361  $ CALL dswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
362  CALL dswap( kp-k-1, a( k+1, k ), 1, a( kp, k+1 ), lda )
363  temp = a( k, k )
364  a( k, k ) = a( kp, kp )
365  a( kp, kp ) = temp
366  IF( kstep.EQ.2 ) THEN
367  temp = a( k, k-1 )
368  a( k, k-1 ) = a( kp, k-1 )
369  a( kp, k-1 ) = temp
370  END IF
371  END IF
372 *
373  k = k - kstep
374  go to 50
375  60 continue
376  END IF
377 *
378  return
379 *
380 * End of DSYTRI
381 *
382  END