LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ csytri()

subroutine csytri ( character  UPLO,
integer  N,
complex, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
complex, dimension( * )  WORK,
integer  INFO 
)

CSYTRI

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

Purpose:
 CSYTRI computes the inverse of a complex symmetric indefinite matrix
 A using the factorization A = U*D*U**T or A = L*D*L**T computed by
 CSYTRF.
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.

          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.
[out]WORK
          WORK is COMPLEX array, dimension (2*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
December 2016

Definition at line 116 of file csytri.f.

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