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