LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ dsptri()

subroutine dsptri ( character  UPLO,
integer  N,
double precision, dimension( * )  AP,
integer, dimension( * )  IPIV,
double precision, dimension( * )  WORK,
integer  INFO 
)

DSPTRI

Download DSPTRI + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 DSPTRI computes the inverse of a real symmetric indefinite matrix
 A in packed storage using the factorization A = U*D*U**T or
 A = L*D*L**T computed by DSPTRF.
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 DOUBLE PRECISION 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 DSPTRF,
          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 DSPTRF.
[out]WORK
          WORK is DOUBLE PRECISION 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 dsptri.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  DOUBLE PRECISION AP( * ), WORK( * )
121 * ..
122 *
123 * =====================================================================
124 *
125 * .. Parameters ..
126  DOUBLE PRECISION ONE, ZERO
127  parameter( one = 1.0d+0, zero = 0.0d+0 )
128 * ..
129 * .. Local Scalars ..
130  LOGICAL UPPER
131  INTEGER J, K, KC, KCNEXT, KP, KPC, KSTEP, KX, NPP
132  DOUBLE PRECISION AK, AKKP1, AKP1, D, T, TEMP
133 * ..
134 * .. External Functions ..
135  LOGICAL LSAME
136  DOUBLE PRECISION DDOT
137  EXTERNAL lsame, ddot
138 * ..
139 * .. External Subroutines ..
140  EXTERNAL dcopy, dspmv, dswap, xerbla
141 * ..
142 * .. Intrinsic Functions ..
143  INTRINSIC abs
144 * ..
145 * .. Executable Statements ..
146 *
147 * Test the input parameters.
148 *
149  info = 0
150  upper = lsame( uplo, 'U' )
151  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
152  info = -1
153  ELSE IF( n.LT.0 ) THEN
154  info = -2
155  END IF
156  IF( info.NE.0 ) THEN
157  CALL xerbla( 'DSPTRI', -info )
158  RETURN
159  END IF
160 *
161 * Quick return if possible
162 *
163  IF( n.EQ.0 )
164  $ RETURN
165 *
166 * Check that the diagonal matrix D is nonsingular.
167 *
168  IF( upper ) THEN
169 *
170 * Upper triangular storage: examine D from bottom to top
171 *
172  kp = n*( n+1 ) / 2
173  DO 10 info = n, 1, -1
174  IF( ipiv( info ).GT.0 .AND. ap( kp ).EQ.zero )
175  $ RETURN
176  kp = kp - info
177  10 CONTINUE
178  ELSE
179 *
180 * Lower triangular storage: examine D from top to bottom.
181 *
182  kp = 1
183  DO 20 info = 1, n
184  IF( ipiv( info ).GT.0 .AND. ap( kp ).EQ.zero )
185  $ RETURN
186  kp = kp + n - info + 1
187  20 CONTINUE
188  END IF
189  info = 0
190 *
191  IF( upper ) THEN
192 *
193 * Compute inv(A) from the factorization A = U*D*U**T.
194 *
195 * K is the main loop index, increasing from 1 to N in steps of
196 * 1 or 2, depending on the size of the diagonal blocks.
197 *
198  k = 1
199  kc = 1
200  30 CONTINUE
201 *
202 * If K > N, exit from loop.
203 *
204  IF( k.GT.n )
205  $ GO TO 50
206 *
207  kcnext = kc + k
208  IF( ipiv( k ).GT.0 ) THEN
209 *
210 * 1 x 1 diagonal block
211 *
212 * Invert the diagonal block.
213 *
214  ap( kc+k-1 ) = one / ap( kc+k-1 )
215 *
216 * Compute column K of the inverse.
217 *
218  IF( k.GT.1 ) THEN
219  CALL dcopy( k-1, ap( kc ), 1, work, 1 )
220  CALL dspmv( uplo, k-1, -one, ap, work, 1, zero, ap( kc ),
221  $ 1 )
222  ap( kc+k-1 ) = ap( kc+k-1 ) -
223  $ ddot( k-1, work, 1, ap( kc ), 1 )
224  END IF
225  kstep = 1
226  ELSE
227 *
228 * 2 x 2 diagonal block
229 *
230 * Invert the diagonal block.
231 *
232  t = abs( ap( kcnext+k-1 ) )
233  ak = ap( kc+k-1 ) / t
234  akp1 = ap( kcnext+k ) / t
235  akkp1 = ap( kcnext+k-1 ) / t
236  d = t*( ak*akp1-one )
237  ap( kc+k-1 ) = akp1 / d
238  ap( kcnext+k ) = ak / d
239  ap( kcnext+k-1 ) = -akkp1 / d
240 *
241 * Compute columns K and K+1 of the inverse.
242 *
243  IF( k.GT.1 ) THEN
244  CALL dcopy( k-1, ap( kc ), 1, work, 1 )
245  CALL dspmv( uplo, k-1, -one, ap, work, 1, zero, ap( kc ),
246  $ 1 )
247  ap( kc+k-1 ) = ap( kc+k-1 ) -
248  $ ddot( k-1, work, 1, ap( kc ), 1 )
249  ap( kcnext+k-1 ) = ap( kcnext+k-1 ) -
250  $ ddot( k-1, ap( kc ), 1, ap( kcnext ),
251  $ 1 )
252  CALL dcopy( k-1, ap( kcnext ), 1, work, 1 )
253  CALL dspmv( uplo, k-1, -one, ap, work, 1, zero,
254  $ ap( kcnext ), 1 )
255  ap( kcnext+k ) = ap( kcnext+k ) -
256  $ ddot( k-1, work, 1, ap( kcnext ), 1 )
257  END IF
258  kstep = 2
259  kcnext = kcnext + k + 1
260  END IF
261 *
262  kp = abs( ipiv( k ) )
263  IF( kp.NE.k ) THEN
264 *
265 * Interchange rows and columns K and KP in the leading
266 * submatrix A(1:k+1,1:k+1)
267 *
268  kpc = ( kp-1 )*kp / 2 + 1
269  CALL dswap( kp-1, ap( kc ), 1, ap( kpc ), 1 )
270  kx = kpc + kp - 1
271  DO 40 j = kp + 1, k - 1
272  kx = kx + j - 1
273  temp = ap( kc+j-1 )
274  ap( kc+j-1 ) = ap( kx )
275  ap( kx ) = temp
276  40 CONTINUE
277  temp = ap( kc+k-1 )
278  ap( kc+k-1 ) = ap( kpc+kp-1 )
279  ap( kpc+kp-1 ) = temp
280  IF( kstep.EQ.2 ) THEN
281  temp = ap( kc+k+k-1 )
282  ap( kc+k+k-1 ) = ap( kc+k+kp-1 )
283  ap( kc+k+kp-1 ) = temp
284  END IF
285  END IF
286 *
287  k = k + kstep
288  kc = kcnext
289  GO TO 30
290  50 CONTINUE
291 *
292  ELSE
293 *
294 * Compute inv(A) from the factorization A = L*D*L**T.
295 *
296 * K is the main loop index, increasing from 1 to N in steps of
297 * 1 or 2, depending on the size of the diagonal blocks.
298 *
299  npp = n*( n+1 ) / 2
300  k = n
301  kc = npp
302  60 CONTINUE
303 *
304 * If K < 1, exit from loop.
305 *
306  IF( k.LT.1 )
307  $ GO TO 80
308 *
309  kcnext = kc - ( n-k+2 )
310  IF( ipiv( k ).GT.0 ) THEN
311 *
312 * 1 x 1 diagonal block
313 *
314 * Invert the diagonal block.
315 *
316  ap( kc ) = one / ap( kc )
317 *
318 * Compute column K of the inverse.
319 *
320  IF( k.LT.n ) THEN
321  CALL dcopy( n-k, ap( kc+1 ), 1, work, 1 )
322  CALL dspmv( uplo, n-k, -one, ap( kc+n-k+1 ), work, 1,
323  $ zero, ap( kc+1 ), 1 )
324  ap( kc ) = ap( kc ) - ddot( n-k, work, 1, ap( kc+1 ), 1 )
325  END IF
326  kstep = 1
327  ELSE
328 *
329 * 2 x 2 diagonal block
330 *
331 * Invert the diagonal block.
332 *
333  t = abs( ap( kcnext+1 ) )
334  ak = ap( kcnext ) / t
335  akp1 = ap( kc ) / t
336  akkp1 = ap( kcnext+1 ) / t
337  d = t*( ak*akp1-one )
338  ap( kcnext ) = akp1 / d
339  ap( kc ) = ak / d
340  ap( kcnext+1 ) = -akkp1 / d
341 *
342 * Compute columns K-1 and K of the inverse.
343 *
344  IF( k.LT.n ) THEN
345  CALL dcopy( n-k, ap( kc+1 ), 1, work, 1 )
346  CALL dspmv( uplo, n-k, -one, ap( kc+( n-k+1 ) ), work, 1,
347  $ zero, ap( kc+1 ), 1 )
348  ap( kc ) = ap( kc ) - ddot( n-k, work, 1, ap( kc+1 ), 1 )
349  ap( kcnext+1 ) = ap( kcnext+1 ) -
350  $ ddot( n-k, ap( kc+1 ), 1,
351  $ ap( kcnext+2 ), 1 )
352  CALL dcopy( n-k, ap( kcnext+2 ), 1, work, 1 )
353  CALL dspmv( uplo, n-k, -one, ap( kc+( n-k+1 ) ), work, 1,
354  $ zero, ap( kcnext+2 ), 1 )
355  ap( kcnext ) = ap( kcnext ) -
356  $ ddot( n-k, work, 1, ap( kcnext+2 ), 1 )
357  END IF
358  kstep = 2
359  kcnext = kcnext - ( n-k+3 )
360  END IF
361 *
362  kp = abs( ipiv( k ) )
363  IF( kp.NE.k ) THEN
364 *
365 * Interchange rows and columns K and KP in the trailing
366 * submatrix A(k-1:n,k-1:n)
367 *
368  kpc = npp - ( n-kp+1 )*( n-kp+2 ) / 2 + 1
369  IF( kp.LT.n )
370  $ CALL dswap( n-kp, ap( kc+kp-k+1 ), 1, ap( kpc+1 ), 1 )
371  kx = kc + kp - k
372  DO 70 j = k + 1, kp - 1
373  kx = kx + n - j + 1
374  temp = ap( kc+j-k )
375  ap( kc+j-k ) = ap( kx )
376  ap( kx ) = temp
377  70 CONTINUE
378  temp = ap( kc )
379  ap( kc ) = ap( kpc )
380  ap( kpc ) = temp
381  IF( kstep.EQ.2 ) THEN
382  temp = ap( kc-n+k-1 )
383  ap( kc-n+k-1 ) = ap( kc-n+kp-1 )
384  ap( kc-n+kp-1 ) = temp
385  END IF
386  END IF
387 *
388  k = k - kstep
389  kc = kcnext
390  GO TO 60
391  80 CONTINUE
392  END IF
393 *
394  RETURN
395 *
396 * End of DSPTRI
397 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
double precision function ddot(N, DX, INCX, DY, INCY)
DDOT
Definition: ddot.f:82
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:82
subroutine dspmv(UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY)
DSPMV
Definition: dspmv.f:147
Here is the call graph for this function:
Here is the caller graph for this function: