1 /* ../netlib/dsptrs.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 doublereal c_b7 = -1.;
5 static integer c__1 = 1;
6 static doublereal c_b19 = 1.;
7 /* > \brief \b DSPTRS */
8 /* =========== DOCUMENTATION =========== */
9 /* Online html documentation available at */
10 /* http://www.netlib.org/lapack/explore-html/ */
11 /* > \htmlonly */
12 /* > Download DSPTRS + dependencies */
13 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsptrs. f"> */
14 /* > [TGZ]</a> */
15 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsptrs. f"> */
16 /* > [ZIP]</a> */
17 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsptrs. f"> */
18 /* > [TXT]</a> */
19 /* > \endhtmlonly */
20 /* Definition: */
21 /* =========== */
22 /* SUBROUTINE DSPTRS( UPLO, N, NRHS, AP, IPIV, B, LDB, INFO ) */
23 /* .. Scalar Arguments .. */
24 /* CHARACTER UPLO */
25 /* INTEGER INFO, LDB, N, NRHS */
26 /* .. */
27 /* .. Array Arguments .. */
28 /* INTEGER IPIV( * ) */
29 /* DOUBLE PRECISION AP( * ), B( LDB, * ) */
30 /* .. */
31 /* > \par Purpose: */
32 /* ============= */
33 /* > */
34 /* > \verbatim */
35 /* > */
36 /* > DSPTRS solves a system of linear equations A*X = B with a real */
37 /* > symmetric matrix A stored in packed format using the factorization */
38 /* > A = U*D*U**T or A = L*D*L**T computed by DSPTRF. */
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] NRHS */
59 /* > \verbatim */
60 /* > NRHS is INTEGER */
61 /* > The number of right hand sides, i.e., the number of columns */
62 /* > of the matrix B. NRHS >= 0. */
63 /* > \endverbatim */
64 /* > */
65 /* > \param[in] AP */
66 /* > \verbatim */
67 /* > AP is DOUBLE PRECISION array, dimension (N*(N+1)/2) */
68 /* > The block diagonal matrix D and the multipliers used to */
69 /* > obtain the factor U or L as computed by DSPTRF, stored as a */
70 /* > packed triangular matrix. */
71 /* > \endverbatim */
72 /* > */
73 /* > \param[in] IPIV */
74 /* > \verbatim */
75 /* > IPIV is INTEGER array, dimension (N) */
76 /* > Details of the interchanges and the block structure of D */
77 /* > as determined by DSPTRF. */
78 /* > \endverbatim */
79 /* > */
80 /* > \param[in,out] B */
81 /* > \verbatim */
82 /* > B is DOUBLE PRECISION array, dimension (LDB,NRHS) */
83 /* > On entry, the right hand side matrix B. */
84 /* > On exit, the solution matrix X. */
85 /* > \endverbatim */
86 /* > */
87 /* > \param[in] LDB */
88 /* > \verbatim */
89 /* > LDB is INTEGER */
90 /* > The leading dimension of the array B. LDB >= max(1,N). */
91 /* > \endverbatim */
92 /* > */
93 /* > \param[out] INFO */
94 /* > \verbatim */
95 /* > INFO is INTEGER */
96 /* > = 0: successful exit */
97 /* > < 0: if INFO = -i, the i-th argument had an illegal value */
98 /* > \endverbatim */
99 /* Authors: */
100 /* ======== */
101 /* > \author Univ. of Tennessee */
102 /* > \author Univ. of California Berkeley */
103 /* > \author Univ. of Colorado Denver */
104 /* > \author NAG Ltd. */
105 /* > \date November 2011 */
106 /* > \ingroup doubleOTHERcomputational */
107 /* ===================================================================== */
108 /* Subroutine */
dsptrs_(char * uplo,integer * n,integer * nrhs,doublereal * ap,integer * ipiv,doublereal * b,integer * ldb,integer * info)109 int dsptrs_(char *uplo, integer *n, integer *nrhs, doublereal *ap, integer *ipiv, doublereal *b, integer *ldb, integer * info)
110 {
111 /* System generated locals */
112 integer b_dim1, b_offset, i__1;
113 doublereal d__1;
114 /* Local variables */
115 integer j, k;
116 doublereal ak, bk;
117 integer kc, kp;
118 doublereal akm1, bkm1;
119 extern /* Subroutine */
120 int dger_(integer *, integer *, doublereal *, doublereal *, integer *, doublereal *, integer *, doublereal *, integer *);
121 doublereal akm1k;
122 extern /* Subroutine */
123 int dscal_(integer *, doublereal *, doublereal *, integer *);
124 extern logical lsame_(char *, char *);
125 doublereal denom;
126 extern /* Subroutine */
127 int dgemv_(char *, integer *, integer *, doublereal *, doublereal *, integer *, doublereal *, integer *, doublereal *, doublereal *, integer *), dswap_(integer *, doublereal *, integer *, doublereal *, integer *);
128 logical upper;
129 extern /* Subroutine */
130 int xerbla_(char *, integer *);
131 /* -- LAPACK computational routine (version 3.4.0) -- */
132 /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
133 /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
134 /* November 2011 */
135 /* .. Scalar Arguments .. */
136 /* .. */
137 /* .. Array Arguments .. */
138 /* .. */
139 /* ===================================================================== */
140 /* .. Parameters .. */
141 /* .. */
142 /* .. Local Scalars .. */
143 /* .. */
144 /* .. External Functions .. */
145 /* .. */
146 /* .. External Subroutines .. */
147 /* .. */
148 /* .. Intrinsic Functions .. */
149 /* .. */
150 /* .. Executable Statements .. */
151 /* Parameter adjustments */
152 --ap;
153 --ipiv;
154 b_dim1 = *ldb;
155 b_offset = 1 + b_dim1;
156 b -= b_offset;
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 (*nrhs < 0)
169 {
170 *info = -3;
171 }
172 else if (*ldb < max(1,*n))
173 {
174 *info = -7;
175 }
176 if (*info != 0)
177 {
178 i__1 = -(*info);
179 xerbla_("DSPTRS", &i__1);
180 return 0;
181 }
182 /* Quick return if possible */
183 if (*n == 0 || *nrhs == 0)
184 {
185 return 0;
186 }
187 if (upper)
188 {
189 /* Solve A*X = B, where A = U*D*U**T. */
190 /* First solve U*D*X = B, overwriting B with X. */
191 /* K is the main loop index, decreasing from N to 1 in steps of */
192 /* 1 or 2, depending on the size of the diagonal blocks. */
193 k = *n;
194 kc = *n * (*n + 1) / 2 + 1;
195 L10: /* If K < 1, exit from loop. */
196 if (k < 1)
197 {
198 goto L30;
199 }
200 kc -= k;
201 if (ipiv[k] > 0)
202 {
203 /* 1 x 1 diagonal block */
204 /* Interchange rows K and IPIV(K). */
205 kp = ipiv[k];
206 if (kp != k)
207 {
208 dswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], ldb);
209 }
210 /* Multiply by inv(U(K)), where U(K) is the transformation */
211 /* stored in column K of A. */
212 i__1 = k - 1;
213 dger_(&i__1, nrhs, &c_b7, &ap[kc], &c__1, &b[k + b_dim1], ldb, &b[ b_dim1 + 1], ldb);
214 /* Multiply by the inverse of the diagonal block. */
215 d__1 = 1. / ap[kc + k - 1];
216 dscal_(nrhs, &d__1, &b[k + b_dim1], ldb);
217 --k;
218 }
219 else
220 {
221 /* 2 x 2 diagonal block */
222 /* Interchange rows K-1 and -IPIV(K). */
223 kp = -ipiv[k];
224 if (kp != k - 1)
225 {
226 dswap_(nrhs, &b[k - 1 + b_dim1], ldb, &b[kp + b_dim1], ldb);
227 }
228 /* Multiply by inv(U(K)), where U(K) is the transformation */
229 /* stored in columns K-1 and K of A. */
230 i__1 = k - 2;
231 dger_(&i__1, nrhs, &c_b7, &ap[kc], &c__1, &b[k + b_dim1], ldb, &b[ b_dim1 + 1], ldb);
232 i__1 = k - 2;
233 dger_(&i__1, nrhs, &c_b7, &ap[kc - (k - 1)], &c__1, &b[k - 1 + b_dim1], ldb, &b[b_dim1 + 1], ldb);
234 /* Multiply by the inverse of the diagonal block. */
235 akm1k = ap[kc + k - 2];
236 akm1 = ap[kc - 1] / akm1k;
237 ak = ap[kc + k - 1] / akm1k;
238 denom = akm1 * ak - 1.;
239 i__1 = *nrhs;
240 for (j = 1;
241 j <= i__1;
242 ++j)
243 {
244 bkm1 = b[k - 1 + j * b_dim1] / akm1k;
245 bk = b[k + j * b_dim1] / akm1k;
246 b[k - 1 + j * b_dim1] = (ak * bkm1 - bk) / denom;
247 b[k + j * b_dim1] = (akm1 * bk - bkm1) / denom;
248 /* L20: */
249 }
250 kc = kc - k + 1;
251 k += -2;
252 }
253 goto L10;
254 L30: /* Next solve U**T*X = B, overwriting B with X. */
255 /* K is the main loop index, increasing from 1 to N in steps of */
256 /* 1 or 2, depending on the size of the diagonal blocks. */
257 k = 1;
258 kc = 1;
259 L40: /* If K > N, exit from loop. */
260 if (k > *n)
261 {
262 goto L50;
263 }
264 if (ipiv[k] > 0)
265 {
266 /* 1 x 1 diagonal block */
267 /* Multiply by inv(U**T(K)), where U(K) is the transformation */
268 /* stored in column K of A. */
269 i__1 = k - 1;
270 dgemv_("Transpose", &i__1, nrhs, &c_b7, &b[b_offset], ldb, &ap[kc] , &c__1, &c_b19, &b[k + b_dim1], ldb);
271 /* Interchange rows K and IPIV(K). */
272 kp = ipiv[k];
273 if (kp != k)
274 {
275 dswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], ldb);
276 }
277 kc += k;
278 ++k;
279 }
280 else
281 {
282 /* 2 x 2 diagonal block */
283 /* Multiply by inv(U**T(K+1)), where U(K+1) is the transformation */
284 /* stored in columns K and K+1 of A. */
285 i__1 = k - 1;
286 dgemv_("Transpose", &i__1, nrhs, &c_b7, &b[b_offset], ldb, &ap[kc] , &c__1, &c_b19, &b[k + b_dim1], ldb);
287 i__1 = k - 1;
288 dgemv_("Transpose", &i__1, nrhs, &c_b7, &b[b_offset], ldb, &ap[kc + k], &c__1, &c_b19, &b[k + 1 + b_dim1], ldb);
289 /* Interchange rows K and -IPIV(K). */
290 kp = -ipiv[k];
291 if (kp != k)
292 {
293 dswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], ldb);
294 }
295 kc = kc + (k << 1) + 1;
296 k += 2;
297 }
298 goto L40;
299 L50:
300 ;
301 }
302 else
303 {
304 /* Solve A*X = B, where A = L*D*L**T. */
305 /* First solve L*D*X = B, overwriting B with X. */
306 /* K is the main loop index, increasing from 1 to N in steps of */
307 /* 1 or 2, depending on the size of the diagonal blocks. */
308 k = 1;
309 kc = 1;
310 L60: /* If K > N, exit from loop. */
311 if (k > *n)
312 {
313 goto L80;
314 }
315 if (ipiv[k] > 0)
316 {
317 /* 1 x 1 diagonal block */
318 /* Interchange rows K and IPIV(K). */
319 kp = ipiv[k];
320 if (kp != k)
321 {
322 dswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], ldb);
323 }
324 /* Multiply by inv(L(K)), where L(K) is the transformation */
325 /* stored in column K of A. */
326 if (k < *n)
327 {
328 i__1 = *n - k;
329 dger_(&i__1, nrhs, &c_b7, &ap[kc + 1], &c__1, &b[k + b_dim1], ldb, &b[k + 1 + b_dim1], ldb);
330 }
331 /* Multiply by the inverse of the diagonal block. */
332 d__1 = 1. / ap[kc];
333 dscal_(nrhs, &d__1, &b[k + b_dim1], ldb);
334 kc = kc + *n - k + 1;
335 ++k;
336 }
337 else
338 {
339 /* 2 x 2 diagonal block */
340 /* Interchange rows K+1 and -IPIV(K). */
341 kp = -ipiv[k];
342 if (kp != k + 1)
343 {
344 dswap_(nrhs, &b[k + 1 + b_dim1], ldb, &b[kp + b_dim1], ldb);
345 }
346 /* Multiply by inv(L(K)), where L(K) is the transformation */
347 /* stored in columns K and K+1 of A. */
348 if (k < *n - 1)
349 {
350 i__1 = *n - k - 1;
351 dger_(&i__1, nrhs, &c_b7, &ap[kc + 2], &c__1, &b[k + b_dim1], ldb, &b[k + 2 + b_dim1], ldb);
352 i__1 = *n - k - 1;
353 dger_(&i__1, nrhs, &c_b7, &ap[kc + *n - k + 2], &c__1, &b[k + 1 + b_dim1], ldb, &b[k + 2 + b_dim1], ldb);
354 }
355 /* Multiply by the inverse of the diagonal block. */
356 akm1k = ap[kc + 1];
357 akm1 = ap[kc] / akm1k;
358 ak = ap[kc + *n - k + 1] / akm1k;
359 denom = akm1 * ak - 1.;
360 i__1 = *nrhs;
361 for (j = 1;
362 j <= i__1;
363 ++j)
364 {
365 bkm1 = b[k + j * b_dim1] / akm1k;
366 bk = b[k + 1 + j * b_dim1] / akm1k;
367 b[k + j * b_dim1] = (ak * bkm1 - bk) / denom;
368 b[k + 1 + j * b_dim1] = (akm1 * bk - bkm1) / denom;
369 /* L70: */
370 }
371 kc = kc + (*n - k << 1) + 1;
372 k += 2;
373 }
374 goto L60;
375 L80: /* Next solve L**T*X = B, overwriting B with X. */
376 /* K is the main loop index, decreasing from N to 1 in steps of */
377 /* 1 or 2, depending on the size of the diagonal blocks. */
378 k = *n;
379 kc = *n * (*n + 1) / 2 + 1;
380 L90: /* If K < 1, exit from loop. */
381 if (k < 1)
382 {
383 goto L100;
384 }
385 kc -= *n - k + 1;
386 if (ipiv[k] > 0)
387 {
388 /* 1 x 1 diagonal block */
389 /* Multiply by inv(L**T(K)), where L(K) is the transformation */
390 /* stored in column K of A. */
391 if (k < *n)
392 {
393 i__1 = *n - k;
394 dgemv_("Transpose", &i__1, nrhs, &c_b7, &b[k + 1 + b_dim1], ldb, &ap[kc + 1], &c__1, &c_b19, &b[k + b_dim1], ldb);
395 }
396 /* Interchange rows K and IPIV(K). */
397 kp = ipiv[k];
398 if (kp != k)
399 {
400 dswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], ldb);
401 }
402 --k;
403 }
404 else
405 {
406 /* 2 x 2 diagonal block */
407 /* Multiply by inv(L**T(K-1)), where L(K-1) is the transformation */
408 /* stored in columns K-1 and K of A. */
409 if (k < *n)
410 {
411 i__1 = *n - k;
412 dgemv_("Transpose", &i__1, nrhs, &c_b7, &b[k + 1 + b_dim1], ldb, &ap[kc + 1], &c__1, &c_b19, &b[k + b_dim1], ldb);
413 i__1 = *n - k;
414 dgemv_("Transpose", &i__1, nrhs, &c_b7, &b[k + 1 + b_dim1], ldb, &ap[kc - (*n - k)], &c__1, &c_b19, &b[k - 1 + b_dim1], ldb);
415 }
416 /* Interchange rows K and -IPIV(K). */
417 kp = -ipiv[k];
418 if (kp != k)
419 {
420 dswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], ldb);
421 }
422 kc -= *n - k + 2;
423 k += -2;
424 }
425 goto L90;
426 L100:
427 ;
428 }
429 return 0;
430 /* End of DSPTRS */
431 }
432 /* dsptrs_ */
433