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