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