LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
zsptri.f
Go to the documentation of this file.
1*> \brief \b ZSPTRI
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsptri.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsptri.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsptri.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE ZSPTRI( 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*16 AP( * ), WORK( * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> ZSPTRI 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 ZSPTRF.
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*16 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 ZSPTRF,
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 ZSPTRF.
81*> \endverbatim
82*>
83*> \param[out] WORK
84*> \verbatim
85*> WORK is COMPLEX*16 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 complex16OTHERcomputational
106*
107* =====================================================================
108 SUBROUTINE zsptri( 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*16 AP( * ), WORK( * )
121* ..
122*
123* =====================================================================
124*
125* .. Parameters ..
126 COMPLEX*16 ONE, ZERO
127 parameter( one = ( 1.0d+0, 0.0d+0 ),
128 \$ zero = ( 0.0d+0, 0.0d+0 ) )
129* ..
130* .. Local Scalars ..
131 LOGICAL UPPER
132 INTEGER J, K, KC, KCNEXT, KP, KPC, KSTEP, KX, NPP
133 COMPLEX*16 AK, AKKP1, AKP1, D, T, TEMP
134* ..
135* .. External Functions ..
136 LOGICAL LSAME
137 COMPLEX*16 ZDOTU
138 EXTERNAL lsame, zdotu
139* ..
140* .. External Subroutines ..
141 EXTERNAL xerbla, zcopy, zspmv, zswap
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( 'ZSPTRI', -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 zcopy( k-1, ap( kc ), 1, work, 1 )
221 CALL zspmv( uplo, k-1, -one, ap, work, 1, zero, ap( kc ),
222 \$ 1 )
223 ap( kc+k-1 ) = ap( kc+k-1 ) -
224 \$ zdotu( 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 zcopy( k-1, ap( kc ), 1, work, 1 )
246 CALL zspmv( uplo, k-1, -one, ap, work, 1, zero, ap( kc ),
247 \$ 1 )
248 ap( kc+k-1 ) = ap( kc+k-1 ) -
249 \$ zdotu( k-1, work, 1, ap( kc ), 1 )
250 ap( kcnext+k-1 ) = ap( kcnext+k-1 ) -
251 \$ zdotu( k-1, ap( kc ), 1, ap( kcnext ),
252 \$ 1 )
253 CALL zcopy( k-1, ap( kcnext ), 1, work, 1 )
254 CALL zspmv( uplo, k-1, -one, ap, work, 1, zero,
255 \$ ap( kcnext ), 1 )
256 ap( kcnext+k ) = ap( kcnext+k ) -
257 \$ zdotu( 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 zswap( 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 zcopy( n-k, ap( kc+1 ), 1, work, 1 )
323 CALL zspmv( uplo, n-k, -one, ap( kc+n-k+1 ), work, 1,
324 \$ zero, ap( kc+1 ), 1 )
325 ap( kc ) = ap( kc ) - zdotu( 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 zcopy( n-k, ap( kc+1 ), 1, work, 1 )
348 CALL zspmv( uplo, n-k, -one, ap( kc+( n-k+1 ) ), work, 1,
349 \$ zero, ap( kc+1 ), 1 )
350 ap( kc ) = ap( kc ) - zdotu( n-k, work, 1, ap( kc+1 ),
351 \$ 1 )
352 ap( kcnext+1 ) = ap( kcnext+1 ) -
353 \$ zdotu( n-k, ap( kc+1 ), 1,
354 \$ ap( kcnext+2 ), 1 )
355 CALL zcopy( n-k, ap( kcnext+2 ), 1, work, 1 )
356 CALL zspmv( uplo, n-k, -one, ap( kc+( n-k+1 ) ), work, 1,
357 \$ zero, ap( kcnext+2 ), 1 )
358 ap( kcnext ) = ap( kcnext ) -
359 \$ zdotu( 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 zswap( 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 ZSPTRI
400*
401 END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:81
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:81
subroutine zspmv(UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY)
ZSPMV computes a matrix-vector product for complex vectors using a complex symmetric packed matrix
Definition: zspmv.f:151
subroutine zsptri(UPLO, N, AP, IPIV, WORK, INFO)
ZSPTRI
Definition: zsptri.f:109