LAPACK  3.8.0 LAPACK: Linear Algebra PACKage

## ◆ ssptri()

 subroutine ssptri ( character UPLO, integer N, real, dimension( * ) AP, integer, dimension( * ) IPIV, real, dimension( * ) WORK, integer INFO )

SSPTRI

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.```
Date
December 2016

Definition at line 111 of file ssptri.f.

111 *
112 * -- LAPACK computational routine (version 3.7.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 * December 2016
116 *
117 * .. Scalar Arguments ..
118  CHARACTER uplo
119  INTEGER info, n
120 * ..
121 * .. Array Arguments ..
122  INTEGER ipiv( * )
123  REAL ap( * ), work( * )
124 * ..
125 *
126 * =====================================================================
127 *
128 * .. Parameters ..
129  REAL one, zero
130  parameter( one = 1.0e+0, zero = 0.0e+0 )
131 * ..
132 * .. Local Scalars ..
133  LOGICAL upper
134  INTEGER j, k, kc, kcnext, kp, kpc, kstep, kx, npp
135  REAL ak, akkp1, akp1, d, t, temp
136 * ..
137 * .. External Functions ..
138  LOGICAL lsame
139  REAL sdot
140  EXTERNAL lsame, sdot
141 * ..
142 * .. External Subroutines ..
143  EXTERNAL scopy, sspmv, sswap, xerbla
144 * ..
145 * .. Intrinsic Functions ..
146  INTRINSIC abs
147 * ..
148 * .. Executable Statements ..
149 *
150 * Test the input parameters.
151 *
152  info = 0
153  upper = lsame( uplo, 'U' )
154  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
155  info = -1
156  ELSE IF( n.LT.0 ) THEN
157  info = -2
158  END IF
159  IF( info.NE.0 ) THEN
160  CALL xerbla( 'SSPTRI', -info )
161  RETURN
162  END IF
163 *
164 * Quick return if possible
165 *
166  IF( n.EQ.0 )
167  \$ RETURN
168 *
169 * Check that the diagonal matrix D is nonsingular.
170 *
171  IF( upper ) THEN
172 *
173 * Upper triangular storage: examine D from bottom to top
174 *
175  kp = n*( n+1 ) / 2
176  DO 10 info = n, 1, -1
177  IF( ipiv( info ).GT.0 .AND. ap( kp ).EQ.zero )
178  \$ RETURN
179  kp = kp - info
180  10 CONTINUE
181  ELSE
182 *
183 * Lower triangular storage: examine D from top to bottom.
184 *
185  kp = 1
186  DO 20 info = 1, n
187  IF( ipiv( info ).GT.0 .AND. ap( kp ).EQ.zero )
188  \$ RETURN
189  kp = kp + n - info + 1
190  20 CONTINUE
191  END IF
192  info = 0
193 *
194  IF( upper ) THEN
195 *
196 * Compute inv(A) from the factorization A = U*D*U**T.
197 *
198 * K is the main loop index, increasing from 1 to N in steps of
199 * 1 or 2, depending on the size of the diagonal blocks.
200 *
201  k = 1
202  kc = 1
203  30 CONTINUE
204 *
205 * If K > N, exit from loop.
206 *
207  IF( k.GT.n )
208  \$ GO TO 50
209 *
210  kcnext = kc + k
211  IF( ipiv( k ).GT.0 ) THEN
212 *
213 * 1 x 1 diagonal block
214 *
215 * Invert the diagonal block.
216 *
217  ap( kc+k-1 ) = one / ap( kc+k-1 )
218 *
219 * Compute column K of the inverse.
220 *
221  IF( k.GT.1 ) THEN
222  CALL scopy( k-1, ap( kc ), 1, work, 1 )
223  CALL sspmv( uplo, k-1, -one, ap, work, 1, zero, ap( kc ),
224  \$ 1 )
225  ap( kc+k-1 ) = ap( kc+k-1 ) -
226  \$ sdot( k-1, work, 1, ap( kc ), 1 )
227  END IF
228  kstep = 1
229  ELSE
230 *
231 * 2 x 2 diagonal block
232 *
233 * Invert the diagonal block.
234 *
235  t = abs( ap( kcnext+k-1 ) )
236  ak = ap( kc+k-1 ) / t
237  akp1 = ap( kcnext+k ) / t
238  akkp1 = ap( kcnext+k-1 ) / t
239  d = t*( ak*akp1-one )
240  ap( kc+k-1 ) = akp1 / d
241  ap( kcnext+k ) = ak / d
242  ap( kcnext+k-1 ) = -akkp1 / d
243 *
244 * Compute columns K and K+1 of the inverse.
245 *
246  IF( k.GT.1 ) THEN
247  CALL scopy( k-1, ap( kc ), 1, work, 1 )
248  CALL sspmv( uplo, k-1, -one, ap, work, 1, zero, ap( kc ),
249  \$ 1 )
250  ap( kc+k-1 ) = ap( kc+k-1 ) -
251  \$ sdot( k-1, work, 1, ap( kc ), 1 )
252  ap( kcnext+k-1 ) = ap( kcnext+k-1 ) -
253  \$ sdot( k-1, ap( kc ), 1, ap( kcnext ),
254  \$ 1 )
255  CALL scopy( k-1, ap( kcnext ), 1, work, 1 )
256  CALL sspmv( uplo, k-1, -one, ap, work, 1, zero,
257  \$ ap( kcnext ), 1 )
258  ap( kcnext+k ) = ap( kcnext+k ) -
259  \$ sdot( k-1, work, 1, ap( kcnext ), 1 )
260  END IF
261  kstep = 2
262  kcnext = kcnext + k + 1
263  END IF
264 *
265  kp = abs( ipiv( k ) )
266  IF( kp.NE.k ) THEN
267 *
268 * Interchange rows and columns K and KP in the leading
269 * submatrix A(1:k+1,1:k+1)
270 *
271  kpc = ( kp-1 )*kp / 2 + 1
272  CALL sswap( kp-1, ap( kc ), 1, ap( kpc ), 1 )
273  kx = kpc + kp - 1
274  DO 40 j = kp + 1, k - 1
275  kx = kx + j - 1
276  temp = ap( kc+j-1 )
277  ap( kc+j-1 ) = ap( kx )
278  ap( kx ) = temp
279  40 CONTINUE
280  temp = ap( kc+k-1 )
281  ap( kc+k-1 ) = ap( kpc+kp-1 )
282  ap( kpc+kp-1 ) = temp
283  IF( kstep.EQ.2 ) THEN
284  temp = ap( kc+k+k-1 )
285  ap( kc+k+k-1 ) = ap( kc+k+kp-1 )
286  ap( kc+k+kp-1 ) = temp
287  END IF
288  END IF
289 *
290  k = k + kstep
291  kc = kcnext
292  GO TO 30
293  50 CONTINUE
294 *
295  ELSE
296 *
297 * Compute inv(A) from the factorization A = L*D*L**T.
298 *
299 * K is the main loop index, increasing from 1 to N in steps of
300 * 1 or 2, depending on the size of the diagonal blocks.
301 *
302  npp = n*( n+1 ) / 2
303  k = n
304  kc = npp
305  60 CONTINUE
306 *
307 * If K < 1, exit from loop.
308 *
309  IF( k.LT.1 )
310  \$ GO TO 80
311 *
312  kcnext = kc - ( n-k+2 )
313  IF( ipiv( k ).GT.0 ) THEN
314 *
315 * 1 x 1 diagonal block
316 *
317 * Invert the diagonal block.
318 *
319  ap( kc ) = one / ap( kc )
320 *
321 * Compute column K of the inverse.
322 *
323  IF( k.LT.n ) THEN
324  CALL scopy( n-k, ap( kc+1 ), 1, work, 1 )
325  CALL sspmv( uplo, n-k, -one, ap( kc+n-k+1 ), work, 1,
326  \$ zero, ap( kc+1 ), 1 )
327  ap( kc ) = ap( kc ) - sdot( n-k, work, 1, ap( kc+1 ), 1 )
328  END IF
329  kstep = 1
330  ELSE
331 *
332 * 2 x 2 diagonal block
333 *
334 * Invert the diagonal block.
335 *
336  t = abs( ap( kcnext+1 ) )
337  ak = ap( kcnext ) / t
338  akp1 = ap( kc ) / t
339  akkp1 = ap( kcnext+1 ) / t
340  d = t*( ak*akp1-one )
341  ap( kcnext ) = akp1 / d
342  ap( kc ) = ak / d
343  ap( kcnext+1 ) = -akkp1 / d
344 *
345 * Compute columns K-1 and K of the inverse.
346 *
347  IF( k.LT.n ) THEN
348  CALL scopy( n-k, ap( kc+1 ), 1, work, 1 )
349  CALL sspmv( uplo, n-k, -one, ap( kc+( n-k+1 ) ), work, 1,
350  \$ zero, ap( kc+1 ), 1 )
351  ap( kc ) = ap( kc ) - sdot( n-k, work, 1, ap( kc+1 ), 1 )
352  ap( kcnext+1 ) = ap( kcnext+1 ) -
353  \$ sdot( n-k, ap( kc+1 ), 1,
354  \$ ap( kcnext+2 ), 1 )
355  CALL scopy( n-k, ap( kcnext+2 ), 1, work, 1 )
356  CALL sspmv( uplo, n-k, -one, ap( kc+( n-k+1 ) ), work, 1,
357  \$ zero, ap( kcnext+2 ), 1 )
358  ap( kcnext ) = ap( kcnext ) -
359  \$ sdot( 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 sswap( 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 SSPTRI
400 *
subroutine sspmv(UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY)
SSPMV
Definition: sspmv.f:149
real function sdot(N, SX, INCX, SY, INCY)
SDOT
Definition: sdot.f:84
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:84
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:84
Here is the call graph for this function:
Here is the caller graph for this function: