33 #ifdef USE_RUNTIME_LOADING
40 #include "ac/getopt.h"
52 struct integration_data {
61 static int open_module(
const char* module);
63 static int method_multistep(
const char*, integration_data*,
void*,
const char*);
64 static int method_hope(
const char*, integration_data*,
void*,
const char*);
65 static int method_cubic(
const char*, integration_data*,
void*,
const char*);
66 static int method_cn(
const char*, integration_data*,
void*,
const char*);
68 void* get_method_data(
int,
const char*);
70 static struct funcs *ff = NULL;
71 static bool print_help =
false;
72 static void* p_data = NULL;
75 main(
int argn,
char *
const argv[])
89 {
"ms", METHOD_MULTISTEP },
90 {
"hope", METHOD_HOPE },
91 {
"cubic", METHOD_CUBIC },
92 {
"crank-nicolson", METHOD_CRANK_NICOLSON },
96 int curr_method = METHOD_UNKNOWN;
98 char* user_defined = 0;
99 void* method_data = 0;
100 integration_data d = { 0., 1., 1.e-3, 1.e-6, 10 };
103 const char optstring[] =
"i:m:t:T:n:r:u:h";
104 const struct option longopts[] = {
105 {
"integrator", required_argument, 0,
int(
'i') },
106 {
"method-data", required_argument, 0,
int(
'm') },
107 {
"timing", required_argument, 0,
int(
't') },
108 {
"tolerance", required_argument, 0,
int(
'T') },
109 {
"iterations", required_argument, 0,
int(
'n') },
110 {
"rho", required_argument, 0,
int(
'r') },
111 {
"user-defined", required_argument, 0,
int(
'u') },
112 {
"help", no_argument, 0,
int(
'h') },
113 {
"print-help", no_argument, 0,
int(
'H') },
115 { 0, no_argument, 0,
int(0) }
121 curropt = getopt_long(argn, argv, optstring, longopts, 0);
123 if (curropt == EOF) {
135 std::cerr <<
"missing parameter for option -"
140 std::cerr <<
"usage: int -[imtTnruh] [module]" << std::endl
142 <<
" -i,--integrator : integration method" << std::endl
143 <<
" -m,--method-data : method-dependent data" << std::endl
144 <<
" -t,--timing : ti:dt:tf" << std::endl
145 <<
" -T,--tolerance" << std::endl
146 <<
" -n,--maxiter" << std::endl
147 <<
" -r,--rho : asymptotic radius" << std::endl
148 <<
" -u,--user : user-defined data" << std::endl
150 <<
" -h,--help : print this message and exit" << std::endl;
160 if (strcmp(p->s,
optarg) == 0) {
167 std::cerr <<
"unknown integrator "
175 method_data = get_method_data(curr_method,
optarg);
179 char* s =
new char[strlen(
optarg)+1];
181 char* p = std::strrchr(s,
':');
183 std::cerr <<
"syntax: ti:dt:tf" << std::endl;
188 p = std::strrchr(s,
':');
190 std::cerr <<
"syntax: ti:dt:tf" << std::endl;
226 if (::ff->
size == 0) {
227 std::cerr <<
"no \"size\" func in module" << std::endl;
231 if (::ff->
func == 0) {
232 std::cerr <<
"no \"func\" func in module" << std::endl;
236 if (::ff->
grad == 0) {
237 std::cerr <<
"no \"grad\" func in module" << std::endl;
241 if (::ff->
read && ::ff->
read(&p_data, user_defined)) {
242 std::cerr <<
"unable to read(" << user_defined <<
") "
243 "data for module \"" << module <<
"\"" << std::endl;
249 ::ff->help(p_data, std::cerr);
251 std::cerr <<
"help not available for module "
252 "\"" << module <<
"\" (ignored)" << std::endl;
257 switch (curr_method) {
258 case METHOD_MULTISTEP :
259 rc = method_multistep(module, &d, method_data, user_defined);
262 rc = method_hope(module, &d, method_data, user_defined);
265 rc = method_cubic(module, &d, method_data, user_defined);
267 case METHOD_CRANK_NICOLSON:
268 rc = method_cn(module, &d, method_data, user_defined);
271 std::cerr <<
"you must select an integration method" << std::endl;
279 get_method_data(
int curr_method,
const char*
optarg)
281 switch (curr_method) {
283 std::cerr <<
"not implemented yet" << std::endl;
291 open_module(
const char* module)
298 std::cerr <<
"lt_dlinit() failed" << std::endl;
302 handle = lt_dlopen(module);
305 std::cerr <<
"lt_dlopen(\"" << module
306 <<
"\") returned \"" << err <<
"\"" << std::endl;
310 struct funcs **pf = (
struct funcs **)lt_dlsym(handle,
"ff");
313 std::cerr <<
"lt_dlsym(\"ff\") returned \"" << err <<
"\""
320 std::cerr <<
"invalid \"ff\" symbol in \"" << module <<
"\""
345 method_multistep(
const char* module, integration_data* d,
346 void* method_data,
const char* user_defined)
350 ::ff->read(&p_data, user_defined);
353 int size = ::ff->size(p_data);
391 int maxiter = d->maxiter;
395 doublereal beta = (4.*rho*rho-(1.-rho)*(1.-rho))/(4.-(1.-rho)*(1.-rho));
396 doublereal delta = (1.-rho)*(1.-rho)/(2.*(4.-(1.-rho)*(1.-rho)));
406 ::ff->init(p_data, *pX);
407 ::ff->func(p_data, *pXP, *pX, t);
408 for (
int k = 1; k <=
size; k++) {
416 std::cout << ti <<
" " << 0. <<
" ";
417 ::ff->out(p_data, std::cout, *pX, *pXP) << std::endl;
419 flip(&pX, &pXP, &pXm1, &pXPm1, &pXm2, &pXPm2);
424 for (
int k = size; k > 0; k--) {
429 doublereal x = -4.*xm1+5.*xm2+dt*(4.*xPm1+2.*xPm2);
431 R.PutCoef(k, a1*xm1+a2*xm2+dt*(b1*xPm1+b2*xPm2));
439 ::ff->func(p_data, *pXP, *pX, t);
440 for (
int k = 1; k <=
size; k++) {
451 std::cerr <<
"current iteration " << j
452 <<
" exceeds max iteration number "
453 << maxiter << std::endl;
460 ::ff->grad(p_data, J, *pX, t);
461 for (
int k = 1; k <=
size; k++) {
462 for (
int l = 1; l <=
size; l++) {
463 Jac.
PutCoef(k, l, -coef*J(k, l));
470 for (
int k = size; k > 0; k--) {
477 std::cout << t <<
" " << test <<
" ";
478 ::ff->out(p_data, std::cout, *pX, *pXP) << std::endl;
480 flip(&pX, &pXP, &pXm1, &pXPm1, &pXm2, &pXPm2);
483 ::ff->destroy(&p_data);
489 method_hope(
const char* module, integration_data* d,
490 void* method_data,
const char* user_defined)
492 std::cerr << __FUNCTION__ <<
"not implemented yet!" << std::endl;
499 method_cubic(
const char* module, integration_data* d,
500 void* method_data,
const char* user_defined)
504 ::ff->read(&p_data, user_defined);
507 int size = ::ff->size(p_data);
547 int maxiter = d->maxiter;
563 ::ff->init(p_data, *pX);
564 ::ff->func(p_data, *pXP, *pX, t);
565 for (
int k = 1; k <=
size; k++) {
573 std::cout << ti <<
" " << 0. <<
" ";
574 ::ff->out(p_data, std::cout, *pX, *pXP) << std::endl;
576 flip(&pX, &pXP, &pXm1, &pXPm1, &pXm2, &pXPm2);
581 for (
int k = 1; k <=
size; k++) {
586 doublereal x = -4.*xm1+5.*xm2+dt*(4.*xPm1+2.*xPm2);
595 ::ff->func(p_data, *pXP, *pX, t);
596 for (
int k = 1; k <=
size; k++) {
601 doublereal xz = m0*x+m1*xm1+dt*(n0*xP+n1*xPm1);
605 ::ff->func(p_data, XPz, Xz, t+z*dt);
606 for (
int k = 1; k <=
size; k++) {
608 w1*pXPm1->operator()(k)
610 + w0*pXP->operator()(k)
613 - pXm1->operator()(k)
623 std::cerr <<
"current iteration " << j
624 <<
" exceeds max iteration number "
625 << maxiter << std::endl;
633 ::ff->grad(p_data, Jz, Xz, t+z*dt);
634 ::ff->grad(p_data, J0, *pX, t);
635 for (
int k = 1; k <=
size; k++) {
636 for (
int l = 1; l <=
size; l++) {
638 for (
int m = 1; m <=
size; m++) {
639 d += Jz(k, m)*J0(m, l);
641 d = -dt*(wz*(Jz(k, l)+dt*n0*d)+w0*J0(k, l));
649 for (
int k = size; k > 0; k--) {
656 std::cout << t <<
" " << test <<
" ";
657 ::ff->out(p_data, std::cout, *pX, *pXP) << std::endl;
659 flip(&pX, &pXP, &pXm1, &pXPm1, &pXm2, &pXPm2);
662 ::ff->destroy(&p_data);
668 method_cn(
const char* module, integration_data* d,
669 void* method_data,
const char* user_defined)
671 std::cerr << __FUNCTION__ <<
"not implemented yet!" << std::endl;
677 #else // ! USE_RUNTIME_LOADING
685 std::cerr <<
"Need dynamic load capabilities" << std::endl;
689 #endif // ! USE_RUNTIME_LOADING
virtual void Reset(void)=0
virtual VectorHandler * pResHdl(void) const =0
virtual void IncCoef(integer iRow, integer iCol, const doublereal &dCoef)
SolutionManager *const GetSolutionManager(integer iNLD, integer iLWS=0) const
virtual void IncCoef(integer iRow, const doublereal &dCoef)=0
virtual doublereal Norm(void) const
virtual MatrixHandler * pMatHdl(void) const =0
virtual void MatrReset(void)=0
virtual void PutCoef(integer iRow, integer iCol, const doublereal &dCoef)
virtual void Solve(void)=0
virtual void PutCoef(integer iRow, const doublereal &dCoef)=0
virtual VectorHandler * pSolHdl(void) const =0
static const std::vector< doublereal > v0