LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ chptri()

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

CHPTRI

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

Purpose:
 CHPTRI 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 CHPTRF.
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 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 CHPTRF,
          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 CHPTRF.
[out]WORK
          WORK is COMPLEX 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
December 2016

Definition at line 111 of file chptri.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 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 ),
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 cswap( 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 = conjg( ap( kc+j-1 ) )
281  ap( kc+j-1 ) = conjg( ap( kx ) )
282  ap( kx ) = temp
283  40 CONTINUE
284  ap( kc+kp-1 ) = conjg( 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 / real( ap( kc ) )
325 *
326 * Compute column K of the inverse.
327 *
328  IF( k.LT.n ) THEN
329  CALL ccopy( n-k, ap( kc+1 ), 1, work, 1 )
330  CALL chpmv( uplo, n-k, -cone, ap( kc+n-k+1 ), work, 1,
331  $ zero, ap( kc+1 ), 1 )
332  ap( kc ) = ap( kc ) - real( cdotc( 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 = real( ap( kcnext ) ) / t
344  akp1 = real( 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 ccopy( n-k, ap( kc+1 ), 1, work, 1 )
355  CALL chpmv( uplo, n-k, -cone, ap( kc+( n-k+1 ) ), work,
356  $ 1, zero, ap( kc+1 ), 1 )
357  ap( kc ) = ap( kc ) - real( cdotc( n-k, work, 1,
358  $ ap( kc+1 ), 1 ) )
359  ap( kcnext+1 ) = ap( kcnext+1 ) -
360  $ cdotc( n-k, ap( kc+1 ), 1,
361  $ ap( kcnext+2 ), 1 )
362  CALL ccopy( n-k, ap( kcnext+2 ), 1, work, 1 )
363  CALL chpmv( uplo, n-k, -cone, ap( kc+( n-k+1 ) ), work,
364  $ 1, zero, ap( kcnext+2 ), 1 )
365  ap( kcnext ) = ap( kcnext ) -
366  $ real( cdotc( 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 cswap( 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 = conjg( ap( kc+j-k ) )
386  ap( kc+j-k ) = conjg( ap( kx ) )
387  ap( kx ) = temp
388  70 CONTINUE
389  ap( kc+kp-k ) = conjg( 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 CHPTRI
409 *
Here is the call graph for this function:
Here is the caller graph for this function:
chpmv
subroutine chpmv(UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY)
CHPMV
Definition: chpmv.f:151
cdotc
complex function cdotc(N, CX, INCX, CY, INCY)
CDOTC
Definition: cdotc.f:85
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
cswap
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:83
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
ccopy
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:83