LAPACK  3.9.1
LAPACK: Linear Algebra PACKage
zhptri.f
Go to the documentation of this file.
1 *> \brief \b ZHPTRI
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZHPTRI + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhptri.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhptri.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhptri.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZHPTRI( UPLO, N, AP, IPIV, WORK, INFO )
22 *
23 * .. Scalar Arguments ..
24 * CHARACTER UPLO
25 * INTEGER INFO, N
26 * ..
27 * .. Array Arguments ..
28 * INTEGER IPIV( * )
29 * COMPLEX*16 AP( * ), WORK( * )
30 * ..
31 *
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *>
38 *> ZHPTRI computes the inverse of a complex Hermitian indefinite matrix
39 *> A in packed storage using the factorization A = U*D*U**H or
40 *> A = L*D*L**H computed by ZHPTRF.
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] AP
62 *> \verbatim
63 *> AP is COMPLEX*16 array, dimension (N*(N+1)/2)
64 *> On entry, the block diagonal matrix D and the multipliers
65 *> used to obtain the factor U or L as computed by ZHPTRF,
66 *> stored as a packed triangular matrix.
67 *>
68 *> On exit, if INFO = 0, the (Hermitian) inverse of the original
69 *> matrix, stored as a packed triangular matrix. The j-th column
70 *> of inv(A) is stored in the array AP as follows:
71 *> if UPLO = 'U', AP(i + (j-1)*j/2) = inv(A)(i,j) for 1<=i<=j;
72 *> if UPLO = 'L',
73 *> AP(i + (j-1)*(2n-j)/2) = inv(A)(i,j) for j<=i<=n.
74 *> \endverbatim
75 *>
76 *> \param[in] IPIV
77 *> \verbatim
78 *> IPIV is INTEGER array, dimension (N)
79 *> Details of the interchanges and the block structure of D
80 *> as determined by ZHPTRF.
81 *> \endverbatim
82 *>
83 *> \param[out] WORK
84 *> \verbatim
85 *> WORK is COMPLEX*16 array, dimension (N)
86 *> \endverbatim
87 *>
88 *> \param[out] INFO
89 *> \verbatim
90 *> INFO is INTEGER
91 *> = 0: successful exit
92 *> < 0: if INFO = -i, the i-th argument had an illegal value
93 *> > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
94 *> inverse could not be computed.
95 *> \endverbatim
96 *
97 * Authors:
98 * ========
99 *
100 *> \author Univ. of Tennessee
101 *> \author Univ. of California Berkeley
102 *> \author Univ. of Colorado Denver
103 *> \author NAG Ltd.
104 *
105 *> \ingroup complex16OTHERcomputational
106 *
107 * =====================================================================
108  SUBROUTINE zhptri( UPLO, N, AP, IPIV, WORK, INFO )
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 *
407  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 zhpmv(UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY)
ZHPMV
Definition: zhpmv.f:149
subroutine zhptri(UPLO, N, AP, IPIV, WORK, INFO)
ZHPTRI
Definition: zhptri.f:109