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