LAPACK  3.10.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.

Definition at line 108 of file chptri.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 AP( * ), WORK( * )
121 * ..
122 *
123 * =====================================================================
124 *
125 * .. Parameters ..
126  REAL ONE
127  COMPLEX CONE, ZERO
128  parameter( one = 1.0e+0, cone = ( 1.0e+0, 0.0e+0 ),
129  $ zero = ( 0.0e+0, 0.0e+0 ) )
130 * ..
131 * .. Local Scalars ..
132  LOGICAL UPPER
133  INTEGER J, K, KC, KCNEXT, KP, KPC, KSTEP, KX, NPP
134  REAL AK, AKP1, D, T
135  COMPLEX AKKP1, TEMP
136 * ..
137 * .. External Functions ..
138  LOGICAL LSAME
139  COMPLEX CDOTC
140  EXTERNAL lsame, cdotc
141 * ..
142 * .. External Subroutines ..
143  EXTERNAL ccopy, chpmv, cswap, xerbla
144 * ..
145 * .. Intrinsic Functions ..
146  INTRINSIC abs, conjg, real
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( 'CHPTRI', -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 / real( ap( kc+k-1 ) )
218 *
219 * Compute column K of the inverse.
220 *
221  IF( k.GT.1 ) THEN
222  CALL ccopy( k-1, ap( kc ), 1, work, 1 )
223  CALL chpmv( uplo, k-1, -cone, ap, work, 1, zero,
224  $ ap( kc ), 1 )
225  ap( kc+k-1 ) = ap( kc+k-1 ) -
226  $ real( cdotc( 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 = real( ap( kc+k-1 ) ) / t
237  akp1 = real( 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 ccopy( k-1, ap( kc ), 1, work, 1 )
248  CALL chpmv( uplo, k-1, -cone, ap, work, 1, zero,
249  $ ap( kc ), 1 )
250  ap( kc+k-1 ) = ap( kc+k-1 ) -
251  $ real( cdotc( k-1, work, 1, ap( kc ), 1 ) )
252  ap( kcnext+k-1 ) = ap( kcnext+k-1 ) -
253  $ cdotc( k-1, ap( kc ), 1, ap( kcnext ),
254  $ 1 )
255  CALL ccopy( k-1, ap( kcnext ), 1, work, 1 )
256  CALL chpmv( uplo, k-1, -cone, ap, work, 1, zero,
257  $ ap( kcnext ), 1 )
258  ap( kcnext+k ) = ap( kcnext+k ) -
259  $ real( cdotc( 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 cswap( 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 = conjg( ap( kc+j-1 ) )
278  ap( kc+j-1 ) = conjg( ap( kx ) )
279  ap( kx ) = temp
280  40 CONTINUE
281  ap( kc+kp-1 ) = conjg( 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 / real( ap( kc ) )
322 *
323 * Compute column K of the inverse.
324 *
325  IF( k.LT.n ) THEN
326  CALL ccopy( n-k, ap( kc+1 ), 1, work, 1 )
327  CALL chpmv( uplo, n-k, -cone, ap( kc+n-k+1 ), work, 1,
328  $ zero, ap( kc+1 ), 1 )
329  ap( kc ) = ap( kc ) - real( cdotc( 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 = real( ap( kcnext ) ) / t
341  akp1 = real( 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 ccopy( n-k, ap( kc+1 ), 1, work, 1 )
352  CALL chpmv( uplo, n-k, -cone, ap( kc+( n-k+1 ) ), work,
353  $ 1, zero, ap( kc+1 ), 1 )
354  ap( kc ) = ap( kc ) - real( cdotc( n-k, work, 1,
355  $ 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 ),
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 cswap( 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 = conjg( ap( kc+j-k ) )
383  ap( kc+j-k ) = conjg( ap( kx ) )
384  ap( kx ) = temp
385  70 CONTINUE
386  ap( kc+kp-k ) = conjg( 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 CHPTRI
406 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
complex function cdotc(N, CX, INCX, CY, INCY)
CDOTC
Definition: cdotc.f:83
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:81
subroutine chpmv(UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY)
CHPMV
Definition: chpmv.f:149
Here is the call graph for this function:
Here is the caller graph for this function: