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