Changeset 12988

Timestamp:

08/10/12 17:09:30 (13 years ago)

Author:

Eric.Larour

Message:

NEW: created wrapper to SolverxSeq solver, so that when running with Adolc support,
we can go from solving for a IssmDouble* vector solution to solving a double* vector
solution, and come back to an IssmDouble* vector.

In addition, added tape start phase and tape statistics at the beginning and end of issm.cpp
respectively.

Also fixed a bug (reCast type) in PddSurfaceMassBalance.

Fixed a bug in solve.m, where cluster.np was set to 1 too late for Adolc
runs.

config-macosx64-larour-ad.sh now defaults to using adolc.

Location:

issm/trunk-jpl

Files:

• configs/config-macosx64-larour-ad.sh (modified) (1 diff)
• src/c/classes/objects/Params/TransientParam.cpp (modified) (1 diff)
• src/c/modules/Solverx/Solverx.h (modified) (1 diff)
• src/c/modules/Solverx/SolverxSeq.cpp (modified) (1 diff)
• src/c/shared/Elements/PddSurfaceMassBalance.cpp (modified) (2 diffs)
• src/c/solutions/issm.cpp (modified) (3 diffs)
• src/m/model/solve.m (modified) (2 diffs)

issm/trunk-jpl/configs/config-macosx64-larour-ad.sh

@@ -17 +17 @@
         --without-groundingline \
         --without-kriging  \
-        --with-gsl-dir=$ISSM_DIR/externalpackages/gsl/install
-        #--with-adolc-dir=$ISSM_DIR/externalpackages/adolc/install\
+        --with-gsl-dir=$ISSM_DIR/externalpackages/gsl/install \
+        --with-adolc-dir=$ISSM_DIR/externalpackages/adolc/install

issm/trunk-jpl/src/c/classes/objects/Params/TransientParam.cpp

@@ -99 +99 @@
 void  TransientParam::GetParameterValue(IssmDouble* pdouble,IssmDouble time){
 
-        double output;
+        IssmDouble output;
         bool   found;
 

issm/trunk-jpl/src/c/modules/Solverx/Solverx.h

@@ -25 +25 @@
 void SolverxSeq(SeqVec** puf,SeqMat* Kff, SeqVec* pf);
 void SolverxSeq(IssmDouble** pX,IssmDouble* A,IssmDouble* B,int n);
+void SolverxSeq(double** pX,double* A,double* B,int n);
 
 #endif  /* _SOLVERX_H */

issm/trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp

@@ -44 +44 @@
 
 }/*}}}*/
-void SolverxSeq(IssmDouble** pX,IssmDouble* A,IssmDouble* B,int n){/*{{{*/
+#ifdef _HAVE_ADOLC_
+void SolverxSeq(IssmDouble** pX,IssmDouble* A,IssmDouble* B,int n){//{{{
+        /* if we use Adol-C then the IssmDouble will be an adouble
+           and the calls to gsl_... will not work. 
+           We therefore call a wrapped solver instead. 
+        */
 
+        /*Output: */
+        IssmDouble* X=NULL;
+
+        /*Intermediary: */
+        int     i;
+        double* doubleA=NULL;
+        double* doubleB=NULL;
+        double* doubleX=NULL;
+
+        /*First, transfer from IssmDouble to double all our matrices and vectors: */
+        doubleA=xNew<double>(n*n);
+        doubleB=xNew<double>(n);
+        for(i=0;i<n*n;i++)A[i]>>=doubleA[i];
+        for(i=0;i<n;i++)B[i]>>=doubleB[i];
+        
+        /*Call wrapped solver: */
+        SolverxSeq(&doubleX,doubleA, doubleB, n);
+
+        /*Transfer solution vector from double to IssmDouble: */
+        X = xNew<IssmDouble>(n);
+        for(i=0;i<n;i++)X[i]<<=doubleX[i];
+
+        /*Free ressources:*/
+        xDelete<double>(doubleA);
+        xDelete<double>(doubleB);
+
+        /*Assign output pointers: */
+        *pX=X;
+}
+/*}}}*/
+#endif
+//void SolverxSeq(double** pX,double* A,double* B,int n){ //{{{
+#ifdef _HAVE_ADOLC_
+void SolverxSeq(double** pX,double* A,double* B,int n){ 
+#else
+void SolverxSeq(IssmDouble** pX,IssmDouble* A,IssmDouble* B,int n){
+#endif
         #ifdef _HAVE_GSL_
-                /*GSL Matrices and vectors: */
-                int              s;
-                gsl_matrix_view  a;
-                gsl_vector_view  b;
-                gsl_vector      *x = NULL;
-                gsl_permutation *p = NULL;
-                #ifdef _HAVE_ADOLC_
-                        // if we use Adol-C then the IssmDouble will be an adouble
-                        // and the calls to gsl_... will not work
-                        // and we should call  a suitable wrapped solve instead
-                        _error2_("SolverxSeq: should not be here with Adol-C");
-                #else
-                        /*A will be modified by LU decomposition. Use copy*/
-                        IssmDouble* Acopy = xNew<IssmDouble>(n*n);
-                        xMemCpy<IssmDouble>(Acopy,A,n*n);
+        /*GSL Matrices and vectors: */
+        int              s;
+        gsl_matrix_view  a;
+        gsl_vector_view  b;
+        gsl_vector      *x = NULL;
+        gsl_permutation *p = NULL;
+        /*A will be modified by LU decomposition. Use copy*/
+        double* Acopy = xNew<double>(n*n);
+        xMemCpy<double>(Acopy,A,n*n);
 
-                        /*Initialize gsl matrices and vectors: */
-                        a = gsl_matrix_view_array (Acopy,n,n);
-                        b = gsl_vector_view_array (B,n);
-                        x = gsl_vector_alloc (n);
+        /*Initialize gsl matrices and vectors: */
+        a = gsl_matrix_view_array (Acopy,n,n);
+        b = gsl_vector_view_array (B,n);
+        x = gsl_vector_alloc (n);
 
-                        /*Run LU and solve: */
-                        p = gsl_permutation_alloc (n);
-                        gsl_linalg_LU_decomp (&a.matrix, p, &s);
-                        gsl_linalg_LU_solve (&a.matrix, p, &b.vector, x);
+        /*Run LU and solve: */
+        p = gsl_permutation_alloc (n);
+        gsl_linalg_LU_decomp (&a.matrix, p, &s);
+        gsl_linalg_LU_solve (&a.matrix, p, &b.vector, x);
 
-                        //printf ("x = \n");
-                        //gsl_vector_fprintf (stdout, x, "%g");
+        //printf ("x = \n");
+        //gsl_vector_fprintf (stdout, x, "%g");
 
-                        /*Copy result*/
-                        IssmDouble* X = xNew<IssmDouble>(n);
-                        memcpy(X,gsl_vector_ptr(x,0),n*sizeof(IssmDouble));
+        /*Copy result*/
+        double* X = xNew<double>(n);
+        memcpy(X,gsl_vector_ptr(x,0),n*sizeof(double));
 
-                        /*Clean up and assign output pointer*/
-                        xDelete<IssmDouble>(Acopy);
-                        gsl_permutation_free(p);
-                        gsl_vector_free(x);
-                        *pX=X;
-                #endif
+        /*Clean up and assign output pointer*/
+        xDelete<double>(Acopy);
+        gsl_permutation_free(p);
+        gsl_vector_free(x);
+        *pX=X;
         #endif
-}/*}}}*/
+}
+/*}}}*/

issm/trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp

@@ -80 +80 @@
       if (tstar < PDup){
         pd = 1.;
-        if (tstar >= -siglimc){ pd = pds[int(tstar/DT + siglim0c)];}}
+        if (tstar >= -siglimc){ pd = pds[reCast<int,IssmDouble>(tstar/DT + siglim0c)];}}
       else { 
         pd = 0.;}
@@ -100 +100 @@
       if (tstar >= siglim) {pdd = pdd + tstar*deltm;}
       else if (tstar> -siglim){
-        pddsig=pdds[int(tstar/DT + siglim0)];
+        pddsig=pdds[reCast<int,IssmDouble>(tstar/DT + siglim0)];
         pdd = pdd + pddsig*deltm;
         frzndd = frzndd - (tstar-pddsig)*deltm;}

issm/trunk-jpl/src/c/solutions/issm.cpp

@@ -14 +14 @@
         bool  dakota_analysis,control_analysis,tao_analysis;
 
+        /*AD: */
+        int   tape_stats[11];
+
         /*FemModel: */
         FemModel *femmodel = NULL;
@@ -41 +44 @@
 
         ISSMBOOT();
+
+        /*If running AD, then initialize the tape: */
+        #ifdef _HAVE_ADOLC_
+        trace_on(1);
+        #endif
 
         /*Initialize environments: Petsc, MPI, etc...: */
@@ -200 +208 @@
         #endif
         #endif
+
+        /*If running AD, close our tape, print statistics: */
+        #ifdef _HAVE_ADOLC_
+        trace_off();
+        tapestats(1,tape_stats); //reading of tape statistics
+        _pprintLine_("   ADOLC statistics: ");
+        _pprintLine_("   "<<setw(45)<<left<<"Number of independents: " <<tape_stats[0]);
+        _pprintLine_("   "<<setw(45)<<left<<"Number of dependents: " <<tape_stats[1]);
+        _pprintLine_("   "<<setw(45)<<left<<"Maximal number of live active variables: " <<tape_stats[2]);
+        _pprintLine_("   "<<setw(45)<<left<<"Size of value stack (number of overwrites): " <<tape_stats[3]);
+        _pprintLine_("   "<<setw(45)<<left<<"Buffer size (a multiple of eight): " <<tape_stats[4]);
+        _pprintLine_("   "<<setw(45)<<left<<"Total number of operations recorded: " <<tape_stats[5]);
+        #endif
         
         /*end module: */

issm/trunk-jpl/src/m/model/solve.m

@@ -68 +68 @@
 end
 
+%we need to make sure we have PETSC support, otherwise, we run with only one cpu: 
+if ~ispetsc,
+        disp('PETSC support not included, running on 1 cpu only!');
+        cluster.np=1;
+end
+
+
 %Wite all input files
 marshall(md);                                          % bin file
@@ -73 +80 @@
 BuildQueueScript(cluster,md.private.runtimename,md.miscellaneous.name,md.private.solution,md.settings.io_gather,md.debug.valgrind,md.debug.gprof); % queue file
 
-%we need to make sure we have PETSC support, otherwise, we run with only one cpu: 
-if ~ispetsc,
-        disp('PETSC support not included, running on 1 cpu only!');
-        cluster.np=1;
-end
 
 %Stop here if batch mode