1 /*
2 * Copyright (C) 1998, 2000-2007, 2010, 2011, 2012, 2013 SINTEF ICT,
3 * Applied Mathematics, Norway.
4 *
5 * Contact information: E-mail: tor.dokken@sintef.no
6 * SINTEF ICT, Department of Applied Mathematics,
7 * P.O. Box 124 Blindern,
8 * 0314 Oslo, Norway.
9 *
10 * This file is part of SISL.
11 *
12 * SISL is free software: you can redistribute it and/or modify
13 * it under the terms of the GNU Affero General Public License as
14 * published by the Free Software Foundation, either version 3 of the
15 * License, or (at your option) any later version.
16 *
17 * SISL is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU Affero General Public License for more details.
21 *
22 * You should have received a copy of the GNU Affero General Public
23 * License along with SISL. If not, see
24 * <http://www.gnu.org/licenses/>.
25 *
26 * In accordance with Section 7(b) of the GNU Affero General Public
27 * License, a covered work must retain the producer line in every data
28 * file that is created or manipulated using SISL.
29 *
30 * Other Usage
31 * You can be released from the requirements of the license by purchasing
32 * a commercial license. Buying such a license is mandatory as soon as you
33 * develop commercial activities involving the SISL library without
34 * disclosing the source code of your own applications.
35 *
36 * This file may be used in accordance with the terms contained in a
37 * written agreement between you and SINTEF ICT.
38 */
39
40 #include "sisl-copyright.h"
41
42 /*
43 *
44 * $Id: s1927.c,v 1.3 2005-02-28 09:04:49 afr Exp $
45 *
46 */
47
48
49 #define S1927
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s1927(double * w1,int nur,int ik,int * ed,double * w2,int nrc,double * w3,int nlr,double * ex[],double * ey,int * jstat)55 s1927 (double *w1, int nur, int ik, int *ed, double *w2, int nrc,
56 double *w3, int nlr, double *ex[], double *ey, int *jstat)
57 #else
58 void
59 s1927 (w1, nur, ik, ed, w2, nrc, w3, nlr, ex, ey, jstat)
60 double *w1;
61 int nur;
62 int ik;
63 int *ed;
64 double *w2;
65 int nrc;
66 double *w3;
67 int nlr;
68 double *ex[];
69 double *ey;
70 int *jstat;
71 #endif
72 /*
73 *********************************************************************
74 *
75 *********************************************************************
76 *
77 * PURPOSE : To solve the equation W*ex = ey where W is an almost
78 * banded matrix whose LU-factorization is contained in
79 * w1, w2, w3.
80 *
81 *
82 * INPUT : w1 - Upper part of LU-factorization of W.
83 * Dimension (1:nur*ik-1).
84 * nur - Number of upper rows of W.
85 * ik - Number of non-zero entries in each upper
86 * row of W.
87 * ed - Pointers to the diagonal elements in w1.
88 * w2 - Right part of LU-factorization of W.
89 * Dimension (1:nur*nrc-1).
90 * nrc - Number of right columns of W.
91 * w3 - Lower part of LU-factorization of W.
92 * Dimension (1:nlr*(nur+nlr)-1).
93 * nlr - Number of lower rows.
94 * ey - Right side of equation.
95 * Dimension (1:nur+nlr-1).
96 * dim - Number of rows of the arrays ey and ex.
97 *
98 * OUTPUT : ex - Solution to W*ex=ey. Dimension (1:nur+nlr-1).
99 * jstat - Output status:
100 * < 0: Error.
101 * = 0: Ok.
102 * > 0: Warning.
103 *
104 * METHOD : A description of the form of the matrix W may be found in the
105 * subroutine s1898.
106 * Since the LU-factorization of W is contained in w1 and w3, the
107 * solution of W*ex=ey proceeds as follows:
108 * First solve L*z=ey, then solve U*ex=z.
109 *
110 * REFERENCES : Fortran version:
111 * E.Aarn[s, CP, 79-01
112 *
113 * CALLS :
114 *
115 * WRITTEN BY : Christophe R. Birkeland, SI, 1991-06
116 * REWRITTEN BY :
117 * REVISED BY :
118 *
119 *********************************************************************
120 */
121 {
122 int kpos = 0;
123 int ii, jj; /* Loop control parameters */
124 int di; /* Pointer to diagonal element of W */
125 int midi; /* Parameter always equal: ii-di */
126 int dim; /* di minus 2: di-2 */
127 int mur; /* Used in calculation of index for w3 */
128 int nn; /* Number of rows/columns in w3 */
129 int nlc; /* Number of left columns in W */
130 double wii; /* Used to store values from matrix W */
131 double sum; /* Stores values for calculation of ex */
132
133 *jstat = 0;
134
135
136 /* Test if legal dimension of interpolatoin problem */
137
138 if (nur < 1 || ik < 1 || nrc < 0 || nlr < 0)
139 goto err160;
140 nn = nur + nlr;
141 nlc = nn - nrc;
142 if (ik > nlc)
143 goto err160;
144
145
146 /* Allocate output array ex */
147
148 *ex = new0array (nn, DOUBLE);
149 if (*ex == SISL_NULL)
150 goto err101;
151
152
153 /* Solve L*z = ey */
154
155 for (ii = 0; ii < nur; ii++)
156 {
157 di = ed[ii];
158 wii = w1[(di - 1) * nur + ii];
159
160
161 /* Test for errors */
162
163 if (ii >= nlc)
164 goto err163;
165 if (di < 1 || ik < di || wii == (double) 0.0)
166 goto err162;
167 sum = ey[ii];
168 if (di > 1)
169 {
170 dim = di - 1;
171 midi = ii - di + 1;
172 for (jj = 0; jj < dim; jj++)
173 sum -= w1[jj * nur + ii] * ((*ex)[jj + midi]);
174 }
175 (*ex)[ii] = sum / wii;
176 }
177
178 /* Solve filled part of L*z = ey */
179
180 for (; ii < nn; ii++)
181 {
182 mur = ii - nur;
183 wii = w3[ii * nlr + mur];
184 if (wii == (double) 0.0)
185 goto err162;
186 sum = ey[ii];
187 if (ii >= 1)
188 {
189 for (jj = 0; jj < ii; jj++)
190 sum -= w3[jj * nlr + mur] * ((*ex)[jj]);
191 }
192 (*ex)[ii] = sum / wii;
193 }
194
195 /* Solve U*ex = z ; Jump if filled part of U is exhausted */
196
197 for (ii = nn - 2; ii >= nur; ii--)
198 {
199 sum = (*ex)[ii];
200 mur = ii - nur;
201 for (jj = ii + 1; jj < nn; jj++)
202 sum -= w3[jj * nlr + mur] * ((*ex)[jj]);
203 (*ex)[ii] = sum;
204 }
205
206 /* Test if w2 contains diagonal elements */
207
208 if (ii >= nlc)
209 goto err163;
210 if (nlc < nn)
211 {
212 for (; ii >= 0; ii--)
213 {
214 sum = (*ex)[ii];
215 for (jj = nlc; jj < nn; jj++)
216 sum -= w2[(jj - nlc) * nur + ii] * ((*ex)[jj]);
217 (*ex)[ii] = sum;
218 }
219 }
220 for (ii = nur - 1; ii >= 0; ii--)
221 {
222 di = ed[ii];
223 if (di < ik)
224 {
225 sum = (*ex)[ii];
226 midi = ii - di + 1;
227 for (jj = di; jj < ik; jj++)
228 sum -= w1[jj * nur + ii] * ((*ex)[jj + midi]);
229 (*ex)[ii] = sum;
230 }
231 }
232
233 goto out;
234
235
236 /* Memory error, array ex not allocated */
237
238 err101:
239 *jstat = -101;
240 s6err ("s1927", *jstat, kpos);
241 goto out;
242
243 /* error in dimension of interpolation problem */
244
245 err160:
246 *jstat = -160;
247 s6err ("s1927", *jstat, kpos);
248 goto out;
249
250 /* W is non-invertible */
251
252 err162:
253 *jstat = -162;
254 s6err ("s1927", *jstat, kpos);
255 goto out;
256
257 /* w2 contains diagonal element */
258
259 err163:
260 *jstat = -163;
261 s6err ("s1927", *jstat, kpos);
262 goto out;
263
264 out:
265 return;
266 }
267