Changeset 13195

Timestamp:

08/30/12 11:35:29 (13 years ago)

Author:

utke

Message:

NEW - to trace the solver call for ADOLC call the wrapped solver through the EDF interface

Location:

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

Files:

• Solverx.cpp (modified) (1 diff)
• Solverx.h (modified) (1 diff)
• SolverxSeq.cpp (modified) (3 diffs)

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

@@ -41 +41 @@
                 #endif
                 case SeqMatType:{
-                        SolverxSeq(&uf->svector,Kff->smatrix,pf->svector);
+                        SolverxSeq(&uf->svector,Kff->smatrix,pf->svector,parameters);
                         break;}
                 default:

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

@@ -23 +23 @@
 #endif
 
-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);
+void SolverxSeq(SeqVec** puf,SeqMat* Kff, SeqVec* pf,Parameters* parameters);
+void SolverxSeq(IssmPDouble** pX,IssmPDouble* A,IssmPDouble* B,int n);
 
 #ifdef _HAVE_ADOLC_
+void SolverxSeq(IssmDouble** pX,IssmDouble* A,IssmDouble* B,int n, Parameters* parameters);
 ADOLC_ext_fct EDF_for_solverx;
 #endif

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

@@ -19 +19 @@
 #endif
 
-void SolverxSeq(SeqVec** puf,SeqMat* Kff, SeqVec* pf){/*{{{*/
+#ifdef _HAVE_ADOLC_
+#include "../../shared/Numerics/adolc_edf.h"
+#endif
+
+void SolverxSeq(SeqVec** puf,SeqMat* Kff, SeqVec* pf, Parameters* parameters){/*{{{*/
 
         #ifdef _HAVE_GSL_
@@ -33 +37 @@
         if(M!=N)_error_("Stiffness matrix should be square!");
 
+#ifdef _HAVE_ADOLC_
+        SolverxSeq(&x,Kff->matrix,pf->vector,N,parameters);
+#else
         SolverxSeq(&x,Kff->matrix,pf->vector,N);
+#endif
         uf=new SeqVec(x,N);
 
@@ -52 +60 @@
 }
 
-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;
-        IssmPDouble* pdoubleA=NULL;
-        IssmPDouble* pdoubleB=NULL;
-        IssmPDouble* pdoubleX=NULL;
-
-        /*First, transfer from IssmDouble to double all our matrices and vectors: */
-        pdoubleA=xNew<IssmPDouble>(n*n);
-        pdoubleB=xNew<IssmPDouble>(n);
-        for(i=0;i<n*n;i++)pdoubleA[i]=reCast<IssmPDouble>(A[i]);
-        for(i=0;i<n;i++)pdoubleB[i]=reCast<IssmPDouble>(B[i]);
-        
-        /*Call wrapped solver: */
-        SolverxSeq(&pdoubleX,pdoubleA, pdoubleB, n);
-
-        /*Transfer solution vector from double to IssmDouble: */
-        X = xNew<IssmDouble>(n);
-        for(i=0;i<n;i++)X[i]=reCast<IssmDouble>(pdoubleX[i]);
-
-        /*Free ressources:*/
-        xDelete<IssmPDouble>(pdoubleA);
-        xDelete<IssmPDouble>(pdoubleB);
-
-        /*Assign output pointers: */
-        *pX=X;
+void SolverxSeq(IssmDouble** pX,IssmDouble* A,IssmDouble* B,int n, Parameters* parameters){//{{{
+        // pack inputs to conform to the EDF-prescribed interface
+        IssmDouble*  adoubleEDF_X=xNew<IssmDouble>(n*(n+1)); // packed inputs, i.e. matrix and right hand side
+        for(int i=0; i<n*n;i++)adoubleEDF_X[i]    =A[i]; // pack matrix
+        for(int i=0; i<n;  i++)adoubleEDF_X[i+n*n]=B[i]; // pack the right hand side
+        IssmPDouble* pdoubleEDF_X=xNew<IssmPDouble>(n*(n+1)); // provide space to transfer inputs during call_ext_fct
+        IssmPDouble* pdoubleEDF_Y=xNew<IssmPDouble>(n);       // provide space to transfer outputs during call_ext_fct
+        // call the wrapped solver through the registry entry we retrieve from parameters
+        call_ext_fct(dynamic_cast<GenericParam<Adolc_edf> * >(parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p,
+                     n*(n+1), pdoubleEDF_X, adoubleEDF_X,
+                     n, pdoubleEDF_Y, B);
+        xDelete(adoubleEDF_X);
+        xDelete(pdoubleEDF_X);
+        xDelete(pdoubleEDF_Y);
 }
 /*}}}*/