LAPACK  3.10.1 LAPACK: Linear Algebra PACKage

## ◆ zsptri()

 subroutine zsptri ( character UPLO, integer N, complex*16, dimension( * ) AP, integer, dimension( * ) IPIV, complex*16, dimension( * ) WORK, integer INFO )

ZSPTRI

Purpose:
``` ZSPTRI 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 ZSPTRF.```
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*16 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 ZSPTRF, 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 ZSPTRF.``` [out] WORK ` WORK is COMPLEX*16 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.```

Definition at line 108 of file zsptri.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*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 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:81
complex *16 function zdotu(N, ZX, INCX, ZY, INCY)
ZDOTU
Definition: zdotu.f:83
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
Here is the call graph for this function:
Here is the caller graph for this function: