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

◆ 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*
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: