LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
ssptri.f
Go to the documentation of this file.
1 *> \brief \b SSPTRI
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download SSPTRI + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssptri.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssptri.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssptri.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SSPTRI( UPLO, N, AP, IPIV, WORK, INFO )
22 *
23 * .. Scalar Arguments ..
24 * CHARACTER UPLO
25 * INTEGER INFO, N
26 * ..
27 * .. Array Arguments ..
28 * INTEGER IPIV( * )
29 * REAL AP( * ), WORK( * )
30 * ..
31 *
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *>
38 *> SSPTRI computes the inverse of a real symmetric indefinite matrix
39 *> A in packed storage using the factorization A = U*D*U**T or
40 *> A = L*D*L**T computed by SSPTRF.
41 *> \endverbatim
42 *
43 * Arguments:
44 * ==========
45 *
46 *> \param[in] UPLO
47 *> \verbatim
48 *> UPLO is CHARACTER*1
49 *> Specifies whether the details of the factorization are stored
50 *> as an upper or lower triangular matrix.
51 *> = 'U': Upper triangular, form is A = U*D*U**T;
52 *> = 'L': Lower triangular, form is A = L*D*L**T.
53 *> \endverbatim
54 *>
55 *> \param[in] N
56 *> \verbatim
57 *> N is INTEGER
58 *> The order of the matrix A. N >= 0.
59 *> \endverbatim
60 *>
61 *> \param[in,out] AP
62 *> \verbatim
63 *> AP is REAL array, dimension (N*(N+1)/2)
64 *> On entry, the block diagonal matrix D and the multipliers
65 *> used to obtain the factor U or L as computed by SSPTRF,
66 *> stored as a packed triangular matrix.
67 *>
68 *> On exit, if INFO = 0, the (symmetric) inverse of the original
69 *> matrix, stored as a packed triangular matrix. The j-th column
70 *> of inv(A) is stored in the array AP as follows:
71 *> if UPLO = 'U', AP(i + (j-1)*j/2) = inv(A)(i,j) for 1<=i<=j;
72 *> if UPLO = 'L',
73 *> AP(i + (j-1)*(2n-j)/2) = inv(A)(i,j) for j<=i<=n.
74 *> \endverbatim
75 *>
76 *> \param[in] IPIV
77 *> \verbatim
78 *> IPIV is INTEGER array, dimension (N)
79 *> Details of the interchanges and the block structure of D
80 *> as determined by SSPTRF.
81 *> \endverbatim
82 *>
83 *> \param[out] WORK
84 *> \verbatim
85 *> WORK is REAL array, dimension (N)
86 *> \endverbatim
87 *>
88 *> \param[out] INFO
89 *> \verbatim
90 *> INFO is INTEGER
91 *> = 0: successful exit
92 *> < 0: if INFO = -i, the i-th argument had an illegal value
93 *> > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
94 *> inverse could not be computed.
95 *> \endverbatim
96 *
97 * Authors:
98 * ========
99 *
100 *> \author Univ. of Tennessee
101 *> \author Univ. of California Berkeley
102 *> \author Univ. of Colorado Denver
103 *> \author NAG Ltd.
104 *
105 *> \ingroup realOTHERcomputational
106 *
107 * =====================================================================
108  SUBROUTINE ssptri( UPLO, N, AP, IPIV, WORK, INFO )
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 *
398  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine ssptri(UPLO, N, AP, IPIV, WORK, INFO)
SSPTRI
Definition: ssptri.f:109
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:82
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
subroutine sspmv(UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY)
SSPMV
Definition: sspmv.f:147