Changeset 21912

Timestamp:

08/03/17 16:05:42 (8 years ago)

Author:

Mathieu Morlighem

Message:

CHG: fixed M1qn3

Location:

issm/trunk-jpl/src/wrappers/M1qn3

Files:

• M1qn3.cpp (modified) (6 diffs)
• M1qn3.h (modified) (2 diffs)

issm/trunk-jpl/src/wrappers/M1qn3/M1qn3.cpp

@@ -48 +48 @@
         int     intn;
         int     priorn;
-        int     maxsteps,maxiter;
         Data    data_struct;
 
@@ -56 +55 @@
         /*intermediary: */
         double f;
-        double dxmin=.01;
-        double gttol=.0001;
+        double dxmin=0.01;
+        double gttol=0.0001;
         long   omode;
 
         /*checks on arguments on the matlab side: */
-        if(nlhs!=NLHS){
-                M1qn3Usage();
-                _error_("M1qn3 usage error");
-        }
-        if(nrhs!=5){
+        if(nlhs!=NLHS || nrhs!=3){
                 M1qn3Usage();
                 _error_("M1qn3 usage error");
@@ -74 +69 @@
         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*/
@@ -97 +84 @@
 
         /*Optimization criterions*/
-        long niter = long(maxsteps); /*Maximum number of iterations*/
-        long nsim  = long(maxiter);/*Maximum number of function calls*/
+        long niter = long(intn); /*Maximum number of iterations*/
+        long nsim  = long(intn); /*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<IssmPDouble>(n); for (int i=0;i<n;i++)G[i]=Gs[i*priorn];
-        IssmPDouble* X = xNew<IssmPDouble>(n); for (int i=0;i<n;i++)X[i]=Xs[i*priorn];
+        IssmPDouble* G = xNew<IssmPDouble>(n); for (int i=0;i<n;i++)G[i]=Gs[i*priorn+0];
+        IssmPDouble* X = xNew<IssmPDouble>(n); for (int i=0;i<n;i++)X[i]=Xs[i*priorn+0];
         f = Js[0];
         
-        /*_printf_("X: " << X[0] << "\n");
-        _printf_("G: " << G[0] << "\n");
-        _printf_("J: " << f << "\n");
-        _printf_("n: " << n << "\n");*/
-
-
         /*Allocate m1qn3 working arrays (see doc)*/
         long      m   = 100;
@@ -133 +112 @@
 
         switch(int(omode)){
-                case 0: /*_printf0_("   Stop requested (indic = 0)\n"); */ break;
+                case 0:  _printf0_("   Stop requested\n"); break;
                 case 1:  _printf0_("   Convergence reached (gradient satisfies stopping criterion)\n"); break;
                 case 2:  _printf0_("   Bad initialization\n"); break;
@@ -180 +159 @@
         }
         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];
+                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];
         }
 
+        printf("J [i] = %g\n",*pf);
 }
 #endif

issm/trunk-jpl/src/wrappers/M1qn3/M1qn3.h

@@ -30 +30 @@
 #define GHANDLE       prhs[1]
 #define JHANDLE       prhs[2]
-#define MAXSTEPSHANDLE prhs[3]
-#define MAXITERHANDLE prhs[4]
 /* serial output macros: */
 #define XOUT (mxArray**)&plhs[0]
@@ -39 +37 @@
 /* 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)
+#define GHANDLE   PyTuple_GetItem(args,1)
+#define JHANDLE   PyTuple_GetItem(args,2)
 /* serial output macros: */
 #define XOUT output,0