LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ ssytri_rook()

subroutine ssytri_rook ( character  uplo,
integer  n,
real, dimension( lda, * )  a,
integer  lda,
integer, dimension( * )  ipiv,
real, dimension( * )  work,
integer  info 
)

SSYTRI_ROOK

Download SSYTRI_ROOK + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 SSYTRI_ROOK computes the inverse of a real symmetric
 matrix A using the factorization A = U*D*U**T or A = L*D*L**T
 computed by SSYTRF_ROOK.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          Specifies whether the details of the factorization are stored
          as an upper or lower triangular matrix.
          = 'U':  Upper triangular, form is A = U*D*U**T;
          = 'L':  Lower triangular, form is A = L*D*L**T.
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in,out]A
          A is REAL array, dimension (LDA,N)
          On entry, the block diagonal matrix D and the multipliers
          used to obtain the factor U or L as computed by SSYTRF_ROOK.

          On exit, if INFO = 0, the (symmetric) inverse of the original
          matrix.  If UPLO = 'U', the upper triangular part of the
          inverse is formed and the part of A below the diagonal is not
          referenced; if UPLO = 'L' the lower triangular part of the
          inverse is formed and the part of A above the diagonal is
          not referenced.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[in]IPIV
          IPIV is INTEGER array, dimension (N)
          Details of the interchanges and the block structure of D
          as determined by SSYTRF_ROOK.
[out]WORK
          WORK is REAL array, dimension (N)
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value
          > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
               inverse could not be computed.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
   April 2012, Igor Kozachenko,
                  Computer Science Division,
                  University of California, Berkeley

  September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
                  School of Mathematics,
                  University of Manchester

Definition at line 128 of file ssytri_rook.f.

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 REAL A( LDA, * ), WORK( * )
141* ..
142*
143* =====================================================================
144*
145* .. Parameters ..
146 REAL ONE, ZERO
147 parameter( one = 1.0e+0, zero = 0.0e+0 )
148* ..
149* .. Local Scalars ..
150 LOGICAL UPPER
151 INTEGER K, KP, KSTEP
152 REAL AK, AKKP1, AKP1, D, T, TEMP
153* ..
154* .. External Functions ..
155 LOGICAL LSAME
156 REAL SDOT
157 EXTERNAL lsame, sdot
158* ..
159* .. External Subroutines ..
160 EXTERNAL scopy, sswap, ssymv, 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( 'SSYTRI_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 scopy( k-1, a( 1, k ), 1, work, 1 )
236 CALL ssymv( uplo, k-1, -one, a, lda, work, 1, zero,
237 $ a( 1, k ), 1 )
238 a( k, k ) = a( k, k ) - sdot( 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 scopy( k-1, a( 1, k ), 1, work, 1 )
261 CALL ssymv( uplo, k-1, -one, a, lda, work, 1, zero,
262 $ a( 1, k ), 1 )
263 a( k, k ) = a( k, k ) - sdot( k-1, work, 1, a( 1, k ),
264 $ 1 )
265 a( k, k+1 ) = a( k, k+1 ) -
266 $ sdot( k-1, a( 1, k ), 1, a( 1, k+1 ), 1 )
267 CALL scopy( k-1, a( 1, k+1 ), 1, work, 1 )
268 CALL ssymv( 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 $ sdot( 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 sswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
285 CALL sswap( 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 sswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
299 CALL sswap( 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 sswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
314 CALL sswap( 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 scopy( n-k, a( k+1, k ), 1, work, 1 )
352 CALL ssymv( 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 ) - sdot( 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 scopy( n-k, a( k+1, k ), 1, work, 1 )
377 CALL ssymv( 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 ) - sdot( n-k, work, 1, a( k+1, k ),
380 $ 1 )
381 a( k, k-1 ) = a( k, k-1 ) -
382 $ sdot( n-k, a( k+1, k ), 1, a( k+1, k-1 ),
383 $ 1 )
384 CALL scopy( n-k, a( k+1, k-1 ), 1, work, 1 )
385 CALL ssymv( 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 $ sdot( 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 sswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
402 CALL sswap( 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 sswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
416 CALL sswap( 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 sswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
431 CALL sswap( 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 SSYTRI_ROOK
446*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
sdot
real function sdot(n, sx, incx, sy, incy)
SDOT
Definition sdot.f:82
ssymv
subroutine ssymv(uplo, n, alpha, a, lda, x, incx, beta, y, incy)
SSYMV
Definition ssymv.f:152
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
sswap
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82
Here is the call graph for this function:
Here is the caller graph for this function: