LAPACK  3.10.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.
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 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:81
complex *16 function zdotc(N, ZX, INCX, ZY, INCY)
ZDOTC
Definition: zdotc.f:83
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:81
subroutine zhemv(UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZHEMV
Definition: zhemv.f:154
Here is the call graph for this function:
Here is the caller graph for this function: