LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ zhetri()

subroutine zhetri ( character  UPLO,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
complex*16, dimension( * )  WORK,
integer  INFO 
)

ZHETRI

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

Purpose:
 ZHETRI 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.
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.

          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.
[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.

Definition at line 113 of file zhetri.f.

114 *
115 * -- LAPACK computational routine --
116 * -- LAPACK is a software package provided by Univ. of Tennessee, --
117 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
118 *
119 * .. Scalar Arguments ..
120  CHARACTER UPLO
121  INTEGER INFO, LDA, N
122 * ..
123 * .. Array Arguments ..
124  INTEGER IPIV( * )
125  COMPLEX*16 A( LDA, * ), WORK( * )
126 * ..
127 *
128 * =====================================================================
129 *
130 * .. Parameters ..
131  DOUBLE PRECISION ONE
132  COMPLEX*16 CONE, ZERO
133  parameter( one = 1.0d+0, cone = ( 1.0d+0, 0.0d+0 ),
134  $ zero = ( 0.0d+0, 0.0d+0 ) )
135 * ..
136 * .. Local Scalars ..
137  LOGICAL UPPER
138  INTEGER J, K, KP, KSTEP
139  DOUBLE PRECISION AK, AKP1, D, T
140  COMPLEX*16 AKKP1, TEMP
141 * ..
142 * .. External Functions ..
143  LOGICAL LSAME
144  COMPLEX*16 ZDOTC
145  EXTERNAL lsame, zdotc
146 * ..
147 * .. External Subroutines ..
148  EXTERNAL xerbla, zcopy, zhemv, zswap
149 * ..
150 * .. Intrinsic Functions ..
151  INTRINSIC abs, dble, dconjg, max
152 * ..
153 * .. Executable Statements ..
154 *
155 * Test the input parameters.
156 *
157  info = 0
158  upper = lsame( uplo, 'U' )
159  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
160  info = -1
161  ELSE IF( n.LT.0 ) THEN
162  info = -2
163  ELSE IF( lda.LT.max( 1, n ) ) THEN
164  info = -4
165  END IF
166  IF( info.NE.0 ) THEN
167  CALL xerbla( 'ZHETRI', -info )
168  RETURN
169  END IF
170 *
171 * Quick return if possible
172 *
173  IF( n.EQ.0 )
174  $ RETURN
175 *
176 * Check that the diagonal matrix D is nonsingular.
177 *
178  IF( upper ) THEN
179 *
180 * Upper triangular storage: examine D from bottom to top
181 *
182  DO 10 info = n, 1, -1
183  IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.zero )
184  $ RETURN
185  10 CONTINUE
186  ELSE
187 *
188 * Lower triangular storage: examine D from top to bottom.
189 *
190  DO 20 info = 1, n
191  IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.zero )
192  $ RETURN
193  20 CONTINUE
194  END IF
195  info = 0
196 *
197  IF( upper ) THEN
198 *
199 * Compute inv(A) from the factorization A = U*D*U**H.
200 *
201 * K is the main loop index, increasing from 1 to N in steps of
202 * 1 or 2, depending on the size of the diagonal blocks.
203 *
204  k = 1
205  30 CONTINUE
206 *
207 * If K > N, exit from loop.
208 *
209  IF( k.GT.n )
210  $ GO TO 50
211 *
212  IF( ipiv( k ).GT.0 ) THEN
213 *
214 * 1 x 1 diagonal block
215 *
216 * Invert the diagonal block.
217 *
218  a( k, k ) = one / dble( a( k, k ) )
219 *
220 * Compute column K of the inverse.
221 *
222  IF( k.GT.1 ) THEN
223  CALL zcopy( k-1, a( 1, k ), 1, work, 1 )
224  CALL zhemv( uplo, k-1, -cone, a, lda, work, 1, zero,
225  $ a( 1, k ), 1 )
226  a( k, k ) = a( k, k ) - dble( zdotc( k-1, work, 1, a( 1,
227  $ k ), 1 ) )
228  END IF
229  kstep = 1
230  ELSE
231 *
232 * 2 x 2 diagonal block
233 *
234 * Invert the diagonal block.
235 *
236  t = abs( a( k, k+1 ) )
237  ak = dble( a( k, k ) ) / t
238  akp1 = dble( a( k+1, k+1 ) ) / t
239  akkp1 = a( k, k+1 ) / t
240  d = t*( ak*akp1-one )
241  a( k, k ) = akp1 / d
242  a( k+1, k+1 ) = ak / d
243  a( k, k+1 ) = -akkp1 / d
244 *
245 * Compute columns K and K+1 of the inverse.
246 *
247  IF( k.GT.1 ) THEN
248  CALL zcopy( k-1, a( 1, k ), 1, work, 1 )
249  CALL zhemv( uplo, k-1, -cone, a, lda, work, 1, zero,
250  $ a( 1, k ), 1 )
251  a( k, k ) = a( k, k ) - dble( zdotc( k-1, work, 1, a( 1,
252  $ k ), 1 ) )
253  a( k, k+1 ) = a( k, k+1 ) -
254  $ zdotc( k-1, a( 1, k ), 1, a( 1, k+1 ), 1 )
255  CALL zcopy( k-1, a( 1, k+1 ), 1, work, 1 )
256  CALL zhemv( uplo, k-1, -cone, a, lda, work, 1, zero,
257  $ a( 1, k+1 ), 1 )
258  a( k+1, k+1 ) = a( k+1, k+1 ) -
259  $ dble( zdotc( k-1, work, 1, a( 1, k+1 ),
260  $ 1 ) )
261  END IF
262  kstep = 2
263  END IF
264 *
265  kp = abs( ipiv( k ) )
266  IF( kp.NE.k ) THEN
267 *
268 * Interchange rows and columns K and KP in the leading
269 * submatrix A(1:k+1,1:k+1)
270 *
271  CALL zswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
272  DO 40 j = kp + 1, k - 1
273  temp = dconjg( a( j, k ) )
274  a( j, k ) = dconjg( a( kp, j ) )
275  a( kp, j ) = temp
276  40 CONTINUE
277  a( kp, k ) = dconjg( a( kp, k ) )
278  temp = a( k, k )
279  a( k, k ) = a( kp, kp )
280  a( kp, kp ) = temp
281  IF( kstep.EQ.2 ) THEN
282  temp = a( k, k+1 )
283  a( k, k+1 ) = a( kp, k+1 )
284  a( kp, k+1 ) = temp
285  END IF
286  END IF
287 *
288  k = k + kstep
289  GO TO 30
290  50 CONTINUE
291 *
292  ELSE
293 *
294 * Compute inv(A) from the factorization A = L*D*L**H.
295 *
296 * K is the main loop index, increasing from 1 to N in steps of
297 * 1 or 2, depending on the size of the diagonal blocks.
298 *
299  k = n
300  60 CONTINUE
301 *
302 * If K < 1, exit from loop.
303 *
304  IF( k.LT.1 )
305  $ GO TO 80
306 *
307  IF( ipiv( k ).GT.0 ) THEN
308 *
309 * 1 x 1 diagonal block
310 *
311 * Invert the diagonal block.
312 *
313  a( k, k ) = one / dble( a( k, k ) )
314 *
315 * Compute column K of the inverse.
316 *
317  IF( k.LT.n ) THEN
318  CALL zcopy( n-k, a( k+1, k ), 1, work, 1 )
319  CALL zhemv( uplo, n-k, -cone, a( k+1, k+1 ), lda, work,
320  $ 1, zero, a( k+1, k ), 1 )
321  a( k, k ) = a( k, k ) - dble( zdotc( n-k, work, 1,
322  $ a( k+1, k ), 1 ) )
323  END IF
324  kstep = 1
325  ELSE
326 *
327 * 2 x 2 diagonal block
328 *
329 * Invert the diagonal block.
330 *
331  t = abs( a( k, k-1 ) )
332  ak = dble( a( k-1, k-1 ) ) / t
333  akp1 = dble( a( k, k ) ) / t
334  akkp1 = a( k, k-1 ) / t
335  d = t*( ak*akp1-one )
336  a( k-1, k-1 ) = akp1 / d
337  a( k, k ) = ak / d
338  a( k, k-1 ) = -akkp1 / d
339 *
340 * Compute columns K-1 and K of the inverse.
341 *
342  IF( k.LT.n ) THEN
343  CALL zcopy( n-k, a( k+1, k ), 1, work, 1 )
344  CALL zhemv( uplo, n-k, -cone, a( k+1, k+1 ), lda, work,
345  $ 1, zero, a( k+1, k ), 1 )
346  a( k, k ) = a( k, k ) - dble( zdotc( n-k, work, 1,
347  $ a( k+1, k ), 1 ) )
348  a( k, k-1 ) = a( k, k-1 ) -
349  $ zdotc( n-k, a( k+1, k ), 1, a( k+1, k-1 ),
350  $ 1 )
351  CALL zcopy( n-k, a( k+1, k-1 ), 1, work, 1 )
352  CALL zhemv( uplo, n-k, -cone, a( k+1, k+1 ), lda, work,
353  $ 1, zero, a( k+1, k-1 ), 1 )
354  a( k-1, k-1 ) = a( k-1, k-1 ) -
355  $ dble( zdotc( n-k, work, 1, a( k+1, k-1 ),
356  $ 1 ) )
357  END IF
358  kstep = 2
359  END IF
360 *
361  kp = abs( ipiv( k ) )
362  IF( kp.NE.k ) THEN
363 *
364 * Interchange rows and columns K and KP in the trailing
365 * submatrix A(k-1:n,k-1:n)
366 *
367  IF( kp.LT.n )
368  $ CALL zswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
369  DO 70 j = k + 1, kp - 1
370  temp = dconjg( a( j, k ) )
371  a( j, k ) = dconjg( a( kp, j ) )
372  a( kp, j ) = temp
373  70 CONTINUE
374  a( kp, k ) = dconjg( a( kp, k ) )
375  temp = a( k, k )
376  a( k, k ) = a( kp, kp )
377  a( kp, kp ) = temp
378  IF( kstep.EQ.2 ) THEN
379  temp = a( k, k-1 )
380  a( k, k-1 ) = a( kp, k-1 )
381  a( kp, k-1 ) = temp
382  END IF
383  END IF
384 *
385  k = k - kstep
386  GO TO 60
387  80 CONTINUE
388  END IF
389 *
390  RETURN
391 *
392 * End of ZHETRI
393 *
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: