LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ 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.
Date
November 2013
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 130 of file zhetri_rook.f.

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