LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ 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 108 of file ssptri.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  REAL AP( * ), WORK( * )
121 * ..
122 *
123 * =====================================================================
124 *
125 * .. Parameters ..
126  REAL ONE, ZERO
127  parameter( one = 1.0e+0, zero = 0.0e+0 )
128 * ..
129 * .. Local Scalars ..
130  LOGICAL UPPER
131  INTEGER J, K, KC, KCNEXT, KP, KPC, KSTEP, KX, NPP
132  REAL AK, AKKP1, AKP1, D, T, TEMP
133 * ..
134 * .. External Functions ..
135  LOGICAL LSAME
136  REAL SDOT
137  EXTERNAL lsame, sdot
138 * ..
139 * .. External Subroutines ..
140  EXTERNAL scopy, sspmv, sswap, 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( 'SSPTRI', -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 scopy( k-1, ap( kc ), 1, work, 1 )
220  CALL sspmv( uplo, k-1, -one, ap, work, 1, zero, ap( kc ),
221  $ 1 )
222  ap( kc+k-1 ) = ap( kc+k-1 ) -
223  $ sdot( 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 scopy( k-1, ap( kc ), 1, work, 1 )
245  CALL sspmv( uplo, k-1, -one, ap, work, 1, zero, 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, ap( kcnext ),
251  $ 1 )
252  CALL scopy( k-1, ap( kcnext ), 1, work, 1 )
253  CALL sspmv( uplo, k-1, -one, ap, work, 1, zero,
254  $ ap( kcnext ), 1 )
255  ap( kcnext+k ) = ap( kcnext+k ) -
256  $ sdot( 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 sswap( 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 scopy( n-k, ap( kc+1 ), 1, work, 1 )
322  CALL sspmv( uplo, n-k, -one, ap( kc+n-k+1 ), work, 1,
323  $ zero, ap( kc+1 ), 1 )
324  ap( kc ) = ap( kc ) - sdot( 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 scopy( n-k, ap( kc+1 ), 1, work, 1 )
346  CALL sspmv( uplo, n-k, -one, ap( kc+( n-k+1 ) ), work, 1,
347  $ zero, ap( kc+1 ), 1 )
348  ap( kc ) = ap( kc ) - sdot( n-k, work, 1, ap( kc+1 ), 1 )
349  ap( kcnext+1 ) = ap( kcnext+1 ) -
350  $ sdot( n-k, ap( kc+1 ), 1,
351  $ ap( kcnext+2 ), 1 )
352  CALL scopy( n-k, ap( kcnext+2 ), 1, work, 1 )
353  CALL sspmv( uplo, n-k, -one, ap( kc+( n-k+1 ) ), work, 1,
354  $ zero, ap( kcnext+2 ), 1 )
355  ap( kcnext ) = ap( kcnext ) -
356  $ sdot( 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 sswap( 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 SSPTRI
397 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:82
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
real function sdot(N, SX, INCX, SY, INCY)
SDOT
Definition: sdot.f:82
subroutine sspmv(UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY)
SSPMV
Definition: sspmv.f:147
Here is the call graph for this function:
Here is the caller graph for this function: