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

◆ zhetri_rook()

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

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

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

Purpose:
 ZHETRI_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
 ZHETRF_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*16 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 ZHETRF_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 ZHETRF_ROOK.
[out]WORK
          WORK is COMPLEX*16 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 zhetri_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*16 A( LDA, * ), WORK( * )
140* ..
141*
142* =====================================================================
143*
144* .. Parameters ..
145 DOUBLE PRECISION ONE
146 COMPLEX*16 CONE, CZERO
147 parameter( one = 1.0d+0, cone = ( 1.0d+0, 0.0d+0 ),
148 $ czero = ( 0.0d+0, 0.0d+0 ) )
149* ..
150* .. Local Scalars ..
151 LOGICAL UPPER
152 INTEGER J, K, KP, KSTEP
153 DOUBLE PRECISION AK, AKP1, D, T
154 COMPLEX*16 AKKP1, TEMP
155* ..
156* .. External Functions ..
157 LOGICAL LSAME
158 COMPLEX*16 ZDOTC
159 EXTERNAL lsame, zdotc
160* ..
161* .. External Subroutines ..
162 EXTERNAL zcopy, zhemv, zswap, xerbla
163* ..
164* .. Intrinsic Functions ..
165 INTRINSIC abs, dconjg, max, dble
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( 'ZHETRI_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 / dble( a( k, k ) )
233*
234* Compute column K of the inverse.
235*
236 IF( k.GT.1 ) THEN
237 CALL zcopy( k-1, a( 1, k ), 1, work, 1 )
238 CALL zhemv( uplo, k-1, -cone, a, lda, work, 1, czero,
239 $ a( 1, k ), 1 )
240 a( k, k ) = a( k, k ) - dble( zdotc( 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 = dble( a( k, k ) ) / t
252 akp1 = dble( 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 zcopy( k-1, a( 1, k ), 1, work, 1 )
263 CALL zhemv( uplo, k-1, -cone, a, lda, work, 1, czero,
264 $ a( 1, k ), 1 )
265 a( k, k ) = a( k, k ) - dble( zdotc( k-1, work, 1, a( 1,
266 $ k ), 1 ) )
267 a( k, k+1 ) = a( k, k+1 ) -
268 $ zdotc( k-1, a( 1, k ), 1, a( 1, k+1 ), 1 )
269 CALL zcopy( k-1, a( 1, k+1 ), 1, work, 1 )
270 CALL zhemv( 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 $ dble( zdotc( 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 zswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
289*
290 DO 40 j = kp + 1, k - 1
291 temp = dconjg( a( j, k ) )
292 a( j, k ) = dconjg( a( kp, j ) )
293 a( kp, j ) = temp
294 40 CONTINUE
295*
296 a( kp, k ) = dconjg( 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 zswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
314*
315 DO 50 j = kp + 1, k - 1
316 temp = dconjg( a( j, k ) )
317 a( j, k ) = dconjg( a( kp, j ) )
318 a( kp, j ) = temp
319 50 CONTINUE
320*
321 a( kp, k ) = dconjg( 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 zswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
340*
341 DO 60 j = kp + 1, k - 1
342 temp = dconjg( a( j, k ) )
343 a( j, k ) = dconjg( a( kp, j ) )
344 a( kp, j ) = temp
345 60 CONTINUE
346*
347 a( kp, k ) = dconjg( 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 / dble( a( k, k ) )
381*
382* Compute column K of the inverse.
383*
384 IF( k.LT.n ) THEN
385 CALL zcopy( n-k, a( k+1, k ), 1, work, 1 )
386 CALL zhemv( 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 ) - dble( zdotc( 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 = dble( a( k-1, k-1 ) ) / t
400 akp1 = dble( 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 zcopy( n-k, a( k+1, k ), 1, work, 1 )
411 CALL zhemv( 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 ) - dble( zdotc( n-k, work, 1,
414 $ a( k+1, k ), 1 ) )
415 a( k, k-1 ) = a( k, k-1 ) -
416 $ zdotc( n-k, a( k+1, k ), 1, a( k+1, k-1 ),
417 $ 1 )
418 CALL zcopy( n-k, a( k+1, k-1 ), 1, work, 1 )
419 CALL zhemv( 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 $ dble( zdotc( 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 zswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
438*
439 DO 90 j = k + 1, kp - 1
440 temp = dconjg( a( j, k ) )
441 a( j, k ) = dconjg( a( kp, j ) )
442 a( kp, j ) = temp
443 90 CONTINUE
444*
445 a( kp, k ) = dconjg( 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 zswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
463*
464 DO 100 j = k + 1, kp - 1
465 temp = dconjg( a( j, k ) )
466 a( j, k ) = dconjg( a( kp, j ) )
467 a( kp, j ) = temp
468 100 CONTINUE
469*
470 a( kp, k ) = dconjg( 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 zswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
489*
490 DO 110 j = k + 1, kp - 1
491 temp = dconjg( a( j, k ) )
492 a( j, k ) = dconjg( a( kp, j ) )
493 a( kp, j ) = temp
494 110 CONTINUE
495*
496 a( kp, k ) = dconjg( 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 ZHETRI_ROOK
512*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zcopy
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
zdotc
complex *16 function zdotc(n, zx, incx, zy, incy)
ZDOTC
Definition zdotc.f:83
zhemv
subroutine zhemv(uplo, n, alpha, a, lda, x, incx, beta, y, incy)
ZHEMV
Definition zhemv.f:154
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
zswap
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP
Definition zswap.f:81
Here is the call graph for this function:
Here is the caller graph for this function: