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

◆ 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.

Definition at line 111 of file chetri.f.

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