66 fprintf(stderr,
"Error in TimeGLMGEEInitialize(): %s is not a supported ",type);
67 fprintf(stderr,
"multi-stage time integration scheme of class %s.\n",
class);
74 params->
A_yyt = (
double*) calloc (s*s,
sizeof(
double));
75 params->
B_yyt = (
double*) calloc (s*r,
sizeof(
double));
76 params->
C_yyt = (
double*) calloc (s*r,
sizeof(
double));
77 params->
D_yyt = (
double*) calloc (r*r,
sizeof(
double));
78 params->
c_yyt = (
double*) calloc (s ,
sizeof(
double));
85 params->
A_yeps = (
double*) calloc (s*s,
sizeof(
double));
86 params->
B_yeps = (
double*) calloc (s*r,
sizeof(
double));
87 params->
C_yeps = (
double*) calloc (s*r,
sizeof(
double));
88 params->
D_yeps = (
double*) calloc (r*r,
sizeof(
double));
89 params->
c_yeps = (
double*) calloc (s ,
sizeof(
double));
98 params->
A_yeps[1*s+0] = 1.0;
101 params->
B_yeps[0*s+0] = params->
B_yeps[0*s+1] = 1.0/12.0; params->
B_yeps[0*s+2] = 5.0/6.0;
102 params->
B_yeps[1*s+0] = params->
B_yeps[1*s+1] = 1.0/12.0; params->
B_yeps[1*s+2] = -1.0/6.0;
104 params->
C_yeps[0*r+0] = 1.0;
105 params->
C_yeps[1*r+0] = 1.0; params->
C_yeps[1*r+1] = 10.0;
106 params->
C_yeps[2*r+0] = 1.0; params->
C_yeps[2*r+1] = -1.0;
108 params->
D_yeps[0*r+0] = 1.0;
109 params->
D_yeps[1*r+1] = 1.0;
111 params->
A_yyt[1*s+0] = 1.0;
112 params->
A_yyt[2*s+0] = params->
A_yyt[2*s+1] = 0.25;
114 params->
B_yyt[0*s+0] = params->
B_yyt[0*s+1] = 1.0/12.0; params->
B_yyt[0*s+2] = 5.0/6.0;
115 params->
B_yyt[1*s+0] = params->
B_yyt[1*s+1] = 1.0/6.0; params->
B_yyt[1*s+2] = 2.0/3.0;
117 params->
C_yyt[0*r+0] = 1.0;
118 params->
C_yyt[1*r+0] = -9.0; params->
C_yyt[1*r+1] = 10.0;
119 params->
C_yyt[2*r+0] = 2.0; params->
C_yyt[2*r+1] = -1.0;
121 params->
D_yyt[0*r+0] = 1.0;
122 params->
D_yyt[1*r+1] = 1.0;
126 params->
A_yyt[1*s+0] = 0.75;
127 params->
A_yyt[2*s+0] = 0.25;
128 params->
A_yyt[2*s+1] = 29.0/60.0;
129 params->
A_yyt[3*s+0] = -21.0/44.0;
130 params->
A_yyt[3*s+1] = 145.0/44.0;
131 params->
A_yyt[3*s+2] = -20.0/11.0;
133 params->
B_yyt[0*s+0] = 109.0/275.0;
134 params->
B_yyt[0*s+1] = 58.0/75.0;
135 params->
B_yyt[0*s+2] = -37.0/110.0;
136 params->
B_yyt[0*s+3] = 1.0/6.0;
137 params->
B_yyt[1*s+0] = 3.0/11.0;
138 params->
B_yyt[1*s+2] = 75.0/88.0;
139 params->
B_yyt[1*s+3] = -1.0/8.0;
141 params->
C_yyt[0*r+1] = 1.0;
142 params->
C_yyt[1*r+0] = 75.0/58.0;
143 params->
C_yyt[1*r+1] = -17.0/58.0;
144 params->
C_yyt[2*r+1] = params->
C_yyt[3*r+1] = 1.0;
146 params->
D_yyt[0*r+0] = 1.0;
147 params->
D_yyt[1*r+1] = 1.0;
149 double T[r*r],Tinv[r*r];
153 T[1*r+1] = 1.0-params->
gamma;
163 params->
A_yyt[1*s+0]=-0.94079244066783383269;
164 params->
A_yyt[2*s+0]= 0.64228187778301907108;
165 params->
A_yyt[2*s+1]= 0.10915356933958500042;
166 params->
A_yyt[3*s+0]=-0.51764297742287450812;
167 params->
A_yyt[3*s+1]= 0.74414270351096040738;
168 params->
A_yyt[3*s+2]=-0.71404164927824538121;
169 params->
A_yyt[4*s+0]=-0.44696561556825969206;
170 params->
A_yyt[4*s+1]=-0.76768425657590196518;
171 params->
A_yyt[4*s+2]= 0.20111608138142987881;
172 params->
A_yyt[4*s+3]= 0.93828186737840469796;
174 params->
B_yyt[0*s+0]=-0.029309178948150356153;
175 params->
B_yyt[0*s+1]=-0.49671981884013874923;
176 params->
B_yyt[0*s+2]= 0.34275801517650053274;
177 params->
B_yyt[0*s+3]= 0.32941112623949194988;
178 params->
B_yyt[0*s+4]= 0.85385985637229662276;
179 params->
B_yyt[1*s+0]= 0.78133219686062535272;
180 params->
B_yyt[1*s+1]= 0.074238691892675897635;
181 params->
B_yyt[1*s+2]= 0.57957363498384957966;
182 params->
B_yyt[1*s+3]=-0.24638502829674959968;
183 params->
B_yyt[1*s+4]=-0.18875949544040123033;
185 params->
C_yyt[0*r+0]= 0.16911424754448327735;
186 params->
C_yyt[0*r+1]= 0.83088575245551672265;
187 params->
C_yyt[1*r+0]= 0.53638465733199574340;
188 params->
C_yyt[1*r+1]= 0.46361534266800425660;
189 params->
C_yyt[2*r+0]= 0.39901579167169582526;
190 params->
C_yyt[2*r+1]= 0.60098420832830417474;
191 params->
C_yyt[3*r+0]= 0.87689005530618575480;
192 params->
C_yyt[3*r+1]= 0.12310994469381424520;
193 params->
C_yyt[4*r+0]= 0.99056100455550913009;
194 params->
C_yyt[4*r+1]= 0.0094389954444908699092;
196 params->
D_yyt[0*r+0] = 1.0;
197 params->
D_yyt[1*r+1] = 1.0;
199 double T[r*r],Tinv[r*r];
203 T[1*r+1] = 1.0-params->
gamma;
213 params->
A_yyt[1*s+0] = - 2169604947363702313.0 / 24313474998937147335.0;
214 params->
A_yyt[2*s+0] = 46526746497697123895.0 / 94116917485856474137.0;
215 params->
A_yyt[2*s+1] = -10297879244026594958.0 / 49199457603717988219.0;
216 params->
A_yyt[3*s+0] = 23364788935845982499.0 / 87425311444725389446.0;
217 params->
A_yyt[3*s+1] = -79205144337496116638.0 / 148994349441340815519.0;
218 params->
A_yyt[3*s+2] = 40051189859317443782.0 / 36487615018004984309.0;
219 params->
A_yyt[4*s+0] = 42089522664062539205.0 / 124911313006412840286.0;
220 params->
A_yyt[4*s+1] = -15074384760342762939.0 / 137927286865289746282.0;
221 params->
A_yyt[4*s+2] = -62274678522253371016.0 / 125918573676298591413.0;
222 params->
A_yyt[4*s+3] = 13755475729852471739.0 / 79257927066651693390.0;
224 params->
B_yyt[0*s+0] = 61546696837458703723.0 / 56982519523786160813.0;
225 params->
B_yyt[0*s+1] = -55810892792806293355.0 / 206957624151308356511.0;
226 params->
B_yyt[0*s+2] = 24061048952676379087.0 / 158739347956038723465.0;
227 params->
B_yyt[0*s+3] = 3577972206874351339.0 / 7599733370677197135.0;
228 params->
B_yyt[0*s+4] = -59449832954780563947.0 / 137360038685338563670.0;
229 params->
B_yyt[1*s+0] = - 9738262186984159168.0 / 99299082461487742983.0;
230 params->
B_yyt[1*s+1] = -32797097931948613195.0 / 61521565616362163366.0;
231 params->
B_yyt[1*s+2] = 42895514606418420631.0 / 71714201188501437336.0;
232 params->
B_yyt[1*s+3] = 22608567633166065068.0 / 55371917805607957003.0;
233 params->
B_yyt[1*s+4] = 94655809487476459565.0 / 151517167160302729021.0;
235 params->
C_yyt[0*r+0] = 70820309139834661559.0 / 80863923579509469826.0;
236 params->
C_yyt[0*r+1] = 10043614439674808267.0 / 80863923579509469826.0;
237 params->
C_yyt[1*r+0] = 161694774978034105510.0 / 106187653640211060371.0;
238 params->
C_yyt[1*r+1] = -55507121337823045139.0 / 106187653640211060371.0;
239 params->
C_yyt[2*r+0] = 78486094644566264568.0 / 88171030896733822981.0;
240 params->
C_yyt[2*r+1] = 9684936252167558413.0 / 88171030896733822981.0;
241 params->
C_yyt[3*r+0] = 65394922146334854435.0 / 84570853840405479554.0;
242 params->
C_yyt[3*r+1] = 19175931694070625119.0 / 84570853840405479554.0;
243 params->
C_yyt[4*r+0] = 8607282770183754108.0 / 108658046436496925911.0;
244 params->
C_yyt[4*r+1] = 100050763666313171803.0 / 108658046436496925911.0;
246 params->
D_yyt[0*r+0] = 1.0;
247 params->
D_yyt[1*r+1] = 1.0;
249 double T[r*r],Tinv[r*r];
253 T[1*r+1] = 1.0-params->
gamma;
263 params->
A_yeps[1*s+0] = 1.0;
264 params->
A_yeps[3*s+2] = 0.5;
265 params->
A_yeps[4*s+2] = 0.25;
266 params->
A_yeps[4*s+3] = 0.25;
267 params->
A_yeps[5*s+2] = 0.25;
268 params->
A_yeps[5*s+3] = 0.25;
269 params->
A_yeps[5*s+4] = 0.5;
272 params->
B_yeps[1*s+0] = params->
B_yeps[1*s+1] = -2.0/3.0;
274 = params->
B_yeps[1*s+5] = 1.0/3.0;
276 params->
C_yeps[0*r+0] = 1.0;
277 params->
C_yeps[1*r+0] = 1.0;
278 params->
C_yeps[2*r+0] = 1.0; params->
C_yeps[2*r+1] = 0.75;
279 params->
C_yeps[3*r+0] = 1.0; params->
C_yeps[3*r+1] = 0.75;
280 params->
C_yeps[4*r+0] = 1.0; params->
C_yeps[4*r+1] = 0.75;
281 params->
C_yeps[5*r+0] = 1.0; params->
C_yeps[5*r+1] = 0.75;
283 params->
D_yeps[0*r+0] = 1.0;
284 params->
D_yeps[1*r+1] = 1.0;
286 params->
A_yyt[1*s+0] = 1.0;
287 params->
A_yyt[3*s+2] = 0.5;
288 params->
A_yyt[4*s+2] = 0.25;
289 params->
A_yyt[4*s+3] = 0.25;
290 params->
A_yyt[5*s+2] = 0.25;
291 params->
A_yyt[5*s+3] = 0.25;
292 params->
A_yyt[5*s+4] = 0.5;
294 params->
B_yyt[0*s+0] = params->
B_yyt[0*s+1] = 0.5;
296 = params->
B_yyt[1*s+5] = 0.25;
298 params->
C_yyt[0*r+0] = 1.0;
299 params->
C_yyt[1*r+0] = 1.0;
300 params->
C_yyt[2*r+1] = 1.0;
301 params->
C_yyt[3*r+1] = 1.0;
302 params->
C_yyt[4*r+1] = 1.0;
303 params->
C_yyt[5*r+1] = 1.0;
305 params->
D_yyt[0*r+0] = 1.0;
306 params->
D_yyt[1*r+1] = 1.0;
310 params->
A_yeps[1*s+0] = 0.5;
312 params->
A_yeps[2*s+0] = -1.0;
313 params->
A_yeps[2*s+1] = 2.0;
315 params->
A_yeps[3*s+0] = 1.0/6.0;
316 params->
A_yeps[3*s+1] = 2.0/3.0;
317 params->
A_yeps[3*s+2] = 1.0/6.0;
319 params->
A_yeps[5*s+0] = -7.0/24.0;
320 params->
A_yeps[5*s+1] = 1.0/3.0;
321 params->
A_yeps[5*s+2] = 1.0/12.0;
322 params->
A_yeps[5*s+3] = -1.0/8.0;
323 params->
A_yeps[5*s+4] = 0.5;
325 params->
A_yeps[6*s+0] = 7.0/6.0;
326 params->
A_yeps[6*s+1] = -4.0/3.0;
327 params->
A_yeps[6*s+2] = -1.0/3.0;
328 params->
A_yeps[6*s+3] = 0.5;
329 params->
A_yeps[6*s+4] = -1.0;
330 params->
A_yeps[6*s+5] = 2.0;
332 params->
A_yeps[7*s+4] = 1.0/6.0;
333 params->
A_yeps[7*s+5] = 2.0/3.0;
334 params->
A_yeps[7*s+6] = 1.0/6.0;
336 params->
B_yeps[0*s+0] = 1.0/6.0;
337 params->
B_yeps[0*s+1] = 2.0/3.0;
338 params->
B_yeps[0*s+2] = 1.0/6.0;
340 params->
B_yeps[1*s+0] = -1.0/6.0;
341 params->
B_yeps[1*s+1] = -2.0/3.0;
342 params->
B_yeps[1*s+2] = -1.0/6.0;
343 params->
B_yeps[1*s+4] = 1.0/6.0;
344 params->
B_yeps[1*s+5] = 2.0/3.0;
345 params->
B_yeps[1*s+6] = 1.0/6.0;
347 params->
C_yeps[0*r+0] = 1.0;
348 params->
C_yeps[1*r+0] = 1.0;
349 params->
C_yeps[2*r+0] = 1.0;
350 params->
C_yeps[3*r+0] = 1.0;
351 params->
C_yeps[4*r+0] = 1.0;
352 params->
C_yeps[5*r+0] = 1.0;
353 params->
C_yeps[6*r+0] = 1.0;
354 params->
C_yeps[7*r+0] = 1.0;
356 params->
C_yeps[4*r+1] = 1.0;
357 params->
C_yeps[5*r+1] = 1.0;
358 params->
C_yeps[6*r+1] = 1.0;
359 params->
C_yeps[7*r+1] = 1.0;
361 params->
D_yeps[0*r+0] = 1.0;
362 params->
D_yeps[1*r+1] = 1.0;
364 double T[r*r],Tinv[r*r];
368 T[1*r+1] = 1.0-params->
gamma;
378 params->
A_yeps[1*s+0] = 0.585786437626904966;
379 params->
A_yeps[2*s+0] = 0.149999999999999994;
380 params->
A_yeps[2*s+1] = 0.849999999999999978;
381 params->
A_yeps[4*s+3] = 0.292893218813452483;
382 params->
A_yeps[5*s+3] = 0.0749999999999999972;
383 params->
A_yeps[5*s+4] = 0.424999999999999989;
384 params->
A_yeps[6*s+3] = 0.176776695296636893;
385 params->
A_yeps[6*s+4] = 0.176776695296636893;
386 params->
A_yeps[6*s+5] = 0.146446609406726241;
387 params->
A_yeps[7*s+3] = 0.176776695296636893;
388 params->
A_yeps[7*s+4] = 0.176776695296636893;
389 params->
A_yeps[7*s+5] = 0.146446609406726241;
390 params->
A_yeps[7*s+6] = 0.292893218813452483;
391 params->
A_yeps[8*s+3] = 0.176776695296636893;
392 params->
A_yeps[8*s+4] = 0.176776695296636893;
393 params->
A_yeps[8*s+5] = 0.146446609406726241;
394 params->
A_yeps[8*s+6] = 0.0749999999999999972;
395 params->
A_yeps[8*s+7] = 0.424999999999999989;
397 params->
B_yeps[0*s+0] = 0.353553390593273786;
398 params->
B_yeps[0*s+1] = 0.353553390593273786;
399 params->
B_yeps[0*s+2] = 0.292893218813452483;
401 params->
B_yeps[1*s+0] =-0.471404520791031678;
402 params->
B_yeps[1*s+1] =-0.471404520791031678;
403 params->
B_yeps[1*s+2] =-0.390524291751269959;
404 params->
B_yeps[1*s+3] = 0.235702260395515839;
405 params->
B_yeps[1*s+4] = 0.235702260395515839;
406 params->
B_yeps[1*s+5] = 0.195262145875634979;
407 params->
B_yeps[1*s+6] = 0.235702260395515839;
408 params->
B_yeps[1*s+7] = 0.235702260395515839;
409 params->
B_yeps[1*s+8] = 0.195262145875634979;
411 params->
C_yeps[0*r+0] = 1.0;
412 params->
C_yeps[1*r+0] = 1.0;
413 params->
C_yeps[2*r+0] = 1.0;
414 params->
C_yeps[3*r+0] = 1.0; params->
C_yeps[3*r+1] = 0.75;
415 params->
C_yeps[4*r+0] = 1.0; params->
C_yeps[4*r+1] = 0.75;
416 params->
C_yeps[5*r+0] = 1.0; params->
C_yeps[5*r+1] = 0.75;
417 params->
C_yeps[6*r+0] = 1.0; params->
C_yeps[6*r+1] = 0.75;
418 params->
C_yeps[7*r+0] = 1.0; params->
C_yeps[7*r+1] = 0.75;
419 params->
C_yeps[8*r+0] = 1.0; params->
C_yeps[8*r+1] = 0.75;
421 params->
D_yeps[0*r+0] = 1.0;
422 params->
D_yeps[1*r+1] = 1.0;
424 double T[r*r],Tinv[r*r];
428 T[1*r+1] = 1.0-params->
gamma;
438 for (i=0; i<s; i++) {
439 for (j=0; j<s; j++) {
448 in = fopen(
"glm_gee.inp",
"r");
451 printf(
"Reading GLM-GEE method parameters from glm_gee.inp.\n");
453 ferr = fscanf(in,
"%s",word);
if (ferr != 1)
return(1);
454 if (!strcmp(word,
"begin")) {
455 while (strcmp(word,
"end")) {
456 ferr = fscanf(in,
"%s",word);
if (ferr != 1)
return(1);
457 if (!strcmp(word,
"ee_mode")) {
458 ferr = fscanf(in,
"%s",params->
ee_mode);
459 if (ferr != 1)
return(1);
460 }
else if (strcmp(word,
"end")) {
462 ferr = fscanf(in,
"%s",useless);
if (ferr != 1)
return(ferr);
463 printf(
"Warning: keyword %s in file \"glm_gee.inp\" with value %s not ",word,useless);
464 printf(
"recognized or extraneous. Ignoring.\n");
468 fprintf(stderr,
"Error: Illegal format in file \"glm_gee.inp\".\n");
473 fprintf(stderr,
"Error in TimeGLMGEEInitialize(): %s is not a valid value for ",params->
ee_mode);
478 printf(
"GLM-GEE time integration error estimation mode: %s\n",params->
ee_mode);
483 params->
A = params->
A_yyt;
484 params->
B = params->
B_yyt;
485 params->
C = params->
C_yyt;
486 params->
D = params->
D_yyt;
487 params->
c = params->
c_yyt;
497 fprintf(stderr,
"Error in TimeGLMGEEInitialize(): Code should not have ");
498 fprintf(stderr,
"reached here for %s class of time-integrators. This is a ",
class);
499 fprintf(stderr,
"coding mistake.\n");
Contains macros and function definitions for common matrix operations.
MPI related function definitions.
#define _MatrixMultiplyNonSquare_(A, B, C, NRowA, NColA, NColB)
Some basic definitions and macros.
Contains function declarations for time integration.
int MPIBroadcast_character(char *, int, int, void *)
#define _MAX_STRING_SIZE_
char ee_mode[_MAX_STRING_SIZE_]
#define _ArraySetValue_(x, size, value)
int TimeGLMGEEInitialize(char *class, char *type, void *s, void *m)
Structure of MPI-related variables.
Structure containing the parameters for the General Linear Methods with Global Error Estimators (GLM-...
#define _MatrixInvert_(A, B, N)
#define _GLM_GEE_RK285EX_
#define _ArrayCopy1D_(x, y, size)
Contains macros and function definitions for common array operations.