LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
zhetri.f
Go to the documentation of this file.
1 *> \brief \b ZHETRI
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZHETRI + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetri.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetri.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetri.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZHETRI( UPLO, N, A, LDA, IPIV, WORK, INFO )
22 *
23 * .. Scalar Arguments ..
24 * CHARACTER UPLO
25 * INTEGER INFO, LDA, N
26 * ..
27 * .. Array Arguments ..
28 * INTEGER IPIV( * )
29 * COMPLEX*16 A( LDA, * ), WORK( * )
30 * ..
31 *
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *>
38 *> ZHETRI computes the inverse of a complex Hermitian indefinite matrix
39 *> A using the factorization A = U*D*U**H or A = L*D*L**H computed by
40 *> ZHETRF.
41 *> \endverbatim
42 *
43 * Arguments:
44 * ==========
45 *
46 *> \param[in] UPLO
47 *> \verbatim
48 *> UPLO is CHARACTER*1
49 *> Specifies whether the details of the factorization are stored
50 *> as an upper or lower triangular matrix.
51 *> = 'U': Upper triangular, form is A = U*D*U**H;
52 *> = 'L': Lower triangular, form is A = L*D*L**H.
53 *> \endverbatim
54 *>
55 *> \param[in] N
56 *> \verbatim
57 *> N is INTEGER
58 *> The order of the matrix A. N >= 0.
59 *> \endverbatim
60 *>
61 *> \param[in,out] A
62 *> \verbatim
63 *> A is COMPLEX*16 array, dimension (LDA,N)
64 *> On entry, the block diagonal matrix D and the multipliers
65 *> used to obtain the factor U or L as computed by ZHETRF.
66 *>
67 *> On exit, if INFO = 0, the (Hermitian) inverse of the original
68 *> matrix. If UPLO = 'U', the upper triangular part of the
69 *> inverse is formed and the part of A below the diagonal is not
70 *> referenced; if UPLO = 'L' the lower triangular part of the
71 *> inverse is formed and the part of A above the diagonal is
72 *> not referenced.
73 *> \endverbatim
74 *>
75 *> \param[in] LDA
76 *> \verbatim
77 *> LDA is INTEGER
78 *> The leading dimension of the array A. LDA >= max(1,N).
79 *> \endverbatim
80 *>
81 *> \param[in] IPIV
82 *> \verbatim
83 *> IPIV is INTEGER array, dimension (N)
84 *> Details of the interchanges and the block structure of D
85 *> as determined by ZHETRF.
86 *> \endverbatim
87 *>
88 *> \param[out] WORK
89 *> \verbatim
90 *> WORK is COMPLEX*16 array, dimension (N)
91 *> \endverbatim
92 *>
93 *> \param[out] INFO
94 *> \verbatim
95 *> INFO is INTEGER
96 *> = 0: successful exit
97 *> < 0: if INFO = -i, the i-th argument had an illegal value
98 *> > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
99 *> inverse could not be computed.
100 *> \endverbatim
101 *
102 * Authors:
103 * ========
104 *
105 *> \author Univ. of Tennessee
106 *> \author Univ. of California Berkeley
107 *> \author Univ. of Colorado Denver
108 *> \author NAG Ltd.
109 *
110 *> \ingroup complex16HEcomputational
111 *
112 * =====================================================================
113  SUBROUTINE zhetri( UPLO, N, A, LDA, IPIV, WORK, INFO )
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 *
394  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:81
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
subroutine zhetri(UPLO, N, A, LDA, IPIV, WORK, INFO)
ZHETRI
Definition: zhetri.f:114