LAPACK  3.10.1
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

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.

Definition at line 108 of file zhptri.f.

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