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

◆ ssptri()

subroutine ssptri ( character uplo,
integer n,
real, dimension( * ) ap,
integer, dimension( * ) ipiv,
real, dimension( * ) work,
integer info )

SSPTRI

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

Purpose:
!>
!> SSPTRI 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 SSPTRF.
!>
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 REAL 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 SSPTRF,
!>          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 SSPTRF.
!>
[out]WORK
!>          WORK is REAL 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 106 of file ssptri.f.

107*
108* -- LAPACK computational routine --
109* -- LAPACK is a software package provided by Univ. of Tennessee, --
110* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
111*
112* .. Scalar Arguments ..
113 CHARACTER UPLO
114 INTEGER INFO, N
115* ..
116* .. Array Arguments ..
117 INTEGER IPIV( * )
118 REAL AP( * ), WORK( * )
119* ..
120*
121* =====================================================================
122*
123* .. Parameters ..
124 REAL ONE, ZERO
125 parameter( one = 1.0e+0, zero = 0.0e+0 )
126* ..
127* .. Local Scalars ..
128 LOGICAL UPPER
129 INTEGER J, K, KC, KCNEXT, KP, KPC, KSTEP, KX, NPP
130 REAL AK, AKKP1, AKP1, D, T, TEMP
131* ..
132* .. External Functions ..
133 LOGICAL LSAME
134 REAL SDOT
135 EXTERNAL lsame, sdot
136* ..
137* .. External Subroutines ..
138 EXTERNAL scopy, sspmv, sswap, xerbla
139* ..
140* .. Intrinsic Functions ..
141 INTRINSIC abs
142* ..
143* .. Executable Statements ..
144*
145* Test the input parameters.
146*
147 info = 0
148 upper = lsame( uplo, 'U' )
149 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
150 info = -1
151 ELSE IF( n.LT.0 ) THEN
152 info = -2
153 END IF
154 IF( info.NE.0 ) THEN
155 CALL xerbla( 'SSPTRI', -info )
156 RETURN
157 END IF
158*
159* Quick return if possible
160*
161 IF( n.EQ.0 )
162 $ RETURN
163*
164* Check that the diagonal matrix D is nonsingular.
165*
166 IF( upper ) THEN
167*
168* Upper triangular storage: examine D from bottom to top
169*
170 kp = n*( n+1 ) / 2
171 DO 10 info = n, 1, -1
172 IF( ipiv( info ).GT.0 .AND. ap( kp ).EQ.zero )
173 $ RETURN
174 kp = kp - info
175 10 CONTINUE
176 ELSE
177*
178* Lower triangular storage: examine D from top to bottom.
179*
180 kp = 1
181 DO 20 info = 1, n
182 IF( ipiv( info ).GT.0 .AND. ap( kp ).EQ.zero )
183 $ RETURN
184 kp = kp + n - info + 1
185 20 CONTINUE
186 END IF
187 info = 0
188*
189 IF( upper ) THEN
190*
191* Compute inv(A) from the factorization A = U*D*U**T.
192*
193* K is the main loop index, increasing from 1 to N in steps of
194* 1 or 2, depending on the size of the diagonal blocks.
195*
196 k = 1
197 kc = 1
198 30 CONTINUE
199*
200* If K > N, exit from loop.
201*
202 IF( k.GT.n )
203 $ GO TO 50
204*
205 kcnext = kc + k
206 IF( ipiv( k ).GT.0 ) THEN
207*
208* 1 x 1 diagonal block
209*
210* Invert the diagonal block.
211*
212 ap( kc+k-1 ) = one / ap( kc+k-1 )
213*
214* Compute column K of the inverse.
215*
216 IF( k.GT.1 ) THEN
217 CALL scopy( k-1, ap( kc ), 1, work, 1 )
218 CALL sspmv( uplo, k-1, -one, ap, work, 1, zero,
219 $ ap( kc ),
220 $ 1 )
221 ap( kc+k-1 ) = ap( kc+k-1 ) -
222 $ sdot( k-1, work, 1, ap( kc ), 1 )
223 END IF
224 kstep = 1
225 ELSE
226*
227* 2 x 2 diagonal block
228*
229* Invert the diagonal block.
230*
231 t = abs( ap( kcnext+k-1 ) )
232 ak = ap( kc+k-1 ) / t
233 akp1 = ap( kcnext+k ) / t
234 akkp1 = ap( kcnext+k-1 ) / t
235 d = t*( ak*akp1-one )
236 ap( kc+k-1 ) = akp1 / d
237 ap( kcnext+k ) = ak / d
238 ap( kcnext+k-1 ) = -akkp1 / d
239*
240* Compute columns K and K+1 of the inverse.
241*
242 IF( k.GT.1 ) THEN
243 CALL scopy( k-1, ap( kc ), 1, work, 1 )
244 CALL sspmv( uplo, k-1, -one, ap, work, 1, zero,
245 $ ap( kc ),
246 $ 1 )
247 ap( kc+k-1 ) = ap( kc+k-1 ) -
248 $ sdot( k-1, work, 1, ap( kc ), 1 )
249 ap( kcnext+k-1 ) = ap( kcnext+k-1 ) -
250 $ sdot( k-1, ap( kc ), 1,
251 $ ap( kcnext ),
252 $ 1 )
253 CALL scopy( k-1, ap( kcnext ), 1, work, 1 )
254 CALL sspmv( uplo, k-1, -one, ap, work, 1, zero,
255 $ ap( kcnext ), 1 )
256 ap( kcnext+k ) = ap( kcnext+k ) -
257 $ sdot( 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 sswap( 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 scopy( n-k, ap( kc+1 ), 1, work, 1 )
323 CALL sspmv( uplo, n-k, -one, ap( kc+n-k+1 ), work, 1,
324 $ zero, ap( kc+1 ), 1 )
325 ap( kc ) = ap( kc ) - sdot( 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 = abs( 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 scopy( n-k, ap( kc+1 ), 1, work, 1 )
348 CALL sspmv( uplo, n-k, -one, ap( kc+( n-k+1 ) ), work,
349 $ 1,
350 $ zero, ap( kc+1 ), 1 )
351 ap( kc ) = ap( kc ) - sdot( n-k, work, 1, ap( kc+1 ),
352 $ 1 )
353 ap( kcnext+1 ) = ap( kcnext+1 ) -
354 $ sdot( n-k, ap( kc+1 ), 1,
355 $ ap( kcnext+2 ), 1 )
356 CALL scopy( n-k, ap( kcnext+2 ), 1, work, 1 )
357 CALL sspmv( uplo, n-k, -one, ap( kc+( n-k+1 ) ), work,
358 $ 1,
359 $ zero, ap( kcnext+2 ), 1 )
360 ap( kcnext ) = ap( kcnext ) -
361 $ sdot( n-k, work, 1, ap( kcnext+2 ), 1 )
362 END IF
363 kstep = 2
364 kcnext = kcnext - ( n-k+3 )
365 END IF
366*
367 kp = abs( ipiv( k ) )
368 IF( kp.NE.k ) THEN
369*
370* Interchange rows and columns K and KP in the trailing
371* submatrix A(k-1:n,k-1:n)
372*
373 kpc = npp - ( n-kp+1 )*( n-kp+2 ) / 2 + 1
374 IF( kp.LT.n )
375 $ CALL sswap( n-kp, ap( kc+kp-k+1 ), 1, ap( kpc+1 ), 1 )
376 kx = kc + kp - k
377 DO 70 j = k + 1, kp - 1
378 kx = kx + n - j + 1
379 temp = ap( kc+j-k )
380 ap( kc+j-k ) = ap( kx )
381 ap( kx ) = temp
382 70 CONTINUE
383 temp = ap( kc )
384 ap( kc ) = ap( kpc )
385 ap( kpc ) = temp
386 IF( kstep.EQ.2 ) THEN
387 temp = ap( kc-n+k-1 )
388 ap( kc-n+k-1 ) = ap( kc-n+kp-1 )
389 ap( kc-n+kp-1 ) = temp
390 END IF
391 END IF
392*
393 k = k - kstep
394 kc = kcnext
395 GO TO 60
396 80 CONTINUE
397 END IF
398*
399 RETURN
400*
401* End of SSPTRI
402*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
sdot
real function sdot(n, sx, incx, sy, incy)
SDOT
Definition sdot.f:82
sspmv
subroutine sspmv(uplo, n, alpha, ap, x, incx, beta, y, incy)
SSPMV
Definition sspmv.f:147
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
sswap
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82
Here is the call graph for this function:
Here is the caller graph for this function: