Index: /issm/trunk-jpl/src/wrappers/M1qn3/M1qn3.cpp =================================================================== --- /issm/trunk-jpl/src/wrappers/M1qn3/M1qn3.cpp (revision 19011) +++ /issm/trunk-jpl/src/wrappers/M1qn3/M1qn3.cpp (revision 19011) @@ -0,0 +1,178 @@ +/*!\file M1qn3.c + * \brief: data interpolation from a list of (x,y,values) into mesh vertices +*/ + +#include "./M1qn3.h" + +/*m1qn3 prototypes {{{*/ +extern "C" void *ctonbe_; // DIS mode : Conversion +extern "C" void *ctcabe_; // DIS mode : Conversion +extern "C" void *euclid_; // Scalar product +typedef void (*SimulFunc) (long* indic,long* n, double* x, double* pf,double* g,long [],float [],void* dzs); +extern "C" void m1qn3_ (void f(long* indic,long* n, double* x, double* pf,double* g,long [],float [],void* dzs), + void **, void **, void **, + long *, double [], double *, double [], double*, double *, + double *, char [], long *, long *, long *, long *, long *, long *, long [], double [], long *, + long *, long *, long [], float [],void* ); +/*Cost function prototype*/ +void fakesimul(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs); +typedef struct { + int priorn; + int counter; + double* Gs; + double* Xs; + double* Js; +} Data; + +/*}}}*/ +void M1qn3Usage(void){/*{{{*/ + _printf0_(" usage:\n"); + _printf0_(" X=M1qn3(Xs,Gs);\n"); + _printf0_(" where:\n"); + _printf0_(" Xs are the X values (m x n, where m is the number of independents, n the number of evaluations previously carried out on X)\n"); + _printf0_(" Gs are the G (gradient) values (m x n, where m is the number of independents, n the number of evaluations previously carried out on X,G)\n"); + _printf0_(" X - the new direction.\n"); + _printf0_("\n"); +}/*}}}*/ +WRAPPER(M1qn3){ + + /*input: */ + double* Xs=NULL; + double* Gs=NULL; + double* Js=NULL; + int intn; + int priorn; + int maxsteps,maxiter; + Data data_struct; + + /* output: */ + double* X_out = NULL; + + /*intermediary: */ + double* G=NULL; + double* X=NULL; + double f; + double dxmin=.01; + double gttol=.0001; + long omode; + + /*Boot module: */ + MODULEBOOT(); + + /*checks on arguments on the matlab side: */ + if(nlhs!=NLHS){ + M1qn3Usage(); + _error_("M1qn3 usage error"); + } + if(nrhs!=5){ + M1qn3Usage(); + _error_("M1qn3 usage error"); + } + + /*Input datasets: */ + FetchData(&Xs,&intn,&priorn,XHANDLE); + FetchData(&Gs,&intn,&priorn,GHANDLE); + FetchData(&Js,&priorn,JHANDLE); + FetchData(&maxsteps,MAXSTEPSHANDLE); + FetchData(&maxiter,MAXITERHANDLE); + + /*_printf_("Xs: " << Xs[0] << "\n"); + _printf_("Gs: " << Gs[0] << "\n"); + _printf_("Js: " << Js[0] << "\n"); + _printf_("maxiter: " << maxiter << "\n"); + _printf_("maxsteps: " << maxsteps << "\n");*/ + + /*Initialize M1QN3 parameters*/ + SimulFunc costfuncion = &fakesimul; /*Cost function address*/ + void** prosca = &euclid_; /*Dot product function (euclid is the default)*/ + char normtype[] = "dfn"; /*Norm type: dfn = scalar product defined by prosca*/ + long izs[5]; /*Arrays used by m1qn3 subroutines*/ + long iz[5]; /*Integer m1qn3 working array of size 5*/ + float rzs[1]; /*Arrays used by m1qn3 subroutines*/ + long impres = 0; /*verbosity level*/ + long imode[3] = {0}; /*scaling and starting mode, 0 by default*/ + long indic = 4; /*compute f and g*/ + long reverse = 0; /*reverse or direct mode*/ + long io = 6; /*Channel number for the output*/ + + /*Optimization criterions*/ + long niter = long(maxsteps); /*Maximum number of iterations*/ + long nsim = long(maxiter);/*Maximum number of function calls*/ + + /*Get problem dimension and initialize X, G and f: */ + /*_printf_("intn: " << intn << "\n"); + _printf_("priorn: " << priorn << "\n");*/ + long n = long(intn); + IssmPDouble* G = xNew(n); for (int i=0;i(n); for (int i=0;i(ndz); + double f1=f; + + /*Initialize: */ + data_struct.priorn=priorn; + data_struct.counter=0; + data_struct.Gs=Gs; + data_struct.Js=Js; + data_struct.Xs=Xs; + + m1qn3_(costfuncion,prosca,&ctonbe_,&ctcabe_, + &n,X,&f,G,&dxmin,&f1, + >tol,normtype,&impres,&io,imode,&omode,&niter,&nsim,iz,dz,&ndz, + &reverse,&indic,izs,rzs,(void*)&data_struct); + + switch(int(omode)){ + case 0: /*_printf0_(" Stop requested (indic = 0)\n"); */ break; + case 1: _printf0_(" Convergence reached (gradient satisfies stopping criterion)\n"); break; + case 2: _printf0_(" Bad initialization\n"); break; + case 3: _printf0_(" Line search failure\n"); break; + case 4: _printf0_(" Maximum number of iterations exceeded\n");break; + case 5: _printf0_(" Maximum number of function calls exceeded\n"); break; + case 6: _printf0_(" stopped on dxmin during line search\n"); break; + case 7: _printf0_(" > 0 or <0\n"); break; + default: _printf0_(" Unknown end condition\n"); + } + + /*build output: */ + X_out=xNew(n); + for(int i=0;iXs; + double* Gs=ds->Gs; + double* Js=ds->Js; + + /*Are we done? : */ + if(ds->counter+1==ds->priorn){ + *indic=0; + return; + } + else{ + //_printf0_("counter: " << ds->counter << "\n"); + ds->counter++; + *pf=Js[ds->counter]; + for(int i=0;i<*n;i++)X[i]=Xs[ds->priorn*i+ds->counter]; + for(int i=0;i<*n;i++)G[i]=Gs[ds->priorn*i+ds->counter]; + } + +} Index: /issm/trunk-jpl/src/wrappers/M1qn3/M1qn3.h =================================================================== --- /issm/trunk-jpl/src/wrappers/M1qn3/M1qn3.h (revision 19011) +++ /issm/trunk-jpl/src/wrappers/M1qn3/M1qn3.h (revision 19011) @@ -0,0 +1,53 @@ +/*!\file M1qn3.h + * \brief: prototype for Data Interpolation mex module. + */ + +#ifndef _M1QN3_WRAPPER_H +#define _M1QN3_WRAPPER_H + +#ifdef HAVE_CONFIG_H + #include +#else + #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +/*For python modules: needs to come before header files inclusion*/ +#ifdef _HAVE_PYTHON_ +#define PY_ARRAY_UNIQUE_SYMBOL PythonIOSymbol +#endif + +#include "../bindings.h" +#include "../../c/main/globals.h" +#include "../../c/modules/modules.h" +#include "../../c/shared/shared.h" + +#undef __FUNCT__ +#define __FUNCT__ "M1qn3" + +#ifdef _HAVE_MATLAB_MODULES_ +/* serial input macros: */ +#define XHANDLE prhs[0] +#define GHANDLE prhs[1] +#define JHANDLE prhs[2] +#define MAXSTEPSHANDLE prhs[3] +#define MAXITERHANDLE prhs[4] +/* serial output macros: */ +#define XOUT (mxArray**)&plhs[0] +#endif + +#ifdef _HAVE_PYTHON_MODULES_ +/* serial input macros: */ +#define XHANDLE PyTuple_GetItem(args,0) +#define GHANDLE PyTuple_GetItem(args,1) +#define JHANDLE PyTuple_GetItem(args,2) +#define MAXSTEPSHANDLE PyTuple_GetItem(args,3) +#define MAXITERHANDLE PyTuple_GetItem(args,4) +/* serial output macros: */ +#define XOUT output,0 +#endif + +/* serial arg counts: */ +#undef NLHS +#define NLHS 1 + +#endif /* _M1QN3_WRAPPER_H */ Index: /issm/trunk-jpl/src/wrappers/matlab/Makefile.am =================================================================== --- /issm/trunk-jpl/src/wrappers/matlab/Makefile.am (revision 19010) +++ /issm/trunk-jpl/src/wrappers/matlab/Makefile.am (revision 19011) @@ -38,4 +38,5 @@ EnumToString.la\ ExpSimplify.la\ + ExpToLevelSet.la\ InterpFromGridToMesh.la\ InterpFromMeshToMesh2d.la\ @@ -46,4 +47,5 @@ Ll2xy.la\ NodeConnectivity.la\ + M1qn3.la\ MeshPartition.la\ MeshProfileIntersection.la\ @@ -141,4 +143,9 @@ ContourToMesh_la_LIBADD = ${deps} $(MPILIB) $(PETSCLIB) $(MULTITHREADINGLIB) $(GSLLIB) $(PROJ4LIB) +ExpToLevelSet_la_SOURCES = ../ExpToLevelSet/ExpToLevelSet.cpp\ + ../ExpToLevelSet/ExpToLevelSet.h +ExpToLevelSet_la_LIBADD = ${deps} $(MPILIB) $(PETSCLIB) $(MULTITHREADINGLIB) $(GSLLIB) $(PROJ4LIB) + + ContourToNodes_la_SOURCES = ../ContourToNodes/ContourToNodes.cpp\ ../ContourToNodes/ContourToNodes.h @@ -221,4 +228,8 @@ MeshPartition_la_LIBADD = ${deps} $(MPILIB) $(PETSCLIB) $(METISLIB) $(GSLLIB) $(PROJ4LIB) +M1qn3_la_SOURCES = ../M1qn3/M1qn3.cpp\ + ../M1qn3/M1qn3.h +M1qn3_la_LIBADD = ${deps} $(MPILIB) $(PETSCLIB) $(METISLIB) $(GSLLIB) $(PROJ4LIB) + MeshProfileIntersection_la_SOURCES = ../MeshProfileIntersection/MeshProfileIntersection.cpp\ ../MeshProfileIntersection/MeshProfileIntersection.h