LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
ccopy
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
cdotc
complex function cdotc(n, cx, incx, cy, incy)
CDOTC
Definition cdotc.f:83
chpmv
subroutine chpmv(uplo, n, alpha, ap, x, incx, beta, y, incy)
CHPMV
Definition chpmv.f:149
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
cswap
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81
Here is the call graph for this function:
Here is the caller graph for this function: