Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 13184) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 13185) @@ -136,4 +136,5 @@ ./classes/matrix/Vector.cpp\ ./classes/objects/Params/Param.h\ + ./classes/objects/Params/GenericParam.h\ ./classes/objects/Params/BoolParam.cpp\ ./classes/objects/Params/BoolParam.h\ @@ -196,4 +197,5 @@ ./shared/Matrix/matrix.h\ ./shared/Matrix/MatrixUtils.cpp\ + ./shared/Numerics/adolc_edf.h\ ./shared/Numerics/numerics.h\ ./shared/Numerics/Verbosity.h\ Index: /issm/trunk-jpl/src/c/classes/objects/Params/GenericParam.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Params/GenericParam.h (revision 13185) +++ /issm/trunk-jpl/src/c/classes/objects/Params/GenericParam.h (revision 13185) @@ -0,0 +1,102 @@ +/* + * GenericParam.h + * + * Created on: Aug 29, 2012 + * Author: utke + */ + +#ifndef GENERICPARAM_H_ +#define GENERICPARAM_H_ + +/*Headers:*/ +/*{{{*/ +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "./Param.h" +#include "../../../include/include.h" +#include "../../../shared/shared.h" +/*}}}*/ + +/** + * here we have a class that holds an instance of P + * but because it should live side by side with + * the other instances derived from Param it - unfortunately - + * inherits all the accessors that are useless in this context + */ +template class GenericParam: public Param{ + + private: + P myP; + int myEnumVal; + + public: + /*GenericParam constructors, destructors: {{{*/ + GenericParam(int enumVal) : myEnumVal(enumVal){}; + ~GenericParam(){}; + /*}}}*/ + /*Object virtual functions definitions:{{{ */ + // unfortunately having such a printer method implies + // that P must provide the friend << operator + void DeepEcho() { + _printLine_("GenericParam:"); + _printLine_(" enum: " << myEnumVal << " (" << EnumToStringx(myEnumVal) << ")"); + _printLine_(" value: " << myP);; + } + void Echo() {DeepEcho();}; + int Id(){ return -1; }; + int MyRank() { extern int my_rank; return my_rank;} ; + int ObjectEnum() {return AdolcParamEnum;}; + + // the "copy" has to implement the base class abstract function + // but I would prefer to drop this not to hide a "new" in here because + // it does not clarify ownership of the newed up instance... + // use the default copy constructor instead + Object* copy() { return new GenericParam

(*this); }; + /*}}}*/ + /*Param vritual function definitions: {{{*/ + int InstanceEnum(){return myEnumVal;} + void GetParameterName(char**pname) {EnumToStringx(pname,this->myEnumVal);} + + P& GetParameterValue() { return myP;} + const P& GetParameterValue()const { return myP;}; + + // none of these apply ... + void GetParameterValue(bool* pbool){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot return a bool");} + void GetParameterValue(int* pinteger){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot return an integer");} + void GetParameterValue(int** pintarray,int* pM){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot return an array of integers");} + void GetParameterValue(int** pintarray,int* pM,int* pN){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot return an array of integers");} + void GetParameterValue(IssmDouble* pIssmDouble){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot return a IssmDouble");} + void GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot return a IssmDouble for a given time");} + void GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot return a string");} + void GetParameterValue(char*** pstringarray,int* pM){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot return a string array");} + void GetParameterValue(IssmDouble** pIssmDoublearray,int* pM){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot return a IssmDouble array");} + void GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, int* pN){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot return a IssmDouble array");} + void GetParameterValue(IssmDouble*** parray, int* pM,int** pmdims, int** pndims){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot return a matrix array");} + void GetParameterValue(Vector** pvec){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot return a Vec");} + void GetParameterValue(Matrix** pmat){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot return a Mat");} + void GetParameterValue(FILE** pfid){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot return a FILE");} + + void SetValue(bool boolean){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot hold a bool");} + void SetValue(int integer){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot hold an integer");} + void SetValue(int* intarray,int M){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot hold an int array");} + void SetValue(int* intarray,int M,int N){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot hold an int array");} + void SetValue(IssmDouble scalar){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot hold an IssmDouble");} + void SetValue(char* string){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot hold a string");} + void SetValue(char** stringarray,int M){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot hold a string array");} + void SetValue(IssmDouble* IssmDoublearray,int M){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot hold a IssmDouble array");} + void SetValue(IssmDouble* pIssmDoublearray,int M,int N){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot hold a IssmDouble array");} + void SetValue(Vector* vec){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot hold a Vec");} + void SetValue(Matrix* mat){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot hold a Mat");} + void SetValue(FILE* fid){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot hold a FILE");} + void SetValue(IssmDouble** array, int M, int* mdim_array, int* ndim_array){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot hold an array of matrices");} + void UnitConversion(int direction_enum) {/* nothing useful here either */}; + + /*}}}*/ +}; + + +#endif /* GENERICPARAM_H_ */ Index: /issm/trunk-jpl/src/c/classes/objects/objects.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/objects.h (revision 13184) +++ /issm/trunk-jpl/src/c/classes/objects/objects.h (revision 13185) @@ -112,4 +112,5 @@ /*Params: */ +#include "./Params/GenericParam.h" #include "./Params/BoolParam.h" #include "./Params/DoubleMatParam.h" Index: /issm/trunk-jpl/src/c/modules/Solverx/Solverx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/Solverx/Solverx.h (revision 13184) +++ /issm/trunk-jpl/src/c/modules/Solverx/Solverx.h (revision 13185) @@ -27,4 +27,8 @@ void SolverxSeq(double** pX,double* A,double* B,int n); +#ifdef _HAVE_ADOLC_ +ADOLC_ext_fct EDF_for_solverx; +#endif + #endif /* _SOLVERX_H */ Index: /issm/trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp (revision 13184) +++ /issm/trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp (revision 13185) @@ -45,4 +45,11 @@ }/*}}}*/ #ifdef _HAVE_ADOLC_ + +int EDF_for_solverx(int n, IssmPDouble *x, int m, IssmPDouble *y) { + if(m*(m+1)!=n)_error_("Stiffness matrix should be square!"); + SolverxSeq(&y,x, x+m*m, m); + return 0; +} + void SolverxSeq(IssmDouble** pX,IssmDouble* A,IssmDouble* B,int n){//{{{ /* if we use Adol-C then the IssmDouble will be an adouble @@ -61,6 +68,6 @@ /*First, transfer from IssmDouble to double all our matrices and vectors: */ - pdoubleA=xNew(n*n); - pdoubleB=xNew(n); + pdoubleA=xNew(n*n); + pdoubleB=xNew(n); for(i=0;i(A[i]); for(i=0;i(B[i]); Index: /issm/trunk-jpl/src/c/shared/Numerics/adolc_edf.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Numerics/adolc_edf.h (revision 13185) +++ /issm/trunk-jpl/src/c/shared/Numerics/adolc_edf.h (revision 13185) @@ -0,0 +1,28 @@ +/* + * adolc_edf.h + * + * Created on: Jun 26, 2012 + * Author: utke + */ + +#ifndef _ADOLC_EDF_H_ +#define _ADOLC_EDF_H_ + +#include "../../include/include.h" + +#ifdef _HAVE_ADOLC_ + +struct Adolc_edf { + ext_diff_fct *myEDF_for_solverx_p; + Adolc_edf() : myEDF_for_solverx_p(0) {} + inline friend std::ostream& operator << ( ostream&, const Adolc_edf& ); +}; + +std::ostream& operator << ( std::ostream& out, const Adolc_edf& a) { + out << a.myEDF_for_solverx_p; + return out; +} + +#endif + +#endif Index: /issm/trunk-jpl/src/c/solutions/issm.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutions/issm.cpp (revision 13184) +++ /issm/trunk-jpl/src/c/solutions/issm.cpp (revision 13185) @@ -5,4 +5,5 @@ #include "../issm.h" #include "../include/globals.h" +#include "../shared/Numerics/adolc_edf.h" int main(int argc,char **argv){ @@ -119,4 +120,7 @@ /*If running AD, then initialize the tape: */ #ifdef _HAVE_ADOLC_ + GenericParam *theAdolcEDF_p=new GenericParam(AdolcParamEnum); + theAdolcEDF_p->GetParameterValue().myEDF_for_solverx_p=reg_ext_fct(EDF_for_solverx); + femmodel->parameters->AddObject(theAdolcEDF_p); if(autodiff) trace_on(1); #endif