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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
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*
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: