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

## ◆ chetri_rook()

 subroutine chetri_rook ( character uplo, integer n, complex, dimension( lda, * ) a, integer lda, integer, dimension( * ) ipiv, complex, dimension( * ) work, integer info )

CHETRI_ROOK computes the inverse of HE matrix using the factorization obtained with the bounded Bunch-Kaufman ("rook") diagonal pivoting method.

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

Purpose:
``` CHETRI_ROOK 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_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**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_ROOK. 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_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.```
Contributors:
```  November 2013,  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 127 of file chetri_rook.f.

128*
129* -- LAPACK computational routine --
130* -- LAPACK is a software package provided by Univ. of Tennessee, --
131* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
132*
133* .. Scalar Arguments ..
134 CHARACTER UPLO
135 INTEGER INFO, LDA, N
136* ..
137* .. Array Arguments ..
138 INTEGER IPIV( * )
139 COMPLEX A( LDA, * ), WORK( * )
140* ..
141*
142* =====================================================================
143*
144* .. Parameters ..
145 REAL ONE
146 COMPLEX CONE, CZERO
147 parameter( one = 1.0e+0, cone = ( 1.0e+0, 0.0e+0 ),
148 \$ czero = ( 0.0e+0, 0.0e+0 ) )
149* ..
150* .. Local Scalars ..
151 LOGICAL UPPER
152 INTEGER J, K, KP, KSTEP
153 REAL AK, AKP1, D, T
154 COMPLEX AKKP1, TEMP
155* ..
156* .. External Functions ..
157 LOGICAL LSAME
158 COMPLEX CDOTC
159 EXTERNAL lsame, cdotc
160* ..
161* .. External Subroutines ..
162 EXTERNAL ccopy, chemv, cswap, xerbla
163* ..
164* .. Intrinsic Functions ..
165 INTRINSIC abs, conjg, max, real
166* ..
167* .. Executable Statements ..
168*
169* Test the input parameters.
170*
171 info = 0
172 upper = lsame( uplo, 'U' )
173 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
174 info = -1
175 ELSE IF( n.LT.0 ) THEN
176 info = -2
177 ELSE IF( lda.LT.max( 1, n ) ) THEN
178 info = -4
179 END IF
180 IF( info.NE.0 ) THEN
181 CALL xerbla( 'CHETRI_ROOK', -info )
182 RETURN
183 END IF
184*
185* Quick return if possible
186*
187 IF( n.EQ.0 )
188 \$ RETURN
189*
190* Check that the diagonal matrix D is nonsingular.
191*
192 IF( upper ) THEN
193*
194* Upper triangular storage: examine D from bottom to top
195*
196 DO 10 info = n, 1, -1
197 IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.czero )
198 \$ RETURN
199 10 CONTINUE
200 ELSE
201*
202* Lower triangular storage: examine D from top to bottom.
203*
204 DO 20 info = 1, n
205 IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.czero )
206 \$ RETURN
207 20 CONTINUE
208 END IF
209 info = 0
210*
211 IF( upper ) THEN
212*
213* Compute inv(A) from the factorization A = U*D*U**H.
214*
215* K is the main loop index, increasing from 1 to N in steps of
216* 1 or 2, depending on the size of the diagonal blocks.
217*
218 k = 1
219 30 CONTINUE
220*
221* If K > N, exit from loop.
222*
223 IF( k.GT.n )
224 \$ GO TO 70
225*
226 IF( ipiv( k ).GT.0 ) THEN
227*
228* 1 x 1 diagonal block
229*
230* Invert the diagonal block.
231*
232 a( k, k ) = one / real( a( k, k ) )
233*
234* Compute column K of the inverse.
235*
236 IF( k.GT.1 ) THEN
237 CALL ccopy( k-1, a( 1, k ), 1, work, 1 )
238 CALL chemv( uplo, k-1, -cone, a, lda, work, 1, czero,
239 \$ a( 1, k ), 1 )
240 a( k, k ) = a( k, k ) - real( cdotc( k-1, work, 1, a( 1,
241 \$ k ), 1 ) )
242 END IF
243 kstep = 1
244 ELSE
245*
246* 2 x 2 diagonal block
247*
248* Invert the diagonal block.
249*
250 t = abs( a( k, k+1 ) )
251 ak = real( a( k, k ) ) / t
252 akp1 = real( a( k+1, k+1 ) ) / t
253 akkp1 = a( k, k+1 ) / t
254 d = t*( ak*akp1-one )
255 a( k, k ) = akp1 / d
256 a( k+1, k+1 ) = ak / d
257 a( k, k+1 ) = -akkp1 / d
258*
259* Compute columns K and K+1 of the inverse.
260*
261 IF( k.GT.1 ) THEN
262 CALL ccopy( k-1, a( 1, k ), 1, work, 1 )
263 CALL chemv( uplo, k-1, -cone, a, lda, work, 1, czero,
264 \$ a( 1, k ), 1 )
265 a( k, k ) = a( k, k ) - real( cdotc( k-1, work, 1, a( 1,
266 \$ k ), 1 ) )
267 a( k, k+1 ) = a( k, k+1 ) -
268 \$ cdotc( k-1, a( 1, k ), 1, a( 1, k+1 ), 1 )
269 CALL ccopy( k-1, a( 1, k+1 ), 1, work, 1 )
270 CALL chemv( uplo, k-1, -cone, a, lda, work, 1, czero,
271 \$ a( 1, k+1 ), 1 )
272 a( k+1, k+1 ) = a( k+1, k+1 ) -
273 \$ real( cdotc( k-1, work, 1, a( 1, k+1 ),
274 \$ 1 ) )
275 END IF
276 kstep = 2
277 END IF
278*
279 IF( kstep.EQ.1 ) THEN
280*
281* Interchange rows and columns K and IPIV(K) in the leading
282* submatrix A(1:k,1:k)
283*
284 kp = ipiv( k )
285 IF( kp.NE.k ) THEN
286*
287 IF( kp.GT.1 )
288 \$ CALL cswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
289*
290 DO 40 j = kp + 1, k - 1
291 temp = conjg( a( j, k ) )
292 a( j, k ) = conjg( a( kp, j ) )
293 a( kp, j ) = temp
294 40 CONTINUE
295*
296 a( kp, k ) = conjg( a( kp, k ) )
297*
298 temp = a( k, k )
299 a( k, k ) = a( kp, kp )
300 a( kp, kp ) = temp
301 END IF
302 ELSE
303*
304* Interchange rows and columns K and K+1 with -IPIV(K) and
305* -IPIV(K+1) in the leading submatrix A(k+1:n,k+1:n)
306*
307* (1) Interchange rows and columns K and -IPIV(K)
308*
309 kp = -ipiv( k )
310 IF( kp.NE.k ) THEN
311*
312 IF( kp.GT.1 )
313 \$ CALL cswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
314*
315 DO 50 j = kp + 1, k - 1
316 temp = conjg( a( j, k ) )
317 a( j, k ) = conjg( a( kp, j ) )
318 a( kp, j ) = temp
319 50 CONTINUE
320*
321 a( kp, k ) = conjg( a( kp, k ) )
322*
323 temp = a( k, k )
324 a( k, k ) = a( kp, kp )
325 a( kp, kp ) = temp
326*
327 temp = a( k, k+1 )
328 a( k, k+1 ) = a( kp, k+1 )
329 a( kp, k+1 ) = temp
330 END IF
331*
332* (2) Interchange rows and columns K+1 and -IPIV(K+1)
333*
334 k = k + 1
335 kp = -ipiv( k )
336 IF( kp.NE.k ) THEN
337*
338 IF( kp.GT.1 )
339 \$ CALL cswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
340*
341 DO 60 j = kp + 1, k - 1
342 temp = conjg( a( j, k ) )
343 a( j, k ) = conjg( a( kp, j ) )
344 a( kp, j ) = temp
345 60 CONTINUE
346*
347 a( kp, k ) = conjg( a( kp, k ) )
348*
349 temp = a( k, k )
350 a( k, k ) = a( kp, kp )
351 a( kp, kp ) = temp
352 END IF
353 END IF
354*
355 k = k + 1
356 GO TO 30
357 70 CONTINUE
358*
359 ELSE
360*
361* Compute inv(A) from the factorization A = L*D*L**H.
362*
363* K is the main loop index, decreasing from N to 1 in steps of
364* 1 or 2, depending on the size of the diagonal blocks.
365*
366 k = n
367 80 CONTINUE
368*
369* If K < 1, exit from loop.
370*
371 IF( k.LT.1 )
372 \$ GO TO 120
373*
374 IF( ipiv( k ).GT.0 ) THEN
375*
376* 1 x 1 diagonal block
377*
378* Invert the diagonal block.
379*
380 a( k, k ) = one / real( a( k, k ) )
381*
382* Compute column K of the inverse.
383*
384 IF( k.LT.n ) THEN
385 CALL ccopy( n-k, a( k+1, k ), 1, work, 1 )
386 CALL chemv( uplo, n-k, -cone, a( k+1, k+1 ), lda, work,
387 \$ 1, czero, a( k+1, k ), 1 )
388 a( k, k ) = a( k, k ) - real( cdotc( n-k, work, 1,
389 \$ a( k+1, k ), 1 ) )
390 END IF
391 kstep = 1
392 ELSE
393*
394* 2 x 2 diagonal block
395*
396* Invert the diagonal block.
397*
398 t = abs( a( k, k-1 ) )
399 ak = real( a( k-1, k-1 ) ) / t
400 akp1 = real( a( k, k ) ) / t
401 akkp1 = a( k, k-1 ) / t
402 d = t*( ak*akp1-one )
403 a( k-1, k-1 ) = akp1 / d
404 a( k, k ) = ak / d
405 a( k, k-1 ) = -akkp1 / d
406*
407* Compute columns K-1 and K of the inverse.
408*
409 IF( k.LT.n ) THEN
410 CALL ccopy( n-k, a( k+1, k ), 1, work, 1 )
411 CALL chemv( uplo, n-k, -cone, a( k+1, k+1 ), lda, work,
412 \$ 1, czero, a( k+1, k ), 1 )
413 a( k, k ) = a( k, k ) - real( cdotc( n-k, work, 1,
414 \$ a( k+1, k ), 1 ) )
415 a( k, k-1 ) = a( k, k-1 ) -
416 \$ cdotc( n-k, a( k+1, k ), 1, a( k+1, k-1 ),
417 \$ 1 )
418 CALL ccopy( n-k, a( k+1, k-1 ), 1, work, 1 )
419 CALL chemv( uplo, n-k, -cone, a( k+1, k+1 ), lda, work,
420 \$ 1, czero, a( k+1, k-1 ), 1 )
421 a( k-1, k-1 ) = a( k-1, k-1 ) -
422 \$ real( cdotc( n-k, work, 1, a( k+1, k-1 ),
423 \$ 1 ) )
424 END IF
425 kstep = 2
426 END IF
427*
428 IF( kstep.EQ.1 ) THEN
429*
430* Interchange rows and columns K and IPIV(K) in the trailing
431* submatrix A(k:n,k:n)
432*
433 kp = ipiv( k )
434 IF( kp.NE.k ) THEN
435*
436 IF( kp.LT.n )
437 \$ CALL cswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
438*
439 DO 90 j = k + 1, kp - 1
440 temp = conjg( a( j, k ) )
441 a( j, k ) = conjg( a( kp, j ) )
442 a( kp, j ) = temp
443 90 CONTINUE
444*
445 a( kp, k ) = conjg( a( kp, k ) )
446*
447 temp = a( k, k )
448 a( k, k ) = a( kp, kp )
449 a( kp, kp ) = temp
450 END IF
451 ELSE
452*
453* Interchange rows and columns K and K-1 with -IPIV(K) and
454* -IPIV(K-1) in the trailing submatrix A(k-1:n,k-1:n)
455*
456* (1) Interchange rows and columns K and -IPIV(K)
457*
458 kp = -ipiv( k )
459 IF( kp.NE.k ) THEN
460*
461 IF( kp.LT.n )
462 \$ CALL cswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
463*
464 DO 100 j = k + 1, kp - 1
465 temp = conjg( a( j, k ) )
466 a( j, k ) = conjg( a( kp, j ) )
467 a( kp, j ) = temp
468 100 CONTINUE
469*
470 a( kp, k ) = conjg( a( kp, k ) )
471*
472 temp = a( k, k )
473 a( k, k ) = a( kp, kp )
474 a( kp, kp ) = temp
475*
476 temp = a( k, k-1 )
477 a( k, k-1 ) = a( kp, k-1 )
478 a( kp, k-1 ) = temp
479 END IF
480*
481* (2) Interchange rows and columns K-1 and -IPIV(K-1)
482*
483 k = k - 1
484 kp = -ipiv( k )
485 IF( kp.NE.k ) THEN
486*
487 IF( kp.LT.n )
488 \$ CALL cswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
489*
490 DO 110 j = k + 1, kp - 1
491 temp = conjg( a( j, k ) )
492 a( j, k ) = conjg( a( kp, j ) )
493 a( kp, j ) = temp
494 110 CONTINUE
495*
496 a( kp, k ) = conjg( a( kp, k ) )
497*
498 temp = a( k, k )
499 a( k, k ) = a( kp, kp )
500 a( kp, kp ) = temp
501 END IF
502 END IF
503*
504 k = k - 1
505 GO TO 80
506 120 CONTINUE
507 END IF
508*
509 RETURN
510*
511* End of CHETRI_ROOK
512*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
complex function cdotc(n, cx, incx, cy, incy)
CDOTC
Definition cdotc.f:83
subroutine chemv(uplo, n, alpha, a, lda, x, incx, beta, y, incy)
CHEMV
Definition chemv.f:154
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
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: