LAPACK  3.6.1 LAPACK: Linear Algebra PACKage
 subroutine csptri ( character UPLO, integer N, complex, dimension( * ) AP, integer, dimension( * ) IPIV, complex, dimension( * ) WORK, integer INFO )

CSPTRI

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.
Date
November 2011

Definition at line 111 of file csptri.f.

111 *
112 * -- LAPACK computational routine (version 3.4.0) --
113 * -- LAPACK is a software package provided by Univ. of Tennessee, --
114 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
115 * November 2011
116 *
117 * .. Scalar Arguments ..
118  CHARACTER uplo
119  INTEGER info, n
120 * ..
121 * .. Array Arguments ..
122  INTEGER ipiv( * )
123  COMPLEX ap( * ), work( * )
124 * ..
125 *
126 * =====================================================================
127 *
128 * .. Parameters ..
129  COMPLEX one, zero
130  parameter ( one = ( 1.0e+0, 0.0e+0 ),
131  \$ zero = ( 0.0e+0, 0.0e+0 ) )
132 * ..
133 * .. Local Scalars ..
134  LOGICAL upper
135  INTEGER j, k, kc, kcnext, kp, kpc, kstep, kx, npp
136  COMPLEX ak, akkp1, akp1, d, t, temp
137 * ..
138 * .. External Functions ..
139  LOGICAL lsame
140  COMPLEX cdotu
141  EXTERNAL lsame, cdotu
142 * ..
143 * .. External Subroutines ..
144  EXTERNAL ccopy, cspmv, cswap, xerbla
145 * ..
146 * .. Intrinsic Functions ..
147  INTRINSIC abs
148 * ..
149 * .. Executable Statements ..
150 *
151 * Test the input parameters.
152 *
153  info = 0
154  upper = lsame( uplo, 'U' )
155  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
156  info = -1
157  ELSE IF( n.LT.0 ) THEN
158  info = -2
159  END IF
160  IF( info.NE.0 ) THEN
161  CALL xerbla( 'CSPTRI', -info )
162  RETURN
163  END IF
164 *
165 * Quick return if possible
166 *
167  IF( n.EQ.0 )
168  \$ RETURN
169 *
170 * Check that the diagonal matrix D is nonsingular.
171 *
172  IF( upper ) THEN
173 *
174 * Upper triangular storage: examine D from bottom to top
175 *
176  kp = n*( n+1 ) / 2
177  DO 10 info = n, 1, -1
178  IF( ipiv( info ).GT.0 .AND. ap( kp ).EQ.zero )
179  \$ RETURN
180  kp = kp - info
181  10 CONTINUE
182  ELSE
183 *
184 * Lower triangular storage: examine D from top to bottom.
185 *
186  kp = 1
187  DO 20 info = 1, n
188  IF( ipiv( info ).GT.0 .AND. ap( kp ).EQ.zero )
189  \$ RETURN
190  kp = kp + n - info + 1
191  20 CONTINUE
192  END IF
193  info = 0
194 *
195  IF( upper ) THEN
196 *
197 * Compute inv(A) from the factorization A = U*D*U**T.
198 *
199 * K is the main loop index, increasing from 1 to N in steps of
200 * 1 or 2, depending on the size of the diagonal blocks.
201 *
202  k = 1
203  kc = 1
204  30 CONTINUE
205 *
206 * If K > N, exit from loop.
207 *
208  IF( k.GT.n )
209  \$ GO TO 50
210 *
211  kcnext = kc + k
212  IF( ipiv( k ).GT.0 ) THEN
213 *
214 * 1 x 1 diagonal block
215 *
216 * Invert the diagonal block.
217 *
218  ap( kc+k-1 ) = one / ap( kc+k-1 )
219 *
220 * Compute column K of the inverse.
221 *
222  IF( k.GT.1 ) THEN
223  CALL ccopy( k-1, ap( kc ), 1, work, 1 )
224  CALL cspmv( uplo, k-1, -one, ap, work, 1, zero, ap( kc ),
225  \$ 1 )
226  ap( kc+k-1 ) = ap( kc+k-1 ) -
227  \$ cdotu( k-1, work, 1, ap( kc ), 1 )
228  END IF
229  kstep = 1
230  ELSE
231 *
232 * 2 x 2 diagonal block
233 *
234 * Invert the diagonal block.
235 *
236  t = ap( kcnext+k-1 )
237  ak = ap( kc+k-1 ) / t
238  akp1 = ap( kcnext+k ) / t
239  akkp1 = ap( kcnext+k-1 ) / t
240  d = t*( ak*akp1-one )
241  ap( kc+k-1 ) = akp1 / d
242  ap( kcnext+k ) = ak / d
243  ap( kcnext+k-1 ) = -akkp1 / d
244 *
245 * Compute columns K and K+1 of the inverse.
246 *
247  IF( k.GT.1 ) THEN
248  CALL ccopy( k-1, ap( kc ), 1, work, 1 )
249  CALL cspmv( uplo, k-1, -one, ap, work, 1, zero, ap( kc ),
250  \$ 1 )
251  ap( kc+k-1 ) = ap( kc+k-1 ) -
252  \$ cdotu( k-1, work, 1, ap( kc ), 1 )
253  ap( kcnext+k-1 ) = ap( kcnext+k-1 ) -
254  \$ cdotu( k-1, ap( kc ), 1, ap( kcnext ),
255  \$ 1 )
256  CALL ccopy( k-1, ap( kcnext ), 1, work, 1 )
257  CALL cspmv( uplo, k-1, -one, ap, work, 1, zero,
258  \$ ap( kcnext ), 1 )
259  ap( kcnext+k ) = ap( kcnext+k ) -
260  \$ cdotu( k-1, work, 1, ap( kcnext ), 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 = ap( kc+j-1 )
278  ap( kc+j-1 ) = ap( kx )
279  ap( kx ) = temp
280  40 CONTINUE
281  temp = ap( kc+k-1 )
282  ap( kc+k-1 ) = ap( kpc+kp-1 )
283  ap( kpc+kp-1 ) = temp
284  IF( kstep.EQ.2 ) THEN
285  temp = ap( kc+k+k-1 )
286  ap( kc+k+k-1 ) = ap( kc+k+kp-1 )
287  ap( kc+k+kp-1 ) = temp
288  END IF
289  END IF
290 *
291  k = k + kstep
292  kc = kcnext
293  GO TO 30
294  50 CONTINUE
295 *
296  ELSE
297 *
298 * Compute inv(A) from the factorization A = L*D*L**T.
299 *
300 * K is the main loop index, increasing from 1 to N in steps of
301 * 1 or 2, depending on the size of the diagonal blocks.
302 *
303  npp = n*( n+1 ) / 2
304  k = n
305  kc = npp
306  60 CONTINUE
307 *
308 * If K < 1, exit from loop.
309 *
310  IF( k.LT.1 )
311  \$ GO TO 80
312 *
313  kcnext = kc - ( n-k+2 )
314  IF( ipiv( k ).GT.0 ) THEN
315 *
316 * 1 x 1 diagonal block
317 *
318 * Invert the diagonal block.
319 *
320  ap( kc ) = one / ap( kc )
321 *
322 * Compute column K of the inverse.
323 *
324  IF( k.LT.n ) THEN
325  CALL ccopy( n-k, ap( kc+1 ), 1, work, 1 )
326  CALL cspmv( uplo, n-k, -one, ap( kc+n-k+1 ), work, 1,
327  \$ zero, ap( kc+1 ), 1 )
328  ap( kc ) = ap( kc ) - cdotu( n-k, work, 1, ap( kc+1 ),
329  \$ 1 )
330  END IF
331  kstep = 1
332  ELSE
333 *
334 * 2 x 2 diagonal block
335 *
336 * Invert the diagonal block.
337 *
338  t = ap( kcnext+1 )
339  ak = ap( kcnext ) / t
340  akp1 = ap( kc ) / t
341  akkp1 = ap( kcnext+1 ) / t
342  d = t*( ak*akp1-one )
343  ap( kcnext ) = akp1 / d
344  ap( kc ) = ak / d
345  ap( kcnext+1 ) = -akkp1 / d
346 *
347 * Compute columns K-1 and K of the inverse.
348 *
349  IF( k.LT.n ) THEN
350  CALL ccopy( n-k, ap( kc+1 ), 1, work, 1 )
351  CALL cspmv( uplo, n-k, -one, ap( kc+( n-k+1 ) ), work, 1,
352  \$ zero, ap( kc+1 ), 1 )
353  ap( kc ) = ap( kc ) - cdotu( n-k, work, 1, ap( kc+1 ),
354  \$ 1 )
355  ap( kcnext+1 ) = ap( kcnext+1 ) -
356  \$ cdotu( n-k, ap( kc+1 ), 1,
357  \$ ap( kcnext+2 ), 1 )
358  CALL ccopy( n-k, ap( kcnext+2 ), 1, work, 1 )
359  CALL cspmv( uplo, n-k, -one, ap( kc+( n-k+1 ) ), work, 1,
360  \$ zero, ap( kcnext+2 ), 1 )
361  ap( kcnext ) = ap( kcnext ) -
362  \$ cdotu( n-k, work, 1, ap( kcnext+2 ), 1 )
363  END IF
364  kstep = 2
365  kcnext = kcnext - ( n-k+3 )
366  END IF
367 *
368  kp = abs( ipiv( k ) )
369  IF( kp.NE.k ) THEN
370 *
371 * Interchange rows and columns K and KP in the trailing
372 * submatrix A(k-1:n,k-1:n)
373 *
374  kpc = npp - ( n-kp+1 )*( n-kp+2 ) / 2 + 1
375  IF( kp.LT.n )
376  \$ CALL cswap( n-kp, ap( kc+kp-k+1 ), 1, ap( kpc+1 ), 1 )
377  kx = kc + kp - k
378  DO 70 j = k + 1, kp - 1
379  kx = kx + n - j + 1
380  temp = ap( kc+j-k )
381  ap( kc+j-k ) = ap( kx )
382  ap( kx ) = temp
383  70 CONTINUE
384  temp = ap( kc )
385  ap( kc ) = ap( kpc )
386  ap( kpc ) = temp
387  IF( kstep.EQ.2 ) THEN
388  temp = ap( kc-n+k-1 )
389  ap( kc-n+k-1 ) = ap( kc-n+kp-1 )
390  ap( kc-n+kp-1 ) = temp
391  END IF
392  END IF
393 *
394  k = k - kstep
395  kc = kcnext
396  GO TO 60
397  80 CONTINUE
398  END IF
399 *
400  RETURN
401 *
402 * End of CSPTRI
403 *
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:153
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
complex function cdotu(N, CX, INCX, CY, INCY)
CDOTU
Definition: cdotu.f:54
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:52
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:52
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function: