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