LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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
November 2011

Definition at line 111 of file chptri.f.

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 
chpmv
subroutine chpmv(UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY)
CHPMV
Definition: chpmv.f:151
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
cdotc
complex function cdotc(N, CX, INCX, CY, INCY)
CDOTC
Definition: cdotc.f:54
ccopy
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:52
cswap
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:52
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function: