LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
chptri.f
Go to the documentation of this file.
1 *> \brief \b CHPTRI
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CHPTRI + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chptri.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chptri.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chptri.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CHPTRI( 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 AP( * ), WORK( * )
30 * ..
31 *
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *>
38 *> CHPTRI 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 CHPTRF.
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 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 CHPTRF,
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 CHPTRF.
81 *> \endverbatim
82 *>
83 *> \param[out] WORK
84 *> \verbatim
85 *> WORK is COMPLEX 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 *> \date November 2011
106 *
107 *> \ingroup complexOTHERcomputational
108 *
109 * =====================================================================
110  SUBROUTINE chptri( UPLO, N, AP, IPIV, WORK, INFO )
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 ap( * ), work( * )
124 * ..
125 *
126 * =====================================================================
127 *
128 * .. Parameters ..
129  REAL one
130  COMPLEX cone, zero
131  parameter( one = 1.0e+0, cone = ( 1.0e+0, 0.0e+0 ),
132  $ zero = ( 0.0e+0, 0.0e+0 ) )
133 * ..
134 * .. Local Scalars ..
135  LOGICAL upper
136  INTEGER j, k, kc, kcnext, kp, kpc, kstep, kx, npp
137  REAL ak, akp1, d, t
138  COMPLEX akkp1, temp
139 * ..
140 * .. External Functions ..
141  LOGICAL lsame
142  COMPLEX cdotc
143  EXTERNAL lsame, cdotc
144 * ..
145 * .. External Subroutines ..
146  EXTERNAL ccopy, chpmv, cswap, xerbla
147 * ..
148 * .. Intrinsic Functions ..
149  INTRINSIC abs, conjg, real
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( 'CHPTRI', -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 / REAL( AP( KC+K-1 ) )
221 *
222 * Compute column K of the inverse.
223 *
224  IF( k.GT.1 ) THEN
225  CALL ccopy( k-1, ap( kc ), 1, work, 1 )
226  CALL chpmv( uplo, k-1, -cone, ap, work, 1, zero,
227  $ ap( kc ), 1 )
228  ap( kc+k-1 ) = ap( kc+k-1 ) -
229  $ REAL( CDOTC( 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 = REAL( AP( KC+K-1 ) ) / t
240  akp1 = REAL( 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 ccopy( k-1, ap( kc ), 1, work, 1 )
251  CALL chpmv( uplo, k-1, -cone, ap, work, 1, zero,
252  $ ap( kc ), 1 )
253  ap( kc+k-1 ) = ap( kc+k-1 ) -
254  $ REAL( CDOTC( K-1, WORK, 1, AP( KC ), 1 ) )
255  ap( kcnext+k-1 ) = ap( kcnext+k-1 ) -
256  $ cdotc( k-1, ap( kc ), 1, ap( kcnext ),
257  $ 1 )
258  CALL ccopy( k-1, ap( kcnext ), 1, work, 1 )
259  CALL chpmv( uplo, k-1, -cone, ap, work, 1, zero,
260  $ ap( kcnext ), 1 )
261  ap( kcnext+k ) = ap( kcnext+k ) -
262  $ REAL( CDOTC( K-1, WORK, 1, AP( KCNEXT ), $ 1 ) )
263  END IF
264  kstep = 2
265  kcnext = kcnext + k + 1
266  END IF
267 *
268  kp = abs( ipiv( k ) )
269  IF( kp.NE.k ) THEN
270 *
271 * Interchange rows and columns K and KP in the leading
272 * submatrix A(1:k+1,1:k+1)
273 *
274  kpc = ( kp-1 )*kp / 2 + 1
275  CALL cswap( kp-1, ap( kc ), 1, ap( kpc ), 1 )
276  kx = kpc + kp - 1
277  DO 40 j = kp + 1, k - 1
278  kx = kx + j - 1
279  temp = conjg( ap( kc+j-1 ) )
280  ap( kc+j-1 ) = conjg( ap( kx ) )
281  ap( kx ) = temp
282  40 continue
283  ap( kc+kp-1 ) = conjg( ap( kc+kp-1 ) )
284  temp = ap( kc+k-1 )
285  ap( kc+k-1 ) = ap( kpc+kp-1 )
286  ap( kpc+kp-1 ) = temp
287  IF( kstep.EQ.2 ) THEN
288  temp = ap( kc+k+k-1 )
289  ap( kc+k+k-1 ) = ap( kc+k+kp-1 )
290  ap( kc+k+kp-1 ) = temp
291  END IF
292  END IF
293 *
294  k = k + kstep
295  kc = kcnext
296  go to 30
297  50 continue
298 *
299  ELSE
300 *
301 * Compute inv(A) from the factorization A = L*D*L**H.
302 *
303 * K is the main loop index, increasing from 1 to N in steps of
304 * 1 or 2, depending on the size of the diagonal blocks.
305 *
306  npp = n*( n+1 ) / 2
307  k = n
308  kc = npp
309  60 continue
310 *
311 * If K < 1, exit from loop.
312 *
313  IF( k.LT.1 )
314  $ go to 80
315 *
316  kcnext = kc - ( n-k+2 )
317  IF( ipiv( k ).GT.0 ) THEN
318 *
319 * 1 x 1 diagonal block
320 *
321 * Invert the diagonal block.
322 *
323  ap( kc ) = one / REAL( AP( KC ) )
324 *
325 * Compute column K of the inverse.
326 *
327  IF( k.LT.n ) THEN
328  CALL ccopy( n-k, ap( kc+1 ), 1, work, 1 )
329  CALL chpmv( uplo, n-k, -cone, ap( kc+n-k+1 ), work, 1,
330  $ zero, ap( kc+1 ), 1 )
331  ap( kc ) = ap( kc ) - REAL( CDOTC( N-K, WORK, 1, $ AP( KC+1 ), 1 ) )
332  END IF
333  kstep = 1
334  ELSE
335 *
336 * 2 x 2 diagonal block
337 *
338 * Invert the diagonal block.
339 *
340  t = abs( ap( kcnext+1 ) )
341  ak = REAL( AP( KCNEXT ) ) / t
342  akp1 = REAL( AP( KC ) ) / t
343  akkp1 = ap( kcnext+1 ) / t
344  d = t*( ak*akp1-one )
345  ap( kcnext ) = akp1 / d
346  ap( kc ) = ak / d
347  ap( kcnext+1 ) = -akkp1 / d
348 *
349 * Compute columns K-1 and K of the inverse.
350 *
351  IF( k.LT.n ) THEN
352  CALL ccopy( n-k, ap( kc+1 ), 1, work, 1 )
353  CALL chpmv( uplo, n-k, -cone, ap( kc+( n-k+1 ) ), work,
354  $ 1, zero, ap( kc+1 ), 1 )
355  ap( kc ) = ap( kc ) - REAL( CDOTC( N-K, WORK, 1, $ AP( KC+1 ), 1 ) )
356  ap( kcnext+1 ) = ap( kcnext+1 ) -
357  $ cdotc( n-k, ap( kc+1 ), 1,
358  $ ap( kcnext+2 ), 1 )
359  CALL ccopy( n-k, ap( kcnext+2 ), 1, work, 1 )
360  CALL chpmv( uplo, n-k, -cone, ap( kc+( n-k+1 ) ), work,
361  $ 1, zero, ap( kcnext+2 ), 1 )
362  ap( kcnext ) = ap( kcnext ) -
363  $ REAL( CDOTC( N-K, WORK, 1, AP( KCNEXT+2 ), $ 1 ) )
364  END IF
365  kstep = 2
366  kcnext = kcnext - ( n-k+3 )
367  END IF
368 *
369  kp = abs( ipiv( k ) )
370  IF( kp.NE.k ) THEN
371 *
372 * Interchange rows and columns K and KP in the trailing
373 * submatrix A(k-1:n,k-1:n)
374 *
375  kpc = npp - ( n-kp+1 )*( n-kp+2 ) / 2 + 1
376  IF( kp.LT.n )
377  $ CALL cswap( n-kp, ap( kc+kp-k+1 ), 1, ap( kpc+1 ), 1 )
378  kx = kc + kp - k
379  DO 70 j = k + 1, kp - 1
380  kx = kx + n - j + 1
381  temp = conjg( ap( kc+j-k ) )
382  ap( kc+j-k ) = conjg( ap( kx ) )
383  ap( kx ) = temp
384  70 continue
385  ap( kc+kp-k ) = conjg( ap( kc+kp-k ) )
386  temp = ap( kc )
387  ap( kc ) = ap( kpc )
388  ap( kpc ) = temp
389  IF( kstep.EQ.2 ) THEN
390  temp = ap( kc-n+k-1 )
391  ap( kc-n+k-1 ) = ap( kc-n+kp-1 )
392  ap( kc-n+kp-1 ) = temp
393  END IF
394  END IF
395 *
396  k = k - kstep
397  kc = kcnext
398  go to 60
399  80 continue
400  END IF
401 *
402  return
403 *
404 * End of CHPTRI
405 *
406  END
407