LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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: