Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 20652) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 20653) @@ -26,12 +26,11 @@ ./classes/gauss/GaussTetra.cpp\ ./classes/gauss/GaussPenta.cpp\ + ./classes/IoModel.cpp\ ./classes/FemModel.cpp\ ./classes/Loads/Friction.cpp\ ./classes/Inputs/TransientInput.cpp\ ./classes/Constraints/SpcTransient.cpp\ - ./classes/IndependentObject.cpp\ ./classes/DependentObject.cpp\ ./classes/DofIndexing.cpp\ - ./classes/IoModel.cpp\ ./classes/Contours.cpp\ ./classes/Nodes.cpp\ Index: /issm/trunk-jpl/src/c/classes/IoModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/IoModel.cpp (revision 20652) +++ /issm/trunk-jpl/src/c/classes/IoModel.cpp (revision 20653) @@ -3,4 +3,12 @@ * into ISSM, from Matlab, or through a binary file opened for reading. */ + +/*CODES: + * 1: boolean constant + * 2: integer constant + * 3: IssmDouble constant + * 5: boolean vector + * 6: int vector + * 7: IssmDouble vector*/ #ifdef HAVE_CONFIG_H @@ -19,11 +27,364 @@ #include "../shared/shared.h" +/*NEW NEW NEW NEW NEW NEW NEW NEW NEW NEW*/ +IoConstant::IoConstant(){/*{{{*/ + this->isindependent = false; + this->data_enum = -1; + this->name = NULL; + this->constant = NULL; +} +/*}}}*/ +IoConstant::~IoConstant(){/*{{{*/ + xDelete(this->name); + delete this->constant; +} +/*}}}*/ +IoConstant::IoConstant(bool value,int enum_in){/*{{{*/ + this->isindependent = false; + this->data_enum = enum_in; + this->name = NULL; + this->constant = new BoolParam(enum_in,value); +} +/*}}}*/ +IoConstant::IoConstant(int value,int enum_in){/*{{{*/ + this->isindependent = false; + this->data_enum = enum_in; + this->name = NULL; + this->constant = new IntParam(enum_in,value); +} +/*}}}*/ +IoConstant::IoConstant(IssmDouble value,int enum_in){/*{{{*/ + this->isindependent = false; + this->data_enum = enum_in; + this->name = NULL; + this->constant = new DoubleParam(enum_in,value); +} +/*}}}*/ +IoConstant::IoConstant(char* value,int enum_in){/*{{{*/ + this->isindependent = false; + this->data_enum = enum_in; + this->name = NULL; + this->constant = new StringParam(enum_in,value); +} +/*}}}*/ + +IoData::IoData(){/*{{{*/ + this->isindependent = false; + this->data_enum = -1; + this->name = NULL; + this->code = -1; + this->layout = -1; + this->M = 0; + this->N = 0; + this->data = NULL; +} +/*}}}*/ +IoData::~IoData(){/*{{{*/ + xDelete(this->name); + xDelete(this->data); +} +/*}}}*/ +IoData::IoData(IssmDouble* matrix,int code_in,int layout_in,int M_in,int N_in,int enum_in){/*{{{*/ + this->isindependent = false; + this->data_enum = enum_in; + this->name = NULL; + this->code = code_in; + this->layout = layout_in; + this->M = M_in; + this->N = N_in; + this->data = matrix; /*do not copy*/ + _assert_(code_in==5 || code_in==6 || code_in==7); +} +/*}}}*/ + +void IoModel::AddConstant(IoConstant* in_constant){/*{{{*/ + + _assert_(in_constant); + + /*Go through dataset of constant and check whether it already exists */ + vector::iterator iter; + + for(iter=constants.begin();iterdata_enum==in_constant->data_enum){ + delete in_constant; + return; + } + } + + this->constants.push_back(in_constant); +} +/*}}}*/ +void IoModel::AddConstantIndependent(IoConstant* in_constant){/*{{{*/ + + _assert_(in_constant); + + /*Set constant as independent*/ + in_constant->isindependent = true; + + /*Add to constnats*/ + this->AddConstant(in_constant); +} +/*}}}*/ +void IoModel::AddData(IoData* in_data){/*{{{*/ + + _assert_(in_data); + + /*Go through dataset of data and check whether it already exists */ + vector::iterator iter; + + for(iter=data.begin();iterdata_enum==in_data->data_enum){ + delete in_data; + return; + } + } + + this->data.push_back(in_data); +} +/*}}}*/ +void IoModel::AddDataIndependent(IoData* in_data){/*{{{*/ + + _assert_(in_data); + + /*Set data as independent*/ + in_data->isindependent = true; + + /*Add to constnats*/ + this->AddData(in_data); +} +/*}}}*/ +void IoModel::FetchIndependentConstant(int* pXcount,IssmPDouble* X,int name){/*{{{*/ + + /*recover my_rank:*/ + int my_rank=IssmComm::GetRank(); + + /*recover Xcount if X is not NULL:*/ + int Xcount = 0; + if(X) Xcount=*pXcount; + + #ifdef _HAVE_ADOLC_ //cannot come here unless you are running AD mode, from DeclaredIndependents: + + /*output: */ + IssmPDouble pscalar; + IssmDouble scalar; //same as pscalar, except it's an ADOLC independent variable + int code; + + /*Set file pointer to beginning of the data: */ + fid=this->SetFilePointerToData(&code,NULL,name); + if(code!=3) _error_("expecting a IssmDouble for enum " << EnumToStringx(name)); + + /*We have to read a scalar from disk. First read the dimensions of the scalar, then the scalar: */ + if(my_rank==0){ + if(fread(&pscalar,sizeof(IssmPDouble),1,fid)!=1)_error_("could not read scalar "); + + /*Now, before we even broadcast this to other nodes, declare the scalar as an independent variable!. If we + *have been supplied an X vector, use it instead of what we just read: */ + if(X){ + scalar<<=X[Xcount]; + } + else{ + scalar<<=pscalar; + } + } + + ISSM_MPI_Bcast(&scalar,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + this->AddConstantIndependent(new IoConstant(scalar,name)); + + /*increment offset into X vector, now that we have read 1 value:*/ + Xcount++; *pXcount=Xcount; + #endif +} +/*}}}*/ +void IoModel::FetchIndependentData(int* pXcount,IssmPDouble* X,int name){/*{{{*/ + + /*recover my_rank:*/ + int my_rank=IssmComm::GetRank(); + + /*recover Xcount if X is not NULL:*/ + int Xcount = 0; + if(X) Xcount=*pXcount; + + #ifdef _HAVE_ADOLC_ //cannot come here unless you are running AD mode, from DeclaredIndependents: + + /*Intermediaries*/ + int M,N; + IssmPDouble* buffer=NULL; //a buffer to read the data from disk + IssmDouble* matrix=NULL; //our independent variable + int code,layout; + + /*Set file pointer to beginning of the data: */ + fid=this->SetFilePointerToData(&code,&layout,name); + if((code!=5) && (code!=6) && (code!=7))_error_("expecting a IssmDouble, integer or boolean matrix for enum " << EnumToStringx(name)); + + /*We have to read a matrix from disk. First read the dimensions of the matrix, then the whole matrix: */ + /*numberofelements: */ + if(my_rank==0){ + if(fread(&M,sizeof(int),1,fid)!=1) _error_("could not read number of rows for matrix "); + } + ISSM_MPI_Bcast(&M,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + + if(my_rank==0){ + if(fread(&N,sizeof(int),1,fid)!=1) _error_("could not read number of columns for matrix "); + } + ISSM_MPI_Bcast(&N,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + + /*Now allocate matrix: */ + if(M*N){ + buffer=xNew(M*N); + matrix=xNew(M*N); + + /*Read matrix on node 0, then broadcast: */ + if(my_rank==0){ + if(fread(buffer,M*N*sizeof(IssmPDouble),1,fid)!=1) _error_("could not read matrix "); + + /*Now, before we even broadcast this to other nodes, declare the whole matrix as a independent variable! + If we have been supplied an X vector, use it instead of what we just read: */ + if(X){ + for(int i=0;i(buffer); + } + else _error_("cannot declare the independent variable " << EnumToStringx(name) << "if it's empty!"); + + /*Add to data as independent*/ + this->AddDataIndependent(new IoData(matrix,code,layout,M,N,name)); + + /*increment offset into X vector, now that we have read M*N values:*/ + Xcount+=M*N; *pXcount=Xcount; + #endif +} +/*}}}*/ +void IoModel::FindConstant(bool* pvalue,int constant_enum){/*{{{*/ + + /*Intermediary*/ + vector::iterator iter; + + for(iter=constants.begin();iterdata_enum==constant_enum){ + ioconstant->constant->GetParameterValue(pvalue); + return; + } + } + + for(vector::iterator iter=constants.begin();iterconstant->Echo(); + _error_("Could not find constant \""<::iterator iter; + + for(iter=constants.begin();iterdata_enum==constant_enum){ + ioconstant->constant->GetParameterValue(pvalue); + return; + } + } + + _error_("Could not find constant \""<::iterator iter; + + for(iter=constants.begin();iterdata_enum==constant_enum){ + ioconstant->constant->GetParameterValue(pvalue); + return; + } + } + + _error_("Could not find constant \""<::iterator iter; + + for(iter=constants.begin();iterdata_enum==constant_enum){ + ioconstant->constant->GetParameterValue(pvalue); + return; + } + } + + _error_("Could not find constant \""<::iterator iter=constants.begin();iterisindependent){ + num_independents+= 1; + } + } + + /*Process data*/ + for(vector::iterator iter=data.begin();iterisindependent){ + num_independents+= (*iter)->M*(*iter)->N; + } + } + + /*return*/ + return num_independents; +} +/*}}}*/ +void IoModel::FillIndependents(IssmDouble* xp){/*{{{*/ + + _assert_(xp); + + /*Initialize local num ind*/ + int local_num_ind = 0; + + /*Process constants*/ + for(vector::iterator iter=constants.begin();iterisindependent){ + (*iter)->constant->GetParameterValue(&xp[local_num_ind]); + local_num_ind += 1; + } + } + + /*Process data*/ + for(vector::iterator iter=data.begin();iterisindependent){ + for(int i=0;i<(*iter)->M*(*iter)->N;i++){ + xp[local_num_ind+i] = (*iter)->data[i]; + } + local_num_ind += (*iter)->M*(*iter)->N; + } + } + + _assert_(local_num_ind == this->NumIndependents()); +} +/*}}}*/ +/*NEW NEW NEW NEW NEW NEW NEW NEW NEW NEW*/ + IoModel::IoModel(){/*{{{*/ + this->fid=NULL; this->solution_enum=-1; - this->data=NULL; - this->independents=NULL; - this->independent_objects=NULL; - this->constants=NULL; this->my_elements=NULL; @@ -65,8 +426,4 @@ this->solution_enum = solution_enum_in; - /*Initialize data: */ - this->data=xNew(MaximumNumberOfDefinitionsEnum); - for(int i=0;idata[i]=NULL; - /*If we are running in AD mode, we need to start the trace and declare our independent variables now, *and prevent them from being erased during successive calls to iomodel->FetchConstants, iomodel->FetchData and @@ -76,21 +433,19 @@ /*Initialize and read constants:*/ - this->constants=new Parameters(); - this->FetchConstants(); /*this routine goes through the input file, and fetches bools, ints, IssmDoubles and strings only, nothing memory intensive*/ - - /*Now some very specific piece of logic. We want to know whether we are going to carry out an autodiff run. There are - *several cases: either trace is true (which bypasses the isautodiff choice), or trace is false but isautodiff is on, in which - *case two cases are possible: are we running a control method or not? To make things simpler, we are going to pin down everyting - *on isautodiff. By the way, we could not do this before we had the constants dataset initialized! {{{*/ - this->constants->FindParam(&autodiff,AutodiffIsautodiffEnum); - this->constants->FindParam(&iscontrol,InversionIscontrolEnum); - - if(trace)autodiff=true; + this->FetchConstants(); /*this routine goes through the input file, and fetches bool, int, IssmDouble and string only, nothing memory intensive*/ + + /*Is this an autodiff run?*/ + this->FindConstant(&autodiff,AutodiffIsautodiffEnum); + this->FindConstant(&iscontrol,InversionIscontrolEnum); + if(trace){ + autodiff=true; + } else{ - if(autodiff && !iscontrol) autodiff=true; - else autodiff=false; - } - this->constants->SetParam(autodiff,AutodiffIsautodiffEnum); - /*}}}*/ + if(autodiff && !iscontrol) + autodiff=true; + else + autodiff=false; + } + this->AddConstant(new IoConstant(autodiff,AutodiffIsautodiffEnum)); /*Initialize permanent data: */ @@ -119,30 +474,22 @@ IoModel::~IoModel(){/*{{{*/ - /*Delete independents*/ - if(this->independents){ - for(int i=0;iindependents[i]){ - IssmDouble* array=this->data[i]; - xDelete(array); - } - } - } - - /*checks in debugging mode*/ - #if defined(_ISSM_DEBUG_) && !defined(_HAVE_ADOLC_) - if(this->data){ - for(int i=0;idata[i]){ - _printf0_("WARNING: previous pointer of " << EnumToStringx(i) << " has not been freed (DeleteData has not been called)\n"); - } - } - } - #endif - - if(this->constants) delete this->constants; - - xDelete(this->data); - xDelete(this->independents); - if(this->independent_objects)delete this->independent_objects; + /*Delete constants*/ + vector::iterator iter1; + for(iter1=constants.begin();iter1constants.clear(); + + /*Delete data*/ + vector::iterator iter2; + for(iter2=data.begin();iter2isindependent){ + _printf0_("WARNING: IoData " << EnumToStringx((*iter2)->data_enum) << " has not been freed (DeleteData has not been called)\n"); + } + #endif + delete *iter2; + } + this->data.clear(); xDelete(this->my_elements); @@ -244,54 +591,32 @@ } /*}}}*/ -void IoModel::FindConstant(bool* poutput,int constant_enum){/*{{{*/ - - _assert_(constant_enum>=0); - _assert_(this->constants); - - this->constants->FindParam(poutput,constant_enum); -} -/*}}}*/ -void IoModel::FindConstant(int* poutput,int constant_enum){/*{{{*/ - - _assert_(constant_enum>=0); - _assert_(this->constants); - - this->constants->FindParam(poutput,constant_enum); -} -/*}}}*/ -void IoModel::FindConstant(IssmDouble* poutput,int constant_enum){/*{{{*/ - - _assert_(constant_enum>=0); - _assert_(this->constants); - - this->constants->FindParam(poutput,constant_enum); -} -/*}}}*/ -void IoModel::FindConstant(char** poutput,int constant_enum){/*{{{*/ - - _assert_(constant_enum>=0); - _assert_(this->constants); - - this->constants->FindParam(poutput,constant_enum); -} -/*}}}*/ Param* IoModel::CopyConstantObject(int constant_enum){/*{{{*/ - _assert_(this->constants); - - /*Find constant*/ - //FIXME _error_("Needs to be FIXED"); - Param* param=this->constants->FindParamObject(constant_enum); - if(!param) _error_("Constant " << EnumToStringx(constant_enum) << " not found in iomodel"); - - return xDynamicCast(param->copy()); + /*Intermediary*/ + vector::iterator iter; + + for(iter=constants.begin();iterdata_enum==constant_enum){ + return ioconstant->constant->copy(); + } + } + + _error_("Constant " << EnumToStringx(constant_enum) << " not found in iomodel"); + return NULL; } /*}}}*/ IssmDouble* IoModel::Data(int data_enum){/*{{{*/ - _assert_(data_enum=0); - - return this->data[data_enum]; + /*Intermediary*/ + vector::iterator iter; + + for(iter=data.begin();iterdata_enum==data_enum) return iodata->data; + } + + return NULL; } /*}}}*/ @@ -338,7 +663,5 @@ void IoModel::DeclareIndependents(bool trace,IssmPDouble* X){/*{{{*/ - int i; - bool autodiff = false; - bool iscontrol = false; + bool autodiff,iscontrol; int num_independent_objects; int Xcount=0; @@ -347,10 +670,5 @@ int *types = NULL; - int dummy; - /*Initialize array detecting whether data[i] is an independent AD mode variable: */ - this->independents=xNew(MaximumNumberOfDefinitionsEnum); - for(i=0;iindependents[i]=false; - this->FetchData(&autodiff,AutodiffIsautodiffEnum); this->FetchData(&iscontrol,InversionIscontrolEnum); @@ -359,23 +677,24 @@ #ifdef _HAVE_ADOLC_ - /*build dataset made of independent objects:*/ - this->independent_objects=new DataSet(); this->FetchData(&num_independent_objects,AutodiffNumIndependentObjectsEnum); if(num_independent_objects){ - this->FetchData(&names,&dummy,&dummy,AutodiffIndependentObjectNamesEnum); - this->FetchData(&types,&dummy,&dummy,AutodiffIndependentObjectTypesEnum); + this->FetchData(&names,NULL,NULL,AutodiffIndependentObjectNamesEnum); + this->FetchData(&types,NULL,NULL,AutodiffIndependentObjectTypesEnum); /*create independent objects, and at the same time, fetch the corresponding independent variables, *and declare them as such in ADOLC: */ - for(i=0;iindependent_objects->AddObject(independent_object); - - /*now go fetch the independent variable: */ - independent_object->FetchIndependent(this,&Xcount,X); //supply the pointer to iomodel. + for(int i=0;iFetchIndependentConstant(&Xcount,X,names[i]); + } + else if(types[i]==1){ + /* vector:*/ + this->FetchIndependentData(&Xcount,X,names[i]); + } + else{ + _error_("Independent cannot be of size " << types[i]); + } } xDelete(names); @@ -387,29 +706,46 @@ #endif } - else this->independent_objects=NULL; - } /*}}}*/ void IoModel::DeleteData(int num,...){/*{{{*/ + /*Intermediaries*/ va_list ap; int dataenum; - - /*Go through the entire list of enums and delete the corresponding data from the iomodel-data dataset: */ + vector::iterator iter; + + /*Go through the entire list of data and delete the corresponding data from the iomodel-data dataset: */ va_start(ap,num); for(int i=0;iindependents[dataenum]) xDelete(this->data[dataenum]); + + for(iter=data.begin();iterdata_enum==dataenum && !iodata->isindependent){ + delete *iter; + this->data.erase(iter); + break; + } + } } va_end(ap); } /*}}}*/ -void IoModel::DeleteData(IssmDouble* vector, int dataenum){/*{{{*/ - - /*do not erase independent variables for the AD mode computations!: */ - if(vector)if (!this->independents[dataenum]) xDelete(vector); - +void IoModel::DeleteData(IssmDouble* vector_in, int dataenum){/*{{{*/ + + vector::iterator iter; + + /*do not do anything if pointer is NULL*/ + if(!vector_in) return; + + /*do not delete if this is an independent variable*/ + for(iter=data.begin();iterdata_enum==dataenum && iodata->isindependent){ + return; + } + } + + /*Go ahead and delete*/ + xDelete(vector_in); } /*}}}*/ void IoModel::DeleteData(char*** pstringarray, int numstrings, int dataenum){/*{{{*/ @@ -429,6 +765,4 @@ void IoModel::FetchConstants(void){/*{{{*/ - int my_rank; - /*record descriptions; */ int record_enum; @@ -437,17 +771,16 @@ /*records: */ - int booleanint=0; - int integer=0; - IssmPDouble pscalar=0; - IssmDouble scalar=0; - char* string=NULL; - int string_size; + int booleanint = 0; + int integer = 0; + IssmPDouble pscalar = 0; + IssmDouble scalar = 0; + char *string = NULL; + int string_size; /*recover my_rank:*/ - my_rank=IssmComm::GetRank(); + int my_rank=IssmComm::GetRank(); /*Check that some fields have been allocated*/ _assert_(this->fid || my_rank); - _assert_(this->constants); /*Go find in the binary file, the position of the data we want to fetch: */ @@ -487,5 +820,5 @@ /*create BoolParam: */ if(record_enum!=MaximumNumberOfDefinitionsEnum && record_enum!=MaximumNumberOfDefinitionsEnum+1) - this->constants->AddObject(new BoolParam(record_enum,(bool)booleanint)); //cast to boolean + this->AddConstant(new IoConstant((bool)booleanint,record_enum)); //cast to boolean break; @@ -496,23 +829,28 @@ /*create IntParam: */ - this->constants->AddObject(new IntParam(record_enum,integer)); + this->AddConstant(new IoConstant(integer,record_enum)); break; case 3: - /*Read the scalar and broadcast it to other cpus. However, if this record has already - * been read in DeclareIndependents, we should not re-read it, but grab it from the iomodel->data - * slots: */ - if(this->independents[record_enum]){ - scalar=*(this->data[record_enum]); - } - else{ - if(fread(&pscalar,sizeof(IssmPDouble),1,this->fid)!=1) _error_("could not read scalar "); - ISSM_MPI_Bcast(&pscalar,1,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); - scalar=pscalar; - } - - /*create DoubleParam: */ - this->constants->AddObject(new DoubleParam(record_enum,scalar)); - + { + /*IssmDouble, check whether it is already there (from "declare independents")*/ + bool exists = false; + vector::iterator iter; + for(iter=constants.begin();iterdata_enum==record_enum){ + exists = true; + break; + } + } + if(!exists){ + if(fread(&pscalar,sizeof(IssmPDouble),1,this->fid)!=1) _error_("could not read scalar "); + ISSM_MPI_Bcast(&pscalar,1,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); + scalar=pscalar; + + /*create DoubleParam: */ + this->AddConstant(new IoConstant(scalar,record_enum)); + } + } break; case 4: @@ -535,36 +873,13 @@ /*Add string to parameters: */ - this->constants->AddObject(new StringParam(record_enum,string)); + this->AddConstant(new IoConstant(string,record_enum)); /*Free string*/ xDelete(string); - break; case 5: - /*We are not interested in this record, too memory intensive. Skip it: */ - /*skip: */ - fseek(fid,-sizeof(int),SEEK_CUR); //backtrak 1 integer - fseek(fid,record_length,SEEK_CUR); - break; case 6: - /*We are not interested in this record, too memory intensive. Skip it: */ - /*skip: */ - fseek(fid,-sizeof(int),SEEK_CUR); //backtrak 1 integer - fseek(fid,record_length,SEEK_CUR); - break; case 7: - /*We are not interested in this record, too memory intensive. Skip it: */ - /*skip: */ - fseek(fid,-sizeof(int),SEEK_CUR); //backtrak 1 integer - fseek(fid,record_length,SEEK_CUR); - break; - case 8: - /*We are not interested in this record, too memory intensive. Skip it: */ - /*skip: */ - fseek(fid,-sizeof(int),SEEK_CUR); //backtrak 1 integer - fseek(fid,record_length,SEEK_CUR); - break; - case 9: /*We are not interested in this record, too memory intensive. Skip it: */ @@ -573,8 +888,7 @@ fseek(fid,record_length,SEEK_CUR); break; - default: _error_("unknown record type:" << record_code); - break;; + break; } } @@ -597,5 +911,5 @@ /*create BoolParam: */ if(record_enum!=MaximumNumberOfDefinitionsEnum && record_enum!=MaximumNumberOfDefinitionsEnum+1) - this->constants->AddObject(new BoolParam(record_enum,(bool)booleanint)); //cast to a boolean + this->AddConstant(new IoConstant((bool)booleanint,record_enum)); //cast to a boolean break; @@ -605,14 +919,27 @@ /*create IntParam: */ - this->constants->AddObject(new IntParam(record_enum,integer)); + this->AddConstant(new IoConstant(integer,record_enum)); break; case 3: /*scalar. get it from cpu 0 */ - ISSM_MPI_Bcast(&pscalar,1,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); - scalar=pscalar; - /*create DoubleParam: */ - this->constants->AddObject(new DoubleParam(record_enum,scalar)); - + { + /*IssmDouble, check whether it is already there (from "declare independents")*/ + bool exists = false; + vector::iterator iter; + for(iter=constants.begin();iterdata_enum==record_enum){ + exists = true; + break; + } + } + if(!exists){ + ISSM_MPI_Bcast(&pscalar,1,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); + scalar=pscalar; + /*create DoubleParam: */ + this->AddConstant(new IoConstant(scalar,record_enum)); + } + } break; case 4: @@ -630,5 +957,5 @@ } /*Add string to parameters: */ - this->constants->AddObject(new StringParam(record_enum,string)); + this->AddConstant(new IoConstant(string,record_enum)); /*Free string*/ @@ -644,5 +971,5 @@ default: _error_("unknown record type:" << record_code); - break;; + break; } @@ -653,6 +980,4 @@ /*}}}*/ void IoModel::FetchData(bool* pboolean,int data_enum){/*{{{*/ - - int my_rank; /*output: */ @@ -661,5 +986,5 @@ /*recover my_rank:*/ - my_rank=IssmComm::GetRank(); + int my_rank=IssmComm::GetRank(); /*Set file pointer to beginning of the data: */ @@ -682,6 +1007,4 @@ void IoModel::FetchData(int* pinteger,int data_enum){/*{{{*/ - int my_rank; - /*output: */ int integer; @@ -689,5 +1012,5 @@ /*recover my_rank:*/ - my_rank=IssmComm::GetRank(); + int my_rank=IssmComm::GetRank(); /*Set file pointer to beginning of the data: */ @@ -709,6 +1032,4 @@ void IoModel::FetchData(IssmDouble* pscalar,int data_enum){/*{{{*/ - int my_rank; - /*output: */ IssmPDouble scalar; @@ -716,5 +1037,5 @@ /*recover my_rank:*/ - my_rank=IssmComm::GetRank(); + int my_rank=IssmComm::GetRank(); /*Set file pointer to beginning of the data: */ @@ -735,6 +1056,4 @@ /*}}}*/ void IoModel::FetchData(char** pstring,int data_enum){/*{{{*/ - - int my_rank; /*output: */ @@ -744,5 +1063,5 @@ /*recover my_rank:*/ - my_rank=IssmComm::GetRank(); + int my_rank=IssmComm::GetRank(); /*Set file pointer to beginning of the data: */ @@ -781,6 +1100,4 @@ /*}}}*/ void IoModel::FetchData(int** pmatrix,int* pM,int* pN,int data_enum){/*{{{*/ - - int my_rank; int i,j; @@ -792,5 +1109,5 @@ /*recover my_rank:*/ - my_rank=IssmComm::GetRank(); + int my_rank=IssmComm::GetRank(); /*Set file pointer to beginning of the data: */ @@ -850,13 +1167,23 @@ void IoModel::FetchData(IssmDouble** pmatrix,int* pM,int* pN,int data_enum){/*{{{*/ - int my_rank; - + /*First, look if has already been loaded (might be an independent variable)*/ + vector::iterator iter; + for(iter=data.begin();iterdata_enum==data_enum){ + *pmatrix=iodata->data; + if(pM) *pM=iodata->M; + if(pN) *pN=iodata->N; + return; + } + } + /*output: */ - int M,N; - IssmPDouble* matrix=NULL; - int code=0; + int M,N; + IssmPDouble *matrix = NULL; + int code = 0; /*recover my_rank:*/ - my_rank=IssmComm::GetRank(); + int my_rank=IssmComm::GetRank(); /*Set file pointer to beginning of the data: */ @@ -868,5 +1195,5 @@ /*We have to read a matrix from disk. First read the dimensions of the matrix, then the whole matrix: */ /*numberofelements: */ - if(my_rank==0){ + if(my_rank==0){ if(fread(&M,sizeof(int),1,fid)!=1) _error_("could not read number of rows for matrix "); } @@ -888,14 +1215,6 @@ ISSM_MPI_Bcast(matrix,M*N,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); - _assert_(this->independents); - if(this->independents[data_enum]){ - /*this data has already been checked out! So cancel all that we've done here, and return - * the data[data_enum] directly: */ - *pmatrix=this->data[data_enum]; - } - else{ - *pmatrix=xNew(M*N); - for(int i=0;i(M*N); + for(int i=0;i(matrix); } @@ -903,11 +1222,9 @@ *pmatrix=NULL; /*Assign output pointers: */ - if (pM)*pM=M; - if (pN)*pN=N; + if(pM) *pM=M; + if(pN) *pN=N; } /*}}}*/ void IoModel::FetchData(char*** pstrings,int* pnumstrings,int data_enum){/*{{{*/ - - int my_rank; int i; @@ -923,5 +1240,5 @@ /*recover my_rank:*/ - my_rank=IssmComm::GetRank(); + int my_rank=IssmComm::GetRank(); /*Set file pointer to beginning of the data: */ @@ -975,7 +1292,4 @@ int i; - - int my_rank; - /*output: */ IssmDouble** matrices=NULL; @@ -990,5 +1304,5 @@ /*recover my_rank:*/ - my_rank=IssmComm::GetRank(); + int my_rank=IssmComm::GetRank(); /*Set file pointer to beginning of the data: */ @@ -1105,8 +1419,9 @@ va_list ap; - int dataenum; + int dataenum,code,layout; IssmDouble *matrix = NULL; int M,N; - int i; + bool exists; + vector::iterator iter; /*Go through the entire list of enums and fetch the corresponding data. Add it to the iomodel->data dataset. Everything @@ -1114,31 +1429,31 @@ va_start(ap,num); - for(i=0; i=0); - - if (this->independents[dataenum]){ + exists = false; + + for(iter=data.begin();iterdata_enum==dataenum){ + /*Already there, no need to fetch it*/ + _assert_(iodata->isindependent); + exists = true; + break; + } + } + + if(exists){ /*this data has already been checked out! Continue: */ continue; } - - /*Some checks in debugging mode*/ - /*{{{*/ - #ifdef _ISSM_DEBUG_ - _assert_(dataenumdata[dataenum]){ - _error_("Info: trying to fetch " << EnumToStringx(dataenum) << " but previous pointer has not been freed (DeleteData has not been called)"); - } - #endif - /*}}}*/ - - /*Add to this->data: */ - this->FetchData(&matrix,&M,&N,dataenum); - this->data[dataenum]=matrix; + else{ + /*Add to this->data: */ + this->SetFilePointerToData(&code,&layout,dataenum); + this->FetchData(&matrix,&M,&N,dataenum); + this->AddData(new IoData(matrix,code,layout,M,N,dataenum)); + } } va_end(ap); - } /*}}}*/ @@ -1373,14 +1688,13 @@ ISSM_MPI_Bcast(pmatrix,M*N,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); - _assert_(this->independents); - if(this->independents[data_enum]){ - /*this data has already been checked out! So cancel all that we've done here, and return - * the data[data_enum] directly: */ - matrix=this->data[data_enum]; - } - else{ + //if(this->independents[data_enum]){ FIXME + // /*this data has already been checked out! So cancel all that we've done here, and return + // * the data[data_enum] directly: */ + // matrix=this->data[data_enum]; + //} + //else{ matrix=xNew(M*N); for (int i=0;i(pmatrix); } @@ -1479,14 +1793,13 @@ ISSM_MPI_Bcast(pmatrix,M*N,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); - _assert_(this->independents); - if(this->independents[data_enum]){ - /*this data has already been checked out! So cancel all that we've done here, and return - * the data[data_enum] directly: */ - matrix=this->data[data_enum]; - for (int i=0;i(matrix[i]); - } - else{ + //if(this->independents[data_enum]){ FIXME + // /*this data has already been checked out! So cancel all that we've done here, and return + // * the data[data_enum] directly: */ + // matrix=this->data[data_enum]; + // for (int i=0;i(matrix[i]); + //} + //else{ for (int i=0;i(pmatrix); } @@ -1525,4 +1838,18 @@ void IoModel::FetchDataToInput(Elements* elements,int vector_enum,IssmDouble default_value){/*{{{*/ + /*First, look whether it is not already loaded in this->data*/ + vector::iterator iter; + for(iter=data.begin();iterdata_enum==vector_enum){ + _assert_(iodata->code==7); + for(int i=0;iSize();i++){ + Element* element=xDynamicCast(elements->GetObjectByOffset(i)); + element->InputCreate(iodata->data,this,iodata->M,iodata->N,iodata->layout,vector_enum,iodata->code);//we need i to index into elements. + } + return; + } + } + /*intermediary: */ int code,vector_layout; @@ -1544,9 +1871,22 @@ } - /*Free ressources. Pay attention to not freeing an AD mode independent variable though!:*/ - if(!this->independents[vector_enum]) xDelete(doublearray); + /*Free ressources*/ + xDelete(doublearray); } /*}}}*/ void IoModel::FetchDataToInput(Elements* elements,int vector_enum){/*{{{*/ + + /*First, look whether it is not already loaded in this->data*/ + vector::iterator iter; + for(iter=data.begin();iterdata_enum==vector_enum){ + for(int i=0;iSize();i++){ + Element* element=xDynamicCast(elements->GetObjectByOffset(i)); + element->InputCreate(iodata->data,this,iodata->M,iodata->N,iodata->layout,vector_enum,iodata->code);//we need i to index into elements. + } + return; + } + } /*intermediary: */ @@ -1615,6 +1955,6 @@ break; } - /*Free ressources. Pay attention to not freeing an AD mode independent variable though!:*/ - if (!this->independents[vector_enum]) xDelete(doublearray); + /*Free ressources*/ + xDelete(doublearray); xDelete(string); } Index: /issm/trunk-jpl/src/c/classes/IoModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/IoModel.h (revision 20652) +++ /issm/trunk-jpl/src/c/classes/IoModel.h (revision 20653) @@ -9,4 +9,5 @@ #include "../shared/Enum/Enum.h" +#include class Parameters; @@ -15,12 +16,45 @@ class Option; +class IoConstant { /*holds single IssmDouble, int, bool and char from input*/ + public: + bool isindependent; + int data_enum; + char* name; + Param* constant; + + ~IoConstant(); + IoConstant(); + IoConstant(bool value,int enum_in); + IoConstant(int value,int enum_in); + IoConstant(IssmDouble value,int enum_in); + IoConstant(char* value,int enum_in); +}; + +class IoData { /*holds temporary data (array), memory intensive*/ + public: + bool isindependent; + int data_enum; + char* name; + int M,N; + int code; + int layout; + IssmDouble* data; + + ~IoData(); + IoData(); + IoData(IssmDouble* matrix,int code,int layout_in,int M,int N,int enum_in); +}; + class IoModel { private: - Parameters *constants; //this dataset holds all IssmDouble, int, bool and char from input + std::vector constants; //this dataset holds all IssmDouble, int, bool and char from input + std::vector data; //this dataset holds temporary data, memory intensive + + /*for AD mode: to keep track of our independent variables we fetch:*/ + //bool *independents; + //DataSet *independent_objects; public: - IssmDouble **data; //this dataset holds temporary data, memory intensive. - /*pointer to input file*/ FILE *fid; @@ -55,8 +89,4 @@ int constraintcounter; //keep track of how many constraints are being created in each analysis - /*for AD mode: to keep track of our independent variables we fetch:*/ - bool *independents; - DataSet *independent_objects; - /*Methods*/ ~IoModel(); @@ -64,10 +94,20 @@ IoModel(FILE* iomodel_handle,int solution_enum_in,bool trace,IssmPDouble* X); + /*NEW*/ + void AddConstant(IoConstant* constant_in); + void AddConstantIndependent(IoConstant* constant_in); + void AddData(IoData* data_in); + void AddDataIndependent(IoData* data_in); + void FindConstant(bool* pvalue,int constant_enum); + void FindConstant(int* pvalue,int constant_enum); + void FindConstant(IssmDouble* pvalue,int constant_enum); + void FindConstant(char **pvalue,int constant_enum); + void FetchIndependentConstant(int* pXcount,IssmPDouble* X,int name); + void FetchIndependentData(int* pXcount,IssmPDouble* X,int name); + int NumIndependents(); + void FillIndependents(IssmDouble* xp); + /*Input/Output*/ void CheckEnumSync(void); - void FindConstant(bool *poutput,int constant_enum); - void FindConstant(int *poutput,int constant_enum); - void FindConstant(IssmDouble *poutput,int constant_enum); - void FindConstant(char **poutput,int constant_enum); Param *CopyConstantObject(int constant_enum); IssmDouble *Data(int dataenum); Index: /issm/trunk-jpl/src/c/main/kriging.cpp =================================================================== --- /issm/trunk-jpl/src/c/main/kriging.cpp (revision 20652) +++ /issm/trunk-jpl/src/c/main/kriging.cpp (revision 20653) @@ -145,5 +145,4 @@ iomodel->fid=fid; iomodel->CheckEnumSync(); - iomodel->independents=xNew(MaximumNumberOfDefinitionsEnum); for(int i=0;iindependents[i]=false; iomodel->FetchData(&x,&M,&N,0); nobs=M*N; iomodel->FetchData(&y,&M,&N,1); _assert_(M*N==nobs); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp (revision 20652) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp (revision 20653) @@ -133,24 +133,13 @@ /*Independents have already been recovered in iomodel->DeclareIndependents. Just do some more processing. *In particular, figure out num_independents, and create the state vector xp, or size num_independents x 1 :*/ - num_ind=0; - for(i=0;iindependent_objects->Size();i++){ - IndependentObject* ind=(IndependentObject*)iomodel->independent_objects->GetObjectByOffset(i); - num_ind+=ind->NumIndependents(); - } + num_ind=iomodel->NumIndependents(); + parameters->AddObject(new IntParam(AutodiffNumIndependentsEnum,num_ind)); + if(num_ind){ xp=xNew(num_ind); - xp_backup=xp; - for(i=0;iindependent_objects->Size();i++){ - IndependentObject* ind=(IndependentObject*)iomodel->independent_objects->GetObjectByOffset(i); - ind->FillIndependents(iomodel->data,xp); - local_num_ind=ind->NumIndependents(); xp=xp+local_num_ind; - } - xp=xp_backup; parameters->AddObject(new DoubleVecParam(AutodiffXpEnum,xp,num_ind)); + iomodel->FillIndependents(xp); + parameters->AddObject(new DoubleVecParam(AutodiffXpEnum,xp,num_ind)); xDelete(xp); } - parameters->AddObject(new IntParam(AutodiffNumIndependentsEnum,num_ind)); - - /*Don't forget to copy iomodel->independent_objects to parameters: */ - parameters->AddObject(new DataSetParam(AutodiffIndependentObjectsEnum,iomodel->independent_objects)); /*}}}*/ }