LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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.
Date
November 2011

Definition at line 111 of file ssptri.f.

111 *
112 * -- LAPACK computational routine (version 3.4.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 * November 2011
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 *
sspmv
subroutine sspmv(UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY)
SSPMV
Definition: sspmv.f:149
sdot
real function sdot(N, SX, INCX, SY, INCY)
SDOT
Definition: sdot.f:53
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
sswap
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:53
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
scopy
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53

Here is the call graph for this function:

Here is the caller graph for this function: