LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
csptri.f
Go to the documentation of this file.
1*> \brief \b CSPTRI
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CSPTRI + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csptri.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csptri.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csptri.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CSPTRI( UPLO, N, AP, IPIV, WORK, INFO )
22*
23* .. Scalar Arguments ..
24* CHARACTER UPLO
25* INTEGER INFO, N
26* ..
27* .. Array Arguments ..
28* INTEGER IPIV( * )
29* COMPLEX AP( * ), WORK( * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> CSPTRI computes the inverse of a complex symmetric indefinite matrix
39*> A in packed storage using the factorization A = U*D*U**T or
40*> A = L*D*L**T computed by CSPTRF.
41*> \endverbatim
42*
43* Arguments:
44* ==========
45*
46*> \param[in] UPLO
47*> \verbatim
48*> UPLO is CHARACTER*1
49*> Specifies whether the details of the factorization are stored
50*> as an upper or lower triangular matrix.
51*> = 'U': Upper triangular, form is A = U*D*U**T;
52*> = 'L': Lower triangular, form is A = L*D*L**T.
53*> \endverbatim
54*>
55*> \param[in] N
56*> \verbatim
57*> N is INTEGER
58*> The order of the matrix A. N >= 0.
59*> \endverbatim
60*>
61*> \param[in,out] AP
62*> \verbatim
63*> AP is COMPLEX array, dimension (N*(N+1)/2)
64*> On entry, the block diagonal matrix D and the multipliers
65*> used to obtain the factor U or L as computed by CSPTRF,
66*> stored as a packed triangular matrix.
67*>
68*> On exit, if INFO = 0, the (symmetric) inverse of the original
69*> matrix, stored as a packed triangular matrix. The j-th column
70*> of inv(A) is stored in the array AP as follows:
71*> if UPLO = 'U', AP(i + (j-1)*j/2) = inv(A)(i,j) for 1<=i<=j;
72*> if UPLO = 'L',
73*> AP(i + (j-1)*(2n-j)/2) = inv(A)(i,j) for j<=i<=n.
74*> \endverbatim
75*>
76*> \param[in] IPIV
77*> \verbatim
78*> IPIV is INTEGER array, dimension (N)
79*> Details of the interchanges and the block structure of D
80*> as determined by CSPTRF.
81*> \endverbatim
82*>
83*> \param[out] WORK
84*> \verbatim
85*> WORK is COMPLEX array, dimension (N)
86*> \endverbatim
87*>
88*> \param[out] INFO
89*> \verbatim
90*> INFO is INTEGER
91*> = 0: successful exit
92*> < 0: if INFO = -i, the i-th argument had an illegal value
93*> > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
94*> inverse could not be computed.
95*> \endverbatim
96*
97* Authors:
98* ========
99*
100*> \author Univ. of Tennessee
101*> \author Univ. of California Berkeley
102*> \author Univ. of Colorado Denver
103*> \author NAG Ltd.
104*
105*> \ingroup hptri
106*
107* =====================================================================
108 SUBROUTINE csptri( UPLO, N, AP, IPIV, WORK, INFO )
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*
401 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
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
subroutine csptri(uplo, n, ap, ipiv, work, info)
CSPTRI
Definition csptri.f:109
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81