Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 23601) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 23602) @@ -239,5 +239,4 @@ ./modules/ResetFSBasalBoundaryConditionx/ResetFSBasalBoundaryConditionx.cpp\ ./modules/Solverx/Solverx.cpp\ - ./modules/VecMergex/VecMergex.cpp\ ./modules/Mergesolutionfromftogx/Mergesolutionfromftogx.cpp\ ./cores/ProcessArguments.cpp\ @@ -352,5 +351,4 @@ ./toolkits/petsc/patches/VecToMPISerial.cpp\ ./toolkits/petsc/patches/MatToSerial.cpp\ - ./toolkits/petsc/patches/VecMerge.cpp\ ./toolkits/petsc/patches/NewVec.cpp\ ./toolkits/petsc/patches/PetscOptionsDetermineSolverType.cpp\ Index: /issm/trunk-jpl/src/c/classes/Node.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Node.cpp (revision 23601) +++ /issm/trunk-jpl/src/c/classes/Node.cpp (revision 23602) @@ -590,52 +590,59 @@ } /*}}}*/ -void Node::VecMerge(Vector* ug, IssmDouble* vector_serial,int setenum){/*{{{*/ - - IssmDouble *values = NULL; - int *indices = NULL; - int count = 0; - int i; - - if(setenum==FsetEnum){ - if(this->indexing.fsize){ - indices=xNew(this->indexing.fsize); - values=xNew(this->indexing.fsize); - - for(i=0;iindexing.gsize;i++){ - if(this->indexing.f_set[i]){ - _assert_(vector_serial); - values[count]=vector_serial[this->indexing.fdoflist[count]]; - indices[count]=this->indexing.gdoflist[i]; - count++; - } - } - - /*Add values into ug: */ - ug->SetValues(this->indexing.fsize,indices,values,INS_VAL); - } - } - else if(setenum==SsetEnum){ - if(this->indexing.ssize){ - indices=xNew(this->indexing.ssize); - values=xNew(this->indexing.ssize); - - for(i=0;iindexing.gsize;i++){ - if(this->indexing.s_set[i]){ - _assert_(vector_serial); - values[count]=vector_serial[this->indexing.sdoflist[count]]; - indices[count]=this->indexing.gdoflist[i]; - count++; - } - } - - /*Add values into ug: */ - ug->SetValues(this->indexing.ssize,indices,values,INS_VAL); - } - } - else _error_("VecMerge can only merge from the s or f-set onto the g-set!"); - - /*Free ressources:*/ - xDelete(values); - xDelete(indices); +void Node::VecMerge(Vector* ug,IssmDouble* local_uf,int* indices_uf,int* pindex_uf,IssmDouble* local_ys,int* indices_ys,int* pindex_ys){/*{{{*/ + + /*Only perform operation if not clone*/ + if(this->IsClone()) return; + + /*Recover indices*/ + int ind_uf = *pindex_uf; + int ind_ys = *pindex_ys; + + if(this->indexing.fsize){ + int* indices = xNew(this->indexing.fsize); + IssmDouble* values = xNew(this->indexing.fsize); + + int count = 0; + for(int i=0;iindexing.gsize;i++){ + if(this->indexing.f_set[i]){ + _assert_(local_uf); + _assert_(this->indexing.fdoflist[count]==indices_uf[ind_uf]); + + values[count]=local_uf[ind_uf]; + indices[count]=this->indexing.gdoflist[i]; + count++; + ind_uf++; + } + } + ug->SetValues(this->indexing.fsize,indices,values,INS_VAL); + /*Free ressources:*/ + xDelete(values); + xDelete(indices); + } + if(this->indexing.ssize){ + int* indices = xNew(this->indexing.ssize); + IssmDouble* values = xNew(this->indexing.ssize); + + int count = 0; + for(int i=0;iindexing.gsize;i++){ + if(this->indexing.s_set[i]){ + _assert_(local_ys); + _assert_(this->indexing.sdoflist[count]==indices_ys[ind_ys]); + + values[count]=local_ys[ind_ys]; + indices[count]=this->indexing.gdoflist[i]; + count++; + ind_ys++; + } + } + ug->SetValues(this->indexing.ssize,indices,values,INS_VAL); + /*Free ressources:*/ + xDelete(values); + xDelete(indices); + } + + /*Update index values*/ + *pindex_uf = ind_uf; + *pindex_ys = ind_ys; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Node.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Node.h (revision 23601) +++ /issm/trunk-jpl/src/c/classes/Node.h (revision 23602) @@ -81,5 +81,5 @@ int Sid(void); void UpdateCloneDofs(int* alltruerows,int setenum); - void VecMerge(Vector* ug, IssmDouble* vector_serial,int setenum); + void VecMerge(Vector* ug,IssmDouble* local_uf,int* indices_uf,int* pindex_uf,IssmDouble* local_ys,int* indices_ys,int* pindex_ys); void VecReduce(Vector* vector, IssmDouble* ug_serial,int setnum); void SetApproximation(int in_approximation); Index: /issm/trunk-jpl/src/c/modules/Mergesolutionfromftogx/Mergesolutionfromftogx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Mergesolutionfromftogx/Mergesolutionfromftogx.cpp (revision 23601) +++ /issm/trunk-jpl/src/c/modules/Mergesolutionfromftogx/Mergesolutionfromftogx.cpp (revision 23602) @@ -3,5 +3,4 @@ */ -#include "../VecMergex/VecMergex.h" #include "../../shared/io/io.h" #include "./Mergesolutionfromftogx.h" @@ -9,17 +8,11 @@ void Mergesolutionfromftogx( Vector** pug, Vector* uf, Vector* ys, Nodes* nodes, Parameters* parameters, bool flag_ys0){ - /*output: */ - Vector* ug=NULL; - - /*intermediary: */ - int gsize,fsize,ssize; - /*Display message*/ if(VerboseModule()) _printf0_(" Merging solution vector from fset to gset\n"); /*first, get gsize, fsize and ssize: */ - gsize=nodes->NumberOfDofs(GsetEnum); - fsize=nodes->NumberOfDofs(FsetEnum); - ssize=nodes->NumberOfDofs(SsetEnum); + int gsize=nodes->NumberOfDofs(GsetEnum); + int fsize=nodes->NumberOfDofs(FsetEnum); + int ssize=nodes->NumberOfDofs(SsetEnum); /*serialize uf and ys: those two vectors will be indexed by the nodes, who are the only ones @@ -31,18 +24,31 @@ } + /*Get local vectors ys and uf*/ + int *indices_ys = NULL; + IssmDouble *local_ys = NULL; + ys->GetLocalVector(&local_ys,&indices_ys); + int *indices_uf = NULL; + IssmDouble *local_uf = NULL; + uf->GetLocalVector(&local_uf,&indices_uf); + /*initialize ug: */ - ug=new Vector(gsize); + Vector* ug=new Vector(gsize); - /*Merge f set back into g set: */ - if(fsize){ - VecMergex(ug,uf,nodes,parameters,FsetEnum); + /*Let nodes figure it out*/ + int index_uf = 0; + int index_ys = 0; + for(int i=0;iSize();i++){ + Node* node=(Node*)nodes->GetObjectByOffset(i); + node->VecMerge(ug,local_uf,indices_uf,&index_uf,local_ys,indices_ys,&index_ys); } - /*Merge s set back into g set: */ - if(ssize){ - VecMergex(ug,ys,nodes,parameters,SsetEnum); - } + /*Assemble vector: */ + ug->Assemble(); - /*Assign correct pointer*/ + /*Cleanup and assounf output pointer*/ + xDelete(indices_uf); + xDelete(indices_ys); + xDelete(local_uf); + xDelete(local_ys); *pug=ug; } Index: /issm/trunk-jpl/src/c/modules/modules.h =================================================================== --- /issm/trunk-jpl/src/c/modules/modules.h (revision 23601) +++ /issm/trunk-jpl/src/c/modules/modules.h (revision 23602) @@ -99,4 +99,3 @@ #include "./UpdateDynamicConstraintsx/UpdateDynamicConstraintsx.h" #include "./VertexCoordinatesx/VertexCoordinatesx.h" -#include "./VecMergex/VecMergex.h" #endif Index: /issm/trunk-jpl/src/c/toolkits/issm/IssmAbsVec.h =================================================================== --- /issm/trunk-jpl/src/c/toolkits/issm/IssmAbsVec.h (revision 23601) +++ /issm/trunk-jpl/src/c/toolkits/issm/IssmAbsVec.h (revision 23602) @@ -37,4 +37,5 @@ virtual void GetSize(int* pM)=0; virtual void GetLocalSize(int* pM)=0; + virtual void GetLocalVector(double** pvector,int** pindices)=0; virtual IssmAbsVec* Duplicate(void)=0; virtual void Set(doubletype value)=0; Index: /issm/trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h =================================================================== --- /issm/trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h (revision 23601) +++ /issm/trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h (revision 23602) @@ -337,4 +337,41 @@ } /*}}}*/ + void GetLocalVector(double** pvector,int** pindices){/*{{{*/ + + _assert_(this->vector); + + /*First, check that vector size is not 0*/ + int vector_size; + this->GetSize(&vector_size); + if(vector_size==0){ + *pvector=NULL; + *pindices=NULL; + return; + } + + /*Get Ownership range*/ + int lower_row,upper_row; + GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm()); + int range=upper_row-lower_row; + + /*return NULL if no range*/ + if(range==0){ + *pvector=NULL; + *pindices=NULL; + return; + } + + /*Build indices*/ + int* indices=xNew(range); + for(int i=0;im); + IssmDouble* values =xNew(range); + xMemCpy(values,this->vector,this->m); + + *pvector = values; + *pindices = indices; + } /*}}}*/ IssmMpiVec* Duplicate(void){/*{{{*/ Index: /issm/trunk-jpl/src/c/toolkits/issm/IssmSeqVec.h =================================================================== --- /issm/trunk-jpl/src/c/toolkits/issm/IssmSeqVec.h (revision 23601) +++ /issm/trunk-jpl/src/c/toolkits/issm/IssmSeqVec.h (revision 23602) @@ -141,9 +141,9 @@ /*}}}*/ void GetLocalSize(int* pM){/*{{{*/ - *pM=this->M; - - } - /*}}}*/ + }/*}}}*/ + void GetLocalVector(double** pvector,int** pindices){/*{{{*/ + _error_("not implemented"); + } /*}}}*/ IssmSeqVec* Duplicate(void){/*{{{*/ Index: /issm/trunk-jpl/src/c/toolkits/issm/IssmVec.h =================================================================== --- /issm/trunk-jpl/src/c/toolkits/issm/IssmVec.h (revision 23601) +++ /issm/trunk-jpl/src/c/toolkits/issm/IssmVec.h (revision 23602) @@ -153,4 +153,7 @@ } /*}}}*/ + void GetLocalVector(double** pvector,int** pindices){/*{{{*/ + vector->GetLocalVector(pvector,pindices); + } /*}}}*/ IssmVec* Duplicate(void){/*{{{*/ Index: /issm/trunk-jpl/src/c/toolkits/objects/Vector.h =================================================================== --- /issm/trunk-jpl/src/c/toolkits/objects/Vector.h (revision 23601) +++ /issm/trunk-jpl/src/c/toolkits/objects/Vector.h (revision 23602) @@ -217,4 +217,15 @@ } /*}}}*/ + void GetLocalVector(double** pvector,int** pindices){_assert_(this);/*{{{*/ + + if(type==PetscVecType){ + #ifdef _HAVE_PETSC_ + this->pvector->GetLocalVector(pvector,pindices); + #endif + } + else this->ivector->GetLocalVector(pvector,pindices); + + } + /*}}}*/ Vector* Duplicate(void){_assert_(this);/*{{{*/ Index: /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.cpp =================================================================== --- /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.cpp (revision 23601) +++ /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.cpp (revision 23602) @@ -113,14 +113,41 @@ _assert_(this->vector); VecGetSize(this->vector,pM); - -} -/*}}}*/ -void PetscVec::GetLocalSize(int* pM){/*{{{*/ - - _assert_(this->vector); - VecGetLocalSize(this->vector,pM); - -} -/*}}}*/ +} +/*}}}*/ +void PetscVec::GetLocalVector(double** pvector,int** pindices){/*{{{*/ + + _assert_(this->vector); + + /*First, check that vector size is not 0*/ + int vector_size; + this->GetSize(&vector_size); + if(vector_size==0){ + *pvector=NULL; + *pindices=NULL; + return; + } + + /*Get Ownership range*/ + PetscInt lower_row,upper_row; + VecGetOwnershipRange(this->vector,&lower_row,&upper_row); + int range=upper_row-lower_row; + + /*return NULL if no range*/ + if(range==0){ + *pvector=NULL; + *pindices=NULL; + return; + } + + /*Build indices*/ + int* indices=xNew(range); + for(int i=0;i(range); + VecGetValues(this->vector,range,indices,values); + + *pvector = values; + *pindices = indices; +} /*}}}*/ PetscVec* PetscVec::Duplicate(void){/*{{{*/ Index: /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.h =================================================================== --- /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.h (revision 23601) +++ /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.h (revision 23602) @@ -45,4 +45,5 @@ void GetSize(int* pM); void GetLocalSize(int* pM); + void GetLocalVector(double** pvector,int** pindices); PetscVec* Duplicate(void); void Set(IssmDouble value); Index: sm/trunk-jpl/src/c/toolkits/petsc/patches/VecMerge.cpp =================================================================== --- /issm/trunk-jpl/src/c/toolkits/petsc/patches/VecMerge.cpp (revision 23601) +++ (revision ) @@ -1,64 +1,0 @@ -/*!\file: VecMerge.cpp - * \brief merge vector B into A using partitioning vector A(row_partition_vector)=B; - */ - -#ifdef HAVE_CONFIG_H - #include -#else -#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" -#endif - -/*Petsc includes: */ -#include -#include -#include - -#include "./petscpatches.h" - -#include "../../../shared/shared.h" - -void VecMerge(Vec A, Vec B, double* row_partition_vector,int row_partition_size){ - - /*Petsc matrix*/ - int lower_row,upper_row,range; - int* idxm=NULL; - double* values=NULL; - - /*Vector sizes: */ - int MB; - - VecGetSize(B,&MB); - - /*If the dimension of the partitioning vector is not the same as that of vector B, we have a problem: */ - if ((row_partition_size !=MB) ){ - _error_("Dimensions of partitioning vector incompatible with dimensions of input vector\n"); - } - - /*Get values from vector B and plug them into vector A, using the partitioning vector*/ - VecGetOwnershipRange(B,&lower_row,&upper_row); - upper_row--; - range=upper_row-lower_row+1; - - if (range){ - /*This node owns rows of vector B, get them*/ - idxm=xNew(range); - values=xNew(range); - for(int i=0;i(idxm); - xDelete(values); -} Index: /issm/trunk-jpl/src/c/toolkits/petsc/patches/petscpatches.h =================================================================== --- /issm/trunk-jpl/src/c/toolkits/petsc/patches/petscpatches.h (revision 23601) +++ /issm/trunk-jpl/src/c/toolkits/petsc/patches/petscpatches.h (revision 23602) @@ -30,5 +30,4 @@ double* col_partition_vector,int col_partition_vector_size); void PetscOptionsDetermineSolverType(int* psolver_type); -void VecMerge(Vec A, Vec B, double* row_partition_vector,int row_partition_size); void MatMultPatch(Mat A,Vec X, Vec AX,ISSM_MPI_Comm comm); void MatToSerial(double** poutmatrix,Mat matrix,ISSM_MPI_Comm comm);