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

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.
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 *
cspmv
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
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
cdotu
complex function cdotu(N, CX, INCX, CY, INCY)
CDOTU
Definition: cdotu.f:54
ccopy
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:52
cswap
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:52
lsame
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: