1 /* ../netlib/ssptri.f -- translated by f2c (version 20100827). You must link the resulting object file with libf2c: on Microsoft Windows system, link with libf2c.lib;
2 on Linux or Unix systems, link with .../path/to/libf2c.a -lm or, if you install libf2c.a in a standard place, with -lf2c -lm -- in that order, at the end of the command line, as in cc *.o -lf2c -lm Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., http://www.netlib.org/f2c/libf2c.zip */
3 #include "FLA_f2c.h" /* Table of constant values */
4 static integer c__1 = 1;
5 static real c_b11 = -1.f;
6 static real c_b13 = 0.f;
7 /* > \brief \b SSPTRI */
8 /* =========== DOCUMENTATION =========== */
9 /* Online html documentation available at */
10 /* http://www.netlib.org/lapack/explore-html/ */
11 /* > \htmlonly */
12 /* > Download SSPTRI + dependencies */
13 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssptri. f"> */
14 /* > [TGZ]</a> */
15 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssptri. f"> */
16 /* > [ZIP]</a> */
17 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssptri. f"> */
18 /* > [TXT]</a> */
19 /* > \endhtmlonly */
20 /* Definition: */
21 /* =========== */
22 /* SUBROUTINE SSPTRI( UPLO, N, AP, IPIV, WORK, INFO ) */
23 /* .. Scalar Arguments .. */
24 /* CHARACTER UPLO */
25 /* INTEGER INFO, N */
26 /* .. */
27 /* .. Array Arguments .. */
28 /* INTEGER IPIV( * ) */
29 /* REAL AP( * ), WORK( * ) */
30 /* .. */
31 /* > \par Purpose: */
32 /* ============= */
33 /* > */
34 /* > \verbatim */
35 /* > */
36 /* > SSPTRI computes the inverse of a real symmetric indefinite matrix */
37 /* > A in packed storage using the factorization A = U*D*U**T or */
38 /* > A = L*D*L**T computed by SSPTRF. */
39 /* > \endverbatim */
40 /* Arguments: */
41 /* ========== */
42 /* > \param[in] UPLO */
43 /* > \verbatim */
44 /* > UPLO is CHARACTER*1 */
45 /* > Specifies whether the details of the factorization are stored */
46 /* > as an upper or lower triangular matrix. */
47 /* > = 'U': Upper triangular, form is A = U*D*U**T;
48 */
49 /* > = 'L': Lower triangular, form is A = L*D*L**T. */
50 /* > \endverbatim */
51 /* > */
52 /* > \param[in] N */
53 /* > \verbatim */
54 /* > N is INTEGER */
55 /* > The order of the matrix A. N >= 0. */
56 /* > \endverbatim */
57 /* > */
58 /* > \param[in,out] AP */
59 /* > \verbatim */
60 /* > AP is REAL array, dimension (N*(N+1)/2) */
61 /* > On entry, the block diagonal matrix D and the multipliers */
62 /* > used to obtain the factor U or L as computed by SSPTRF, */
63 /* > stored as a packed triangular matrix. */
64 /* > */
65 /* > On exit, if INFO = 0, the (symmetric) inverse of the original */
66 /* > matrix, stored as a packed triangular matrix. The j-th column */
67 /* > of inv(A) is stored in the array AP as follows: */
68 /* > if UPLO = 'U', AP(i + (j-1)*j/2) = inv(A)(i,j) for 1<=i<=j;
69 */
70 /* > if UPLO = 'L', */
71 /* > AP(i + (j-1)*(2n-j)/2) = inv(A)(i,j) for j<=i<=n. */
72 /* > \endverbatim */
73 /* > */
74 /* > \param[in] IPIV */
75 /* > \verbatim */
76 /* > IPIV is INTEGER array, dimension (N) */
77 /* > Details of the interchanges and the block structure of D */
78 /* > as determined by SSPTRF. */
79 /* > \endverbatim */
80 /* > */
81 /* > \param[out] WORK */
82 /* > \verbatim */
83 /* > WORK is REAL array, dimension (N) */
84 /* > \endverbatim */
85 /* > */
86 /* > \param[out] INFO */
87 /* > \verbatim */
88 /* > INFO is INTEGER */
89 /* > = 0: successful exit */
90 /* > < 0: if INFO = -i, the i-th argument had an illegal value */
91 /* > > 0: if INFO = i, D(i,i) = 0;
92 the matrix is singular and its */
93 /* > inverse could not be computed. */
94 /* > \endverbatim */
95 /* Authors: */
96 /* ======== */
97 /* > \author Univ. of Tennessee */
98 /* > \author Univ. of California Berkeley */
99 /* > \author Univ. of Colorado Denver */
100 /* > \author NAG Ltd. */
101 /* > \date November 2011 */
102 /* > \ingroup realOTHERcomputational */
103 /* ===================================================================== */
104 /* Subroutine */
ssptri_(char * uplo,integer * n,real * ap,integer * ipiv,real * work,integer * info)105 int ssptri_(char *uplo, integer *n, real *ap, integer *ipiv, real *work, integer *info)
106 {
107 /* System generated locals */
108 integer i__1;
109 real r__1;
110 /* Local variables */
111 real d__;
112 integer j, k;
113 real t, ak;
114 integer kc, kp, kx, kpc, npp;
115 real akp1, temp;
116 extern real sdot_(integer *, real *, integer *, real *, integer *);
117 real akkp1;
118 extern logical lsame_(char *, char *);
119 integer kstep;
120 logical upper;
121 extern /* Subroutine */
122 int scopy_(integer *, real *, integer *, real *, integer *), sswap_(integer *, real *, integer *, real *, integer * ), sspmv_(char *, integer *, real *, real *, real *, integer *, real *, real *, integer *), xerbla_(char *, integer *);
123 integer kcnext;
124 /* -- LAPACK computational routine (version 3.4.0) -- */
125 /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
126 /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
127 /* November 2011 */
128 /* .. Scalar Arguments .. */
129 /* .. */
130 /* .. Array Arguments .. */
131 /* .. */
132 /* ===================================================================== */
133 /* .. Parameters .. */
134 /* .. */
135 /* .. Local Scalars .. */
136 /* .. */
137 /* .. External Functions .. */
138 /* .. */
139 /* .. External Subroutines .. */
140 /* .. */
141 /* .. Intrinsic Functions .. */
142 /* .. */
143 /* .. Executable Statements .. */
144 /* Test the input parameters. */
145 /* Parameter adjustments */
146 --work;
147 --ipiv;
148 --ap;
149 /* Function Body */
150 *info = 0;
151 upper = lsame_(uplo, "U");
152 if (! upper && ! lsame_(uplo, "L"))
153 {
154 *info = -1;
155 }
156 else if (*n < 0)
157 {
158 *info = -2;
159 }
160 if (*info != 0)
161 {
162 i__1 = -(*info);
163 xerbla_("SSPTRI", &i__1);
164 return 0;
165 }
166 /* Quick return if possible */
167 if (*n == 0)
168 {
169 return 0;
170 }
171 /* Check that the diagonal matrix D is nonsingular. */
172 if (upper)
173 {
174 /* Upper triangular storage: examine D from bottom to top */
175 kp = *n * (*n + 1) / 2;
176 for (*info = *n;
177 *info >= 1;
178 --(*info))
179 {
180 if (ipiv[*info] > 0 && ap[kp] == 0.f)
181 {
182 return 0;
183 }
184 kp -= *info;
185 /* L10: */
186 }
187 }
188 else
189 {
190 /* Lower triangular storage: examine D from top to bottom. */
191 kp = 1;
192 i__1 = *n;
193 for (*info = 1;
194 *info <= i__1;
195 ++(*info))
196 {
197 if (ipiv[*info] > 0 && ap[kp] == 0.f)
198 {
199 return 0;
200 }
201 kp = kp + *n - *info + 1;
202 /* L20: */
203 }
204 }
205 *info = 0;
206 if (upper)
207 {
208 /* Compute inv(A) from the factorization A = U*D*U**T. */
209 /* K is the main loop index, increasing from 1 to N in steps of */
210 /* 1 or 2, depending on the size of the diagonal blocks. */
211 k = 1;
212 kc = 1;
213 L30: /* If K > N, exit from loop. */
214 if (k > *n)
215 {
216 goto L50;
217 }
218 kcnext = kc + k;
219 if (ipiv[k] > 0)
220 {
221 /* 1 x 1 diagonal block */
222 /* Invert the diagonal block. */
223 ap[kc + k - 1] = 1.f / ap[kc + k - 1];
224 /* Compute column K of the inverse. */
225 if (k > 1)
226 {
227 i__1 = k - 1;
228 scopy_(&i__1, &ap[kc], &c__1, &work[1], &c__1);
229 i__1 = k - 1;
230 sspmv_(uplo, &i__1, &c_b11, &ap[1], &work[1], &c__1, &c_b13, & ap[kc], &c__1);
231 i__1 = k - 1;
232 ap[kc + k - 1] -= sdot_(&i__1, &work[1], &c__1, &ap[kc], & c__1);
233 }
234 kstep = 1;
235 }
236 else
237 {
238 /* 2 x 2 diagonal block */
239 /* Invert the diagonal block. */
240 t = (r__1 = ap[kcnext + k - 1], f2c_abs(r__1));
241 ak = ap[kc + k - 1] / t;
242 akp1 = ap[kcnext + k] / t;
243 akkp1 = ap[kcnext + k - 1] / t;
244 d__ = t * (ak * akp1 - 1.f);
245 ap[kc + k - 1] = akp1 / d__;
246 ap[kcnext + k] = ak / d__;
247 ap[kcnext + k - 1] = -akkp1 / d__;
248 /* Compute columns K and K+1 of the inverse. */
249 if (k > 1)
250 {
251 i__1 = k - 1;
252 scopy_(&i__1, &ap[kc], &c__1, &work[1], &c__1);
253 i__1 = k - 1;
254 sspmv_(uplo, &i__1, &c_b11, &ap[1], &work[1], &c__1, &c_b13, & ap[kc], &c__1);
255 i__1 = k - 1;
256 ap[kc + k - 1] -= sdot_(&i__1, &work[1], &c__1, &ap[kc], & c__1);
257 i__1 = k - 1;
258 ap[kcnext + k - 1] -= sdot_(&i__1, &ap[kc], &c__1, &ap[kcnext] , &c__1);
259 i__1 = k - 1;
260 scopy_(&i__1, &ap[kcnext], &c__1, &work[1], &c__1);
261 i__1 = k - 1;
262 sspmv_(uplo, &i__1, &c_b11, &ap[1], &work[1], &c__1, &c_b13, & ap[kcnext], &c__1);
263 i__1 = k - 1;
264 ap[kcnext + k] -= sdot_(&i__1, &work[1], &c__1, &ap[kcnext], & c__1);
265 }
266 kstep = 2;
267 kcnext = kcnext + k + 1;
268 }
269 kp = (i__1 = ipiv[k], f2c_abs(i__1));
270 if (kp != k)
271 {
272 /* Interchange rows and columns K and KP in the leading */
273 /* submatrix A(1:k+1,1:k+1) */
274 kpc = (kp - 1) * kp / 2 + 1;
275 i__1 = kp - 1;
276 sswap_(&i__1, &ap[kc], &c__1, &ap[kpc], &c__1);
277 kx = kpc + kp - 1;
278 i__1 = k - 1;
279 for (j = kp + 1;
280 j <= i__1;
281 ++j)
282 {
283 kx = kx + j - 1;
284 temp = ap[kc + j - 1];
285 ap[kc + j - 1] = ap[kx];
286 ap[kx] = temp;
287 /* L40: */
288 }
289 temp = ap[kc + k - 1];
290 ap[kc + k - 1] = ap[kpc + kp - 1];
291 ap[kpc + kp - 1] = temp;
292 if (kstep == 2)
293 {
294 temp = ap[kc + k + k - 1];
295 ap[kc + k + k - 1] = ap[kc + k + kp - 1];
296 ap[kc + k + kp - 1] = temp;
297 }
298 }
299 k += kstep;
300 kc = kcnext;
301 goto L30;
302 L50:
303 ;
304 }
305 else
306 {
307 /* Compute inv(A) from the factorization A = L*D*L**T. */
308 /* K is the main loop index, increasing from 1 to N in steps of */
309 /* 1 or 2, depending on the size of the diagonal blocks. */
310 npp = *n * (*n + 1) / 2;
311 k = *n;
312 kc = npp;
313 L60: /* If K < 1, exit from loop. */
314 if (k < 1)
315 {
316 goto L80;
317 }
318 kcnext = kc - (*n - k + 2);
319 if (ipiv[k] > 0)
320 {
321 /* 1 x 1 diagonal block */
322 /* Invert the diagonal block. */
323 ap[kc] = 1.f / ap[kc];
324 /* Compute column K of the inverse. */
325 if (k < *n)
326 {
327 i__1 = *n - k;
328 scopy_(&i__1, &ap[kc + 1], &c__1, &work[1], &c__1);
329 i__1 = *n - k;
330 sspmv_(uplo, &i__1, &c_b11, &ap[kc + *n - k + 1], &work[1], & c__1, &c_b13, &ap[kc + 1], &c__1);
331 i__1 = *n - k;
332 ap[kc] -= sdot_(&i__1, &work[1], &c__1, &ap[kc + 1], &c__1);
333 }
334 kstep = 1;
335 }
336 else
337 {
338 /* 2 x 2 diagonal block */
339 /* Invert the diagonal block. */
340 t = (r__1 = ap[kcnext + 1], f2c_abs(r__1));
341 ak = ap[kcnext] / t;
342 akp1 = ap[kc] / t;
343 akkp1 = ap[kcnext + 1] / t;
344 d__ = t * (ak * akp1 - 1.f);
345 ap[kcnext] = akp1 / d__;
346 ap[kc] = ak / d__;
347 ap[kcnext + 1] = -akkp1 / d__;
348 /* Compute columns K-1 and K of the inverse. */
349 if (k < *n)
350 {
351 i__1 = *n - k;
352 scopy_(&i__1, &ap[kc + 1], &c__1, &work[1], &c__1);
353 i__1 = *n - k;
354 sspmv_(uplo, &i__1, &c_b11, &ap[kc + (*n - k + 1)], &work[1], &c__1, &c_b13, &ap[kc + 1], &c__1);
355 i__1 = *n - k;
356 ap[kc] -= sdot_(&i__1, &work[1], &c__1, &ap[kc + 1], &c__1);
357 i__1 = *n - k;
358 ap[kcnext + 1] -= sdot_(&i__1, &ap[kc + 1], &c__1, &ap[kcnext + 2], &c__1);
359 i__1 = *n - k;
360 scopy_(&i__1, &ap[kcnext + 2], &c__1, &work[1], &c__1);
361 i__1 = *n - k;
362 sspmv_(uplo, &i__1, &c_b11, &ap[kc + (*n - k + 1)], &work[1], &c__1, &c_b13, &ap[kcnext + 2], &c__1);
363 i__1 = *n - k;
364 ap[kcnext] -= sdot_(&i__1, &work[1], &c__1, &ap[kcnext + 2], & c__1);
365 }
366 kstep = 2;
367 kcnext -= *n - k + 3;
368 }
369 kp = (i__1 = ipiv[k], f2c_abs(i__1));
370 if (kp != k)
371 {
372 /* Interchange rows and columns K and KP in the trailing */
373 /* submatrix A(k-1:n,k-1:n) */
374 kpc = npp - (*n - kp + 1) * (*n - kp + 2) / 2 + 1;
375 if (kp < *n)
376 {
377 i__1 = *n - kp;
378 sswap_(&i__1, &ap[kc + kp - k + 1], &c__1, &ap[kpc + 1], & c__1);
379 }
380 kx = kc + kp - k;
381 i__1 = kp - 1;
382 for (j = k + 1;
383 j <= i__1;
384 ++j)
385 {
386 kx = kx + *n - j + 1;
387 temp = ap[kc + j - k];
388 ap[kc + j - k] = ap[kx];
389 ap[kx] = temp;
390 /* L70: */
391 }
392 temp = ap[kc];
393 ap[kc] = ap[kpc];
394 ap[kpc] = temp;
395 if (kstep == 2)
396 {
397 temp = ap[kc - *n + k - 1];
398 ap[kc - *n + k - 1] = ap[kc - *n + kp - 1];
399 ap[kc - *n + kp - 1] = temp;
400 }
401 }
402 k -= kstep;
403 kc = kcnext;
404 goto L60;
405 L80:
406 ;
407 }
408 return 0;
409 /* End of SSPTRI */
410 }
411 /* ssptri_ */
412