LAPACK  3.8.0 LAPACK: Linear Algebra PACKage

## ◆ zhptri()

 subroutine zhptri ( character UPLO, integer N, complex*16, dimension( * ) AP, integer, dimension( * ) IPIV, complex*16, dimension( * ) WORK, integer INFO )

ZHPTRI

Purpose:
``` ZHPTRI computes the inverse of a complex Hermitian indefinite matrix
A in packed storage using the factorization A = U*D*U**H or
A = L*D*L**H computed by ZHPTRF.```
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] AP ``` AP is COMPLEX*16 array, dimension (N*(N+1)/2) On entry, the block diagonal matrix D and the multipliers used to obtain the factor U or L as computed by ZHPTRF, stored as a packed triangular matrix. On exit, if INFO = 0, the (Hermitian) inverse of the original matrix, stored as a packed triangular matrix. The j-th column of inv(A) is stored in the array AP as follows: if UPLO = 'U', AP(i + (j-1)*j/2) = inv(A)(i,j) for 1<=i<=j; if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = inv(A)(i,j) for j<=i<=n.``` [in] IPIV ``` IPIV is INTEGER array, dimension (N) Details of the interchanges and the block structure of D as determined by ZHPTRF.``` [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.```
Date
December 2016

Definition at line 111 of file zhptri.f.

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