LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine zhptri ( character  UPLO,
integer  N,
complex*16, dimension( * )  AP,
integer, dimension( * )  IPIV,
complex*16, dimension( * )  WORK,
integer  INFO 
)

ZHPTRI

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

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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 111 of file zhptri.f.

111 *
112 * -- LAPACK computational routine (version 3.4.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 * November 2011
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 *
zcopy
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:52
zswap
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:52
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
zdotc
complex *16 function zdotc(N, ZX, INCX, ZY, INCY)
ZDOTC
Definition: zdotc.f:54
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
zhpmv
subroutine zhpmv(UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY)
ZHPMV
Definition: zhpmv.f:151

Here is the call graph for this function:

Here is the caller graph for this function: