LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ chetri()

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

CHETRI

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

Purpose:
 CHETRI computes the inverse of a complex Hermitian indefinite matrix
 A using the factorization A = U*D*U**H or A = L*D*L**H computed by
 CHETRF.
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**H;
          = 'L':  Lower triangular, form is A = L*D*L**H.
[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 CHETRF.

          On exit, if INFO = 0, the (Hermitian) 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 CHETRF.
[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
December 2016

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