LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ csytri_rook()

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