LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine csytri_rook ( character  UPLO,
integer  N,
complex, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
complex, dimension( * )  WORK,
integer  INFO 
)

CSYTRI_ROOK

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

Purpose: 
 CSYTRI_ROOK computes the inverse of a complex symmetric
 matrix A using the factorization A = U*D*U**T or A = L*D*L**T
 computed by CSYTRF_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 COMPLEX 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 CSYTRF_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 CSYTRF_ROOK.
[out]WORK
          WORK is COMPLEX 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.
Date
November 2015
Contributors: 
   November 2015, 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 131 of file csytri_rook.f.

131 *
132 * -- LAPACK computational routine (version 3.6.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 * November 2015
136 *
137 * .. Scalar Arguments ..
138  CHARACTER uplo
139  INTEGER info, lda, n
140 * ..
141 * .. Array Arguments ..
142  INTEGER ipiv( * )
143  COMPLEX a( lda, * ), work( * )
144 * ..
145 *
146 * =====================================================================
147 *
148 * .. Parameters ..
149  COMPLEX cone, czero
150  parameter ( cone = ( 1.0e+0, 0.0e+0 ),
151  $ czero = ( 0.0e+0, 0.0e+0 ) )
152 * ..
153 * .. Local Scalars ..
154  LOGICAL upper
155  INTEGER k, kp, kstep
156  COMPLEX ak, akkp1, akp1, d, t, temp
157 * ..
158 * .. External Functions ..
159  LOGICAL lsame
160  COMPLEX cdotu
161  EXTERNAL lsame, cdotu
162 * ..
163 * .. External Subroutines ..
164  EXTERNAL ccopy, cswap, csymv, xerbla
165 * ..
166 * .. Intrinsic Functions ..
167  INTRINSIC max
168 * ..
169 * .. Executable Statements ..
170 *
171 * Test the input parameters.
172 *
173  info = 0
174  upper = lsame( uplo, 'U' )
175  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
176  info = -1
177  ELSE IF( n.LT.0 ) THEN
178  info = -2
179  ELSE IF( lda.LT.max( 1, n ) ) THEN
180  info = -4
181  END IF
182  IF( info.NE.0 ) THEN
183  CALL xerbla( 'CSYTRI_ROOK', -info )
184  RETURN
185  END IF
186 *
187 * Quick return if possible
188 *
189  IF( n.EQ.0 )
190  $ RETURN
191 *
192 * Check that the diagonal matrix D is nonsingular.
193 *
194  IF( upper ) THEN
195 *
196 * Upper triangular storage: examine D from bottom to top
197 *
198  DO 10 info = n, 1, -1
199  IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.czero )
200  $ RETURN
201  10 CONTINUE
202  ELSE
203 *
204 * Lower triangular storage: examine D from top to bottom.
205 *
206  DO 20 info = 1, n
207  IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.czero )
208  $ RETURN
209  20 CONTINUE
210  END IF
211  info = 0
212 *
213  IF( upper ) THEN
214 *
215 * Compute inv(A) from the factorization A = U*D*U**T.
216 *
217 * K is the main loop index, increasing from 1 to N in steps of
218 * 1 or 2, depending on the size of the diagonal blocks.
219 *
220  k = 1
221  30 CONTINUE
222 *
223 * If K > N, exit from loop.
224 *
225  IF( k.GT.n )
226  $ GO TO 40
227 *
228  IF( ipiv( k ).GT.0 ) THEN
229 *
230 * 1 x 1 diagonal block
231 *
232 * Invert the diagonal block.
233 *
234  a( k, k ) = cone / a( k, k )
235 *
236 * Compute column K of the inverse.
237 *
238  IF( k.GT.1 ) THEN
239  CALL ccopy( k-1, a( 1, k ), 1, work, 1 )
240  CALL csymv( uplo, k-1, -cone, a, lda, work, 1, czero,
241  $ a( 1, k ), 1 )
242  a( k, k ) = a( k, k ) - cdotu( k-1, work, 1, a( 1, k ),
243  $ 1 )
244  END IF
245  kstep = 1
246  ELSE
247 *
248 * 2 x 2 diagonal block
249 *
250 * Invert the diagonal block.
251 *
252  t = a( k, k+1 )
253  ak = a( k, k ) / t
254  akp1 = a( k+1, k+1 ) / t
255  akkp1 = a( k, k+1 ) / t
256  d = t*( ak*akp1-cone )
257  a( k, k ) = akp1 / d
258  a( k+1, k+1 ) = ak / d
259  a( k, k+1 ) = -akkp1 / d
260 *
261 * Compute columns K and K+1 of the inverse.
262 *
263  IF( k.GT.1 ) THEN
264  CALL ccopy( k-1, a( 1, k ), 1, work, 1 )
265  CALL csymv( uplo, k-1, -cone, a, lda, work, 1, czero,
266  $ a( 1, k ), 1 )
267  a( k, k ) = a( k, k ) - cdotu( k-1, work, 1, a( 1, k ),
268  $ 1 )
269  a( k, k+1 ) = a( k, k+1 ) -
270  $ cdotu( k-1, a( 1, k ), 1, a( 1, k+1 ), 1 )
271  CALL ccopy( k-1, a( 1, k+1 ), 1, work, 1 )
272  CALL csymv( uplo, k-1, -cone, a, lda, work, 1, czero,
273  $ a( 1, k+1 ), 1 )
274  a( k+1, k+1 ) = a( k+1, k+1 ) -
275  $ cdotu( k-1, work, 1, a( 1, k+1 ), 1 )
276  END IF
277  kstep = 2
278  END IF
279 *
280  IF( kstep.EQ.1 ) THEN
281 *
282 * Interchange rows and columns K and IPIV(K) in the leading
283 * submatrix A(1:k+1,1:k+1)
284 *
285  kp = ipiv( k )
286  IF( kp.NE.k ) THEN
287  IF( kp.GT.1 )
288  $ CALL cswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
289  CALL cswap( k-kp-1, a( kp+1, k ), 1, a( kp, kp+1 ), lda )
290  temp = a( k, k )
291  a( k, k ) = a( kp, kp )
292  a( kp, kp ) = temp
293  END IF
294  ELSE
295 *
296 * Interchange rows and columns K and K+1 with -IPIV(K) and
297 * -IPIV(K+1)in the leading submatrix A(1:k+1,1:k+1)
298 *
299  kp = -ipiv( k )
300  IF( kp.NE.k ) THEN
301  IF( kp.GT.1 )
302  $ CALL cswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
303  CALL cswap( k-kp-1, a( kp+1, k ), 1, a( kp, kp+1 ), lda )
304 *
305  temp = a( k, k )
306  a( k, k ) = a( kp, kp )
307  a( kp, kp ) = temp
308  temp = a( k, k+1 )
309  a( k, k+1 ) = a( kp, k+1 )
310  a( kp, k+1 ) = temp
311  END IF
312 *
313  k = k + 1
314  kp = -ipiv( k )
315  IF( kp.NE.k ) THEN
316  IF( kp.GT.1 )
317  $ CALL cswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
318  CALL cswap( k-kp-1, a( kp+1, k ), 1, a( kp, kp+1 ), lda )
319  temp = a( k, k )
320  a( k, k ) = a( kp, kp )
321  a( kp, kp ) = temp
322  END IF
323  END IF
324 *
325  k = k + 1
326  GO TO 30
327  40 CONTINUE
328 *
329  ELSE
330 *
331 * Compute inv(A) from the factorization A = L*D*L**T.
332 *
333 * K is the main loop index, increasing from 1 to N in steps of
334 * 1 or 2, depending on the size of the diagonal blocks.
335 *
336  k = n
337  50 CONTINUE
338 *
339 * If K < 1, exit from loop.
340 *
341  IF( k.LT.1 )
342  $ GO TO 60
343 *
344  IF( ipiv( k ).GT.0 ) THEN
345 *
346 * 1 x 1 diagonal block
347 *
348 * Invert the diagonal block.
349 *
350  a( k, k ) = cone / a( k, k )
351 *
352 * Compute column K of the inverse.
353 *
354  IF( k.LT.n ) THEN
355  CALL ccopy( n-k, a( k+1, k ), 1, work, 1 )
356  CALL csymv( uplo, n-k,-cone, a( k+1, k+1 ), lda, work, 1,
357  $ czero, a( k+1, k ), 1 )
358  a( k, k ) = a( k, k ) - cdotu( n-k, work, 1, a( k+1, k ),
359  $ 1 )
360  END IF
361  kstep = 1
362  ELSE
363 *
364 * 2 x 2 diagonal block
365 *
366 * Invert the diagonal block.
367 *
368  t = a( k, k-1 )
369  ak = a( k-1, k-1 ) / t
370  akp1 = a( k, k ) / t
371  akkp1 = a( k, k-1 ) / t
372  d = t*( ak*akp1-cone )
373  a( k-1, k-1 ) = akp1 / d
374  a( k, k ) = ak / d
375  a( k, k-1 ) = -akkp1 / d
376 *
377 * Compute columns K-1 and K of the inverse.
378 *
379  IF( k.LT.n ) THEN
380  CALL ccopy( n-k, a( k+1, k ), 1, work, 1 )
381  CALL csymv( uplo, n-k,-cone, a( k+1, k+1 ), lda, work, 1,
382  $ czero, a( k+1, k ), 1 )
383  a( k, k ) = a( k, k ) - cdotu( n-k, work, 1, a( k+1, k ),
384  $ 1 )
385  a( k, k-1 ) = a( k, k-1 ) -
386  $ cdotu( n-k, a( k+1, k ), 1, a( k+1, k-1 ),
387  $ 1 )
388  CALL ccopy( n-k, a( k+1, k-1 ), 1, work, 1 )
389  CALL csymv( uplo, n-k,-cone, a( k+1, k+1 ), lda, work, 1,
390  $ czero, a( k+1, k-1 ), 1 )
391  a( k-1, k-1 ) = a( k-1, k-1 ) -
392  $ cdotu( n-k, work, 1, a( k+1, k-1 ), 1 )
393  END IF
394  kstep = 2
395  END IF
396 *
397  IF( kstep.EQ.1 ) THEN
398 *
399 * Interchange rows and columns K and IPIV(K) in the trailing
400 * submatrix A(k-1:n,k-1:n)
401 *
402  kp = ipiv( k )
403  IF( kp.NE.k ) THEN
404  IF( kp.LT.n )
405  $ CALL cswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
406  CALL cswap( kp-k-1, a( k+1, k ), 1, a( kp, k+1 ), lda )
407  temp = a( k, k )
408  a( k, k ) = a( kp, kp )
409  a( kp, kp ) = temp
410  END IF
411  ELSE
412 *
413 * Interchange rows and columns K and K-1 with -IPIV(K) and
414 * -IPIV(K-1) in the trailing submatrix A(k-1:n,k-1:n)
415 *
416  kp = -ipiv( k )
417  IF( kp.NE.k ) THEN
418  IF( kp.LT.n )
419  $ CALL cswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
420  CALL cswap( kp-k-1, a( k+1, k ), 1, a( kp, k+1 ), lda )
421 *
422  temp = a( k, k )
423  a( k, k ) = a( kp, kp )
424  a( kp, kp ) = temp
425  temp = a( k, k-1 )
426  a( k, k-1 ) = a( kp, k-1 )
427  a( kp, k-1 ) = temp
428  END IF
429 *
430  k = k - 1
431  kp = -ipiv( k )
432  IF( kp.NE.k ) THEN
433  IF( kp.LT.n )
434  $ CALL cswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
435  CALL cswap( kp-k-1, a( k+1, k ), 1, a( kp, k+1 ), lda )
436  temp = a( k, k )
437  a( k, k ) = a( kp, kp )
438  a( kp, kp ) = temp
439  END IF
440  END IF
441 *
442  k = k - 1
443  GO TO 50
444  60 CONTINUE
445  END IF
446 *
447  RETURN
448 *
449 * End of CSYTRI_ROOK
450 *
csymv
subroutine csymv(UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CSYMV computes a matrix-vector product for a complex symmetric matrix.
Definition: csymv.f:159
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
cdotu
complex function cdotu(N, CX, INCX, CY, INCY)
CDOTU
Definition: cdotu.f:54
ccopy
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:52
cswap
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:52
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function: