1       SUBROUTINE ZHETRI( UPLO, N, A, LDA, IPIV, WORK, INFO )
  2 *
  3 *  -- LAPACK routine (version 3.3.1) --
  4 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  5 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  6 *  -- April 2011                                                      --
  7 *
  8 *     .. Scalar Arguments ..
  9       CHARACTER          UPLO
 10       INTEGER            INFO, LDA, N
 11 *     ..
 12 *     .. Array Arguments ..
 13       INTEGER            IPIV( * )
 14       COMPLEX*16         A( LDA, * ), WORK( * )
 15 *     ..
 16 *
 17 *  Purpose
 18 *  =======
 19 *
 20 *  ZHETRI computes the inverse of a complex Hermitian indefinite matrix
 21 *  A using the factorization A = U*D*U**H or A = L*D*L**H computed by
 22 *  ZHETRF.
 23 *
 24 *  Arguments
 25 *  =========
 26 *
 27 *  UPLO    (input) CHARACTER*1
 28 *          Specifies whether the details of the factorization are stored
 29 *          as an upper or lower triangular matrix.
 30 *          = 'U':  Upper triangular, form is A = U*D*U**H;
 31 *          = 'L':  Lower triangular, form is A = L*D*L**H.
 32 *
 33 *  N       (input) INTEGER
 34 *          The order of the matrix A.  N >= 0.
 35 *
 36 *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
 37 *          On entry, the block diagonal matrix D and the multipliers
 38 *          used to obtain the factor U or L as computed by ZHETRF.
 39 *
 40 *          On exit, if INFO = 0, the (Hermitian) inverse of the original
 41 *          matrix.  If UPLO = 'U', the upper triangular part of the
 42 *          inverse is formed and the part of A below the diagonal is not
 43 *          referenced; if UPLO = 'L' the lower triangular part of the
 44 *          inverse is formed and the part of A above the diagonal is
 45 *          not referenced.
 46 *
 47 *  LDA     (input) INTEGER
 48 *          The leading dimension of the array A.  LDA >= max(1,N).
 49 *
 50 *  IPIV    (input) INTEGER array, dimension (N)
 51 *          Details of the interchanges and the block structure of D
 52 *          as determined by ZHETRF.
 53 *
 54 *  WORK    (workspace) COMPLEX*16 array, dimension (N)
 55 *
 56 *  INFO    (output) INTEGER
 57 *          = 0: successful exit
 58 *          < 0: if INFO = -i, the i-th argument had an illegal value
 59 *          > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
 60 *               inverse could not be computed.
 61 *
 62 *  =====================================================================
 63 *
 64 *     .. Parameters ..
 65       DOUBLE PRECISION   ONE
 66       COMPLEX*16         CONE, ZERO
 67       PARAMETER          ( ONE = 1.0D+0, CONE = ( 1.0D+00.0D+0 ),
 68      $                   ZERO = ( 0.0D+00.0D+0 ) )
 69 *     ..
 70 *     .. Local Scalars ..
 71       LOGICAL            UPPER
 72       INTEGER            J, K, KP, KSTEP
 73       DOUBLE PRECISION   AK, AKP1, D, T
 74       COMPLEX*16         AKKP1, TEMP
 75 *     ..
 76 *     .. External Functions ..
 77       LOGICAL            LSAME
 78       COMPLEX*16         ZDOTC
 79       EXTERNAL           LSAME, ZDOTC
 80 *     ..
 81 *     .. External Subroutines ..
 82       EXTERNAL           XERBLA, ZCOPY, ZHEMV, ZSWAP
 83 *     ..
 84 *     .. Intrinsic Functions ..
 85       INTRINSIC          ABSDBLEDCONJGMAX
 86 *     ..
 87 *     .. Executable Statements ..
 88 *
 89 *     Test the input parameters.
 90 *
 91       INFO = 0
 92       UPPER = LSAME( UPLO, 'U' )
 93       IF.NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 94          INFO = -1
 95       ELSE IF( N.LT.0 ) THEN
 96          INFO = -2
 97       ELSE IF( LDA.LT.MAX1, N ) ) THEN
 98          INFO = -4
 99       END IF
100       IF( INFO.NE.0 ) THEN
101          CALL XERBLA( 'ZHETRI'-INFO )
102          RETURN
103       END IF
104 *
105 *     Quick return if possible
106 *
107       IF( N.EQ.0 )
108      $   RETURN
109 *
110 *     Check that the diagonal matrix D is nonsingular.
111 *
112       IF( UPPER ) THEN
113 *
114 *        Upper triangular storage: examine D from bottom to top
115 *
116          DO 10 INFO = N, 1-1
117             IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
118      $         RETURN
119    10    CONTINUE
120       ELSE
121 *
122 *        Lower triangular storage: examine D from top to bottom.
123 *
124          DO 20 INFO = 1, N
125             IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
126      $         RETURN
127    20    CONTINUE
128       END IF
129       INFO = 0
130 *
131       IF( UPPER ) THEN
132 *
133 *        Compute inv(A) from the factorization A = U*D*U**H.
134 *
135 *        K is the main loop index, increasing from 1 to N in steps of
136 *        1 or 2, depending on the size of the diagonal blocks.
137 *
138          K = 1
139    30    CONTINUE
140 *
141 *        If K > N, exit from loop.
142 *
143          IF( K.GT.N )
144      $      GO TO 50
145 *
146          IF( IPIV( K ).GT.0 ) THEN
147 *
148 *           1 x 1 diagonal block
149 *
150 *           Invert the diagonal block.
151 *
152             A( K, K ) = ONE / DBLE( A( K, K ) )
153 *
154 *           Compute column K of the inverse.
155 *
156             IF( K.GT.1 ) THEN
157                CALL ZCOPY( K-1, A( 1, K ), 1, WORK, 1 )
158                CALL ZHEMV( UPLO, K-1-CONE, A, LDA, WORK, 1, ZERO,
159      $                     A( 1, K ), 1 )
160                A( K, K ) = A( K, K ) - DBLE( ZDOTC( K-1, WORK, 1, A( 1,
161      $                     K ), 1 ) )
162             END IF
163             KSTEP = 1
164          ELSE
165 *
166 *           2 x 2 diagonal block
167 *
168 *           Invert the diagonal block.
169 *
170             T = ABS( A( K, K+1 ) )
171             AK = DBLE( A( K, K ) ) / T
172             AKP1 = DBLE( A( K+1, K+1 ) ) / T
173             AKKP1 = A( K, K+1 ) / T
174             D = T*( AK*AKP1-ONE )
175             A( K, K ) = AKP1 / D
176             A( K+1, K+1 ) = AK / D
177             A( K, K+1 ) = -AKKP1 / D
178 *
179 *           Compute columns K and K+1 of the inverse.
180 *
181             IF( K.GT.1 ) THEN
182                CALL ZCOPY( K-1, A( 1, K ), 1, WORK, 1 )
183                CALL ZHEMV( UPLO, K-1-CONE, A, LDA, WORK, 1, ZERO,
184      $                     A( 1, K ), 1 )
185                A( K, K ) = A( K, K ) - DBLE( ZDOTC( K-1, WORK, 1, A( 1,
186      $                     K ), 1 ) )
187                A( K, K+1 ) = A( K, K+1 ) -
188      $                       ZDOTC( K-1, A( 1, K ), 1, A( 1, K+1 ), 1 )
189                CALL ZCOPY( K-1, A( 1, K+1 ), 1, WORK, 1 )
190                CALL ZHEMV( UPLO, K-1-CONE, A, LDA, WORK, 1, ZERO,
191      $                     A( 1, K+1 ), 1 )
192                A( K+1, K+1 ) = A( K+1, K+1 ) -
193      $                         DBLE( ZDOTC( K-1, WORK, 1, A( 1, K+1 ),
194      $                         1 ) )
195             END IF
196             KSTEP = 2
197          END IF
198 *
199          KP = ABS( IPIV( K ) )
200          IF( KP.NE.K ) THEN
201 *
202 *           Interchange rows and columns K and KP in the leading
203 *           submatrix A(1:k+1,1:k+1)
204 *
205             CALL ZSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
206             DO 40 J = KP + 1, K - 1
207                TEMP = DCONJG( A( J, K ) )
208                A( J, K ) = DCONJG( A( KP, J ) )
209                A( KP, J ) = TEMP
210    40       CONTINUE
211             A( KP, K ) = DCONJG( A( KP, K ) )
212             TEMP = A( K, K )
213             A( K, K ) = A( KP, KP )
214             A( KP, KP ) = TEMP
215             IF( KSTEP.EQ.2 ) THEN
216                TEMP = A( K, K+1 )
217                A( K, K+1 ) = A( KP, K+1 )
218                A( KP, K+1 ) = TEMP
219             END IF
220          END IF
221 *
222          K = K + KSTEP
223          GO TO 30
224    50    CONTINUE
225 *
226       ELSE
227 *
228 *        Compute inv(A) from the factorization A = L*D*L**H.
229 *
230 *        K is the main loop index, increasing from 1 to N in steps of
231 *        1 or 2, depending on the size of the diagonal blocks.
232 *
233          K = N
234    60    CONTINUE
235 *
236 *        If K < 1, exit from loop.
237 *
238          IF( K.LT.1 )
239      $      GO TO 80
240 *
241          IF( IPIV( K ).GT.0 ) THEN
242 *
243 *           1 x 1 diagonal block
244 *
245 *           Invert the diagonal block.
246 *
247             A( K, K ) = ONE / DBLE( A( K, K ) )
248 *
249 *           Compute column K of the inverse.
250 *
251             IF( K.LT.N ) THEN
252                CALL ZCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
253                CALL ZHEMV( UPLO, N-K, -CONE, A( K+1, K+1 ), LDA, WORK,
254      $                     1, ZERO, A( K+1, K ), 1 )
255                A( K, K ) = A( K, K ) - DBLE( ZDOTC( N-K, WORK, 1,
256      $                     A( K+1, K ), 1 ) )
257             END IF
258             KSTEP = 1
259          ELSE
260 *
261 *           2 x 2 diagonal block
262 *
263 *           Invert the diagonal block.
264 *
265             T = ABS( A( K, K-1 ) )
266             AK = DBLE( A( K-1, K-1 ) ) / T
267             AKP1 = DBLE( A( K, K ) ) / T
268             AKKP1 = A( K, K-1 ) / T
269             D = T*( AK*AKP1-ONE )
270             A( K-1, K-1 ) = AKP1 / D
271             A( K, K ) = AK / D
272             A( K, K-1 ) = -AKKP1 / D
273 *
274 *           Compute columns K-1 and K of the inverse.
275 *
276             IF( K.LT.N ) THEN
277                CALL ZCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
278                CALL ZHEMV( UPLO, N-K, -CONE, A( K+1, K+1 ), LDA, WORK,
279      $                     1, ZERO, A( K+1, K ), 1 )
280                A( K, K ) = A( K, K ) - DBLE( ZDOTC( N-K, WORK, 1,
281      $                     A( K+1, K ), 1 ) )
282                A( K, K-1 ) = A( K, K-1 ) -
283      $                       ZDOTC( N-K, A( K+1, K ), 1, A( K+1, K-1 ),
284      $                       1 )
285                CALL ZCOPY( N-K, A( K+1, K-1 ), 1, WORK, 1 )
286                CALL ZHEMV( UPLO, N-K, -CONE, A( K+1, K+1 ), LDA, WORK,
287      $                     1, ZERO, A( K+1, K-1 ), 1 )
288                A( K-1, K-1 ) = A( K-1, K-1 ) -
289      $                         DBLE( ZDOTC( N-K, WORK, 1, A( K+1, K-1 ),
290      $                         1 ) )
291             END IF
292             KSTEP = 2
293          END IF
294 *
295          KP = ABS( IPIV( K ) )
296          IF( KP.NE.K ) THEN
297 *
298 *           Interchange rows and columns K and KP in the trailing
299 *           submatrix A(k-1:n,k-1:n)
300 *
301             IF( KP.LT.N )
302      $         CALL ZSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
303             DO 70 J = K + 1, KP - 1
304                TEMP = DCONJG( A( J, K ) )
305                A( J, K ) = DCONJG( A( KP, J ) )
306                A( KP, J ) = TEMP
307    70       CONTINUE
308             A( KP, K ) = DCONJG( A( KP, K ) )
309             TEMP = A( K, K )
310             A( K, K ) = A( KP, KP )
311             A( KP, KP ) = TEMP
312             IF( KSTEP.EQ.2 ) THEN
313                TEMP = A( K, K-1 )
314                A( K, K-1 ) = A( KP, K-1 )
315                A( KP, K-1 ) = TEMP
316             END IF
317          END IF
318 *
319          K = K - KSTEP
320          GO TO 60
321    80    CONTINUE
322       END IF
323 *
324       RETURN
325 *
326 *     End of ZHETRI
327 *
328       END