LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ csptri()

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.

Definition at line 108 of file csptri.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 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 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:81
complex function cdotu(N, CX, INCX, CY, INCY)
CDOTU
Definition: cdotu.f:83
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
Here is the call graph for this function:
Here is the caller graph for this function: