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

◆ csptri()

subroutine csptri ( character  uplo,
integer  n,
complex, dimension( * )  ap,
integer, dimension( * )  ipiv,
complex, dimension( * )  work,
integer  info 
)

CSPTRI

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

Purpose:
 CSPTRI computes the inverse of a complex symmetric indefinite matrix
 A in packed storage using the factorization A = U*D*U**T or
 A = L*D*L**T computed by CSPTRF.
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**T;
          = 'L':  Lower triangular, form is A = L*D*L**T.
[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 CSPTRF,
          stored as a packed triangular matrix.

          On exit, if INFO = 0, the (symmetric) 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 CSPTRF.
[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 csptri.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 COMPLEX ONE, ZERO
127 parameter( one = ( 1.0e+0, 0.0e+0 ),
128 $ zero = ( 0.0e+0, 0.0e+0 ) )
129* ..
130* .. Local Scalars ..
131 LOGICAL UPPER
132 INTEGER J, K, KC, KCNEXT, KP, KPC, KSTEP, KX, NPP
133 COMPLEX AK, AKKP1, AKP1, D, T, TEMP
134* ..
135* .. External Functions ..
136 LOGICAL LSAME
137 COMPLEX CDOTU
138 EXTERNAL lsame, cdotu
139* ..
140* .. External Subroutines ..
141 EXTERNAL ccopy, cspmv, cswap, xerbla
142* ..
143* .. Intrinsic Functions ..
144 INTRINSIC abs
145* ..
146* .. Executable Statements ..
147*
148* Test the input parameters.
149*
150 info = 0
151 upper = lsame( uplo, 'U' )
152 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
153 info = -1
154 ELSE IF( n.LT.0 ) THEN
155 info = -2
156 END IF
157 IF( info.NE.0 ) THEN
158 CALL xerbla( 'CSPTRI', -info )
159 RETURN
160 END IF
161*
162* Quick return if possible
163*
164 IF( n.EQ.0 )
165 $ RETURN
166*
167* Check that the diagonal matrix D is nonsingular.
168*
169 IF( upper ) THEN
170*
171* Upper triangular storage: examine D from bottom to top
172*
173 kp = n*( n+1 ) / 2
174 DO 10 info = n, 1, -1
175 IF( ipiv( info ).GT.0 .AND. ap( kp ).EQ.zero )
176 $ RETURN
177 kp = kp - info
178 10 CONTINUE
179 ELSE
180*
181* Lower triangular storage: examine D from top to bottom.
182*
183 kp = 1
184 DO 20 info = 1, n
185 IF( ipiv( info ).GT.0 .AND. ap( kp ).EQ.zero )
186 $ RETURN
187 kp = kp + n - info + 1
188 20 CONTINUE
189 END IF
190 info = 0
191*
192 IF( upper ) THEN
193*
194* Compute inv(A) from the factorization A = U*D*U**T.
195*
196* K is the main loop index, increasing from 1 to N in steps of
197* 1 or 2, depending on the size of the diagonal blocks.
198*
199 k = 1
200 kc = 1
201 30 CONTINUE
202*
203* If K > N, exit from loop.
204*
205 IF( k.GT.n )
206 $ GO TO 50
207*
208 kcnext = kc + k
209 IF( ipiv( k ).GT.0 ) THEN
210*
211* 1 x 1 diagonal block
212*
213* Invert the diagonal block.
214*
215 ap( kc+k-1 ) = one / ap( kc+k-1 )
216*
217* Compute column K of the inverse.
218*
219 IF( k.GT.1 ) THEN
220 CALL ccopy( k-1, ap( kc ), 1, work, 1 )
221 CALL cspmv( uplo, k-1, -one, ap, work, 1, zero, ap( kc ),
222 $ 1 )
223 ap( kc+k-1 ) = ap( kc+k-1 ) -
224 $ cdotu( k-1, work, 1, ap( kc ), 1 )
225 END IF
226 kstep = 1
227 ELSE
228*
229* 2 x 2 diagonal block
230*
231* Invert the diagonal block.
232*
233 t = ap( kcnext+k-1 )
234 ak = ap( kc+k-1 ) / t
235 akp1 = ap( kcnext+k ) / t
236 akkp1 = ap( kcnext+k-1 ) / t
237 d = t*( ak*akp1-one )
238 ap( kc+k-1 ) = akp1 / d
239 ap( kcnext+k ) = ak / d
240 ap( kcnext+k-1 ) = -akkp1 / d
241*
242* Compute columns K and K+1 of the inverse.
243*
244 IF( k.GT.1 ) THEN
245 CALL ccopy( k-1, ap( kc ), 1, work, 1 )
246 CALL cspmv( uplo, k-1, -one, ap, work, 1, zero, ap( kc ),
247 $ 1 )
248 ap( kc+k-1 ) = ap( kc+k-1 ) -
249 $ cdotu( k-1, work, 1, ap( kc ), 1 )
250 ap( kcnext+k-1 ) = ap( kcnext+k-1 ) -
251 $ cdotu( k-1, ap( kc ), 1, ap( kcnext ),
252 $ 1 )
253 CALL ccopy( k-1, ap( kcnext ), 1, work, 1 )
254 CALL cspmv( uplo, k-1, -one, ap, work, 1, zero,
255 $ ap( kcnext ), 1 )
256 ap( kcnext+k ) = ap( kcnext+k ) -
257 $ cdotu( k-1, work, 1, ap( kcnext ), 1 )
258 END IF
259 kstep = 2
260 kcnext = kcnext + k + 1
261 END IF
262*
263 kp = abs( ipiv( k ) )
264 IF( kp.NE.k ) THEN
265*
266* Interchange rows and columns K and KP in the leading
267* submatrix A(1:k+1,1:k+1)
268*
269 kpc = ( kp-1 )*kp / 2 + 1
270 CALL cswap( kp-1, ap( kc ), 1, ap( kpc ), 1 )
271 kx = kpc + kp - 1
272 DO 40 j = kp + 1, k - 1
273 kx = kx + j - 1
274 temp = ap( kc+j-1 )
275 ap( kc+j-1 ) = ap( kx )
276 ap( kx ) = temp
277 40 CONTINUE
278 temp = ap( kc+k-1 )
279 ap( kc+k-1 ) = ap( kpc+kp-1 )
280 ap( kpc+kp-1 ) = temp
281 IF( kstep.EQ.2 ) THEN
282 temp = ap( kc+k+k-1 )
283 ap( kc+k+k-1 ) = ap( kc+k+kp-1 )
284 ap( kc+k+kp-1 ) = temp
285 END IF
286 END IF
287*
288 k = k + kstep
289 kc = kcnext
290 GO TO 30
291 50 CONTINUE
292*
293 ELSE
294*
295* Compute inv(A) from the factorization A = L*D*L**T.
296*
297* K is the main loop index, increasing from 1 to N in steps of
298* 1 or 2, depending on the size of the diagonal blocks.
299*
300 npp = n*( n+1 ) / 2
301 k = n
302 kc = npp
303 60 CONTINUE
304*
305* If K < 1, exit from loop.
306*
307 IF( k.LT.1 )
308 $ GO TO 80
309*
310 kcnext = kc - ( n-k+2 )
311 IF( ipiv( k ).GT.0 ) THEN
312*
313* 1 x 1 diagonal block
314*
315* Invert the diagonal block.
316*
317 ap( kc ) = one / ap( kc )
318*
319* Compute column K of the inverse.
320*
321 IF( k.LT.n ) THEN
322 CALL ccopy( n-k, ap( kc+1 ), 1, work, 1 )
323 CALL cspmv( uplo, n-k, -one, ap( kc+n-k+1 ), work, 1,
324 $ zero, ap( kc+1 ), 1 )
325 ap( kc ) = ap( kc ) - cdotu( n-k, work, 1, ap( kc+1 ),
326 $ 1 )
327 END IF
328 kstep = 1
329 ELSE
330*
331* 2 x 2 diagonal block
332*
333* Invert the diagonal block.
334*
335 t = ap( kcnext+1 )
336 ak = ap( kcnext ) / t
337 akp1 = ap( kc ) / t
338 akkp1 = ap( kcnext+1 ) / t
339 d = t*( ak*akp1-one )
340 ap( kcnext ) = akp1 / d
341 ap( kc ) = ak / d
342 ap( kcnext+1 ) = -akkp1 / d
343*
344* Compute columns K-1 and K of the inverse.
345*
346 IF( k.LT.n ) THEN
347 CALL ccopy( n-k, ap( kc+1 ), 1, work, 1 )
348 CALL cspmv( uplo, n-k, -one, ap( kc+( n-k+1 ) ), work, 1,
349 $ zero, ap( kc+1 ), 1 )
350 ap( kc ) = ap( kc ) - cdotu( n-k, work, 1, ap( kc+1 ),
351 $ 1 )
352 ap( kcnext+1 ) = ap( kcnext+1 ) -
353 $ cdotu( n-k, ap( kc+1 ), 1,
354 $ ap( kcnext+2 ), 1 )
355 CALL ccopy( n-k, ap( kcnext+2 ), 1, work, 1 )
356 CALL cspmv( uplo, n-k, -one, ap( kc+( n-k+1 ) ), work, 1,
357 $ zero, ap( kcnext+2 ), 1 )
358 ap( kcnext ) = ap( kcnext ) -
359 $ cdotu( n-k, work, 1, ap( kcnext+2 ), 1 )
360 END IF
361 kstep = 2
362 kcnext = kcnext - ( n-k+3 )
363 END IF
364*
365 kp = abs( ipiv( k ) )
366 IF( kp.NE.k ) THEN
367*
368* Interchange rows and columns K and KP in the trailing
369* submatrix A(k-1:n,k-1:n)
370*
371 kpc = npp - ( n-kp+1 )*( n-kp+2 ) / 2 + 1
372 IF( kp.LT.n )
373 $ CALL cswap( n-kp, ap( kc+kp-k+1 ), 1, ap( kpc+1 ), 1 )
374 kx = kc + kp - k
375 DO 70 j = k + 1, kp - 1
376 kx = kx + n - j + 1
377 temp = ap( kc+j-k )
378 ap( kc+j-k ) = ap( kx )
379 ap( kx ) = temp
380 70 CONTINUE
381 temp = ap( kc )
382 ap( kc ) = ap( kpc )
383 ap( kpc ) = temp
384 IF( kstep.EQ.2 ) THEN
385 temp = ap( kc-n+k-1 )
386 ap( kc-n+k-1 ) = ap( kc-n+kp-1 )
387 ap( kc-n+kp-1 ) = temp
388 END IF
389 END IF
390*
391 k = k - kstep
392 kc = kcnext
393 GO TO 60
394 80 CONTINUE
395 END IF
396*
397 RETURN
398*
399* End of CSPTRI
400*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
complex function cdotu(n, cx, incx, cy, incy)
CDOTU
Definition cdotu.f:83
subroutine cspmv(uplo, n, alpha, ap, x, incx, beta, y, incy)
CSPMV computes a matrix-vector product for complex vectors using a complex symmetric packed matrix
Definition cspmv.f:151
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
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: