Changeset 26117

Timestamp:

03/18/21 20:58:51 (4 years ago)

Author:

Mathieu Morlighem

Message:

BUG: fixing problem with matrix

Location:

issm/trunk-jpl/src/c

Files:

• Makefile.am (modified) (1 diff)
• toolkits/objects/Matrix.h (modified) (3 diffs)
• toolkits/petsc/objects/PetscMat.cpp (modified) (1 diff)
• toolkits/petsc/objects/PetscMat.h (modified) (1 diff)
• toolkits/petsc/patches/MatToMPISerial.cpp (modified) (1 diff)
• toolkits/petsc/patches/MatToSerial.cpp (deleted)
• toolkits/petsc/patches/VecToMPISerial.cpp (modified) (1 diff)
• toolkits/petsc/patches/petscpatches.h (modified) (1 diff)

issm/trunk-jpl/src/c/Makefile.am

@@ -373 +373 @@
         ./toolkits/petsc/patches \
         ./toolkits/petsc/patches/VecToMPISerial.cpp \
-        ./toolkits/petsc/patches/MatToSerial.cpp \
+        ./toolkits/petsc/patches/MatToMPISerial.cpp \
         ./toolkits/petsc/patches/NewVec.cpp \
         ./toolkits/petsc/patches/PetscOptionsDetermineSolverType.cpp \

issm/trunk-jpl/src/c/toolkits/objects/Matrix.h

@@ -266 +266 @@
                 }
                 /*}}}*/
-                doubletype* ToSerial(void){/*{{{*/
+                doubletype* ToMPISerial0(void){/*{{{*/
 
                         doubletype* output=NULL;
@@ -272 +272 @@
                         if(type==PetscMatType){
                                 #ifdef _HAVE_PETSC_
-                                output=this->pmatrix->ToSerial();
-                                #endif
-                        }
-                        else{
-                                output=this->imatrix->ToSerial();
+                                output=this->pmatrix->ToMPISerial0();
+                                #endif
+                        }
+                        else{
+                                output=this->imatrix->ToMPISerial0();
                         }
 
@@ -288 +288 @@
                         if(type==PetscMatType){
                                 #ifdef _HAVE_PETSC_
-                                output=this->pmatrix->ToSerial();
+                                output=this->pmatrix->ToMPISerial();
                                 #endif
                         }

issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscMat.cpp

@@ -171 +171 @@
 }
 /*}}}*/
-IssmDouble* PetscMat::ToSerial(void){/*{{{*/
+IssmDouble* PetscMat::ToMPISerial(void){/*{{{*/
 
          IssmDouble* output=NULL;
-
-         MatToSerial(&output,this->matrix,IssmComm::GetComm());
+         MatToMPISerial(&output,this->matrix,IssmComm::GetComm(),true);
          return output;
 
 }
 /*}}}*/
-IssmDouble* PetscMat::ToMPISerial(void){/*{{{*/
+IssmDouble* PetscMat::ToMPISerial0(void){/*{{{*/
 
          IssmDouble* output=NULL;
-         
-         MatToMPISerial(&output,this->matrix,IssmComm::GetComm());
+         MatToMPISerial(&output,this->matrix,IssmComm::GetComm(),false);
          return output;
 

issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscMat.h

@@ -49 +49 @@
                 void MatMult(PetscVec* X,PetscVec* AX);
                 PetscMat* Duplicate(void);
-                IssmDouble* ToSerial(void);
                 IssmDouble* ToMPISerial(void);
+                IssmDouble* ToMPISerial0(void);
                 void SetValues(int m,int* idxm,int n,int* idxn,IssmDouble* values,InsMode mode);
                 void Convert(MatrixType type);

issm/trunk-jpl/src/c/toolkits/petsc/patches/MatToMPISerial.cpp

@@ -12 +12 @@
 #include "../../../shared/shared.h"
 
-void MatToMPISerial(double** poutmatrix,Mat matrix,ISSM_MPI_Comm comm){
+void MatToMPISerial(double** poutmatrix,Mat matrix,ISSM_MPI_Comm comm,bool broadcast){
 
+        int i;
+        int my_rank;
+        int num_procs;
+
+        /*Petsc variables*/
+        PetscInt lower_row,upper_row; 
+        int range;
+        int M,N; //size of matrix
+        ISSM_MPI_Status status;
+        int* idxm=NULL;
+        int* idxn=NULL; 
+        double* local_matrix=NULL; /*matrix local to each node used for temporary holding matrix values*/
+        int buffer[3];
+
+        /*recover my_rank and num_procs:*/
+        ISSM_MPI_Comm_size(comm,&num_procs);
+        ISSM_MPI_Comm_rank(comm,&my_rank);
+
+        /*Output*/
         double* outmatrix=NULL;
-
-        /*Call MatToSerial:*/
-        MatToSerial(&outmatrix,matrix,comm);
 
         /*get matrix size: */
         MatGetSize(matrix,&M,&N);
 
-        /*Allocate on cpus:*/
-        if (my_rank!=0)outmatrix=xNew<double>(M*N);
+        /*partition: */
+        MatGetOwnershipRange(matrix,&lower_row,&upper_row);    
+        upper_row--; 
+        range=upper_row-lower_row+1;
 
-        /*Broadcast:*/
-        ISSM_MPI_Bcast(outmatrix,M*N,ISSM_MPI_PDOUBLE,0,comm);
+        /*Local and global allocation*/
+        if(broadcast || my_rank==0){ 
+                outmatrix=xNew<double>(M*N);
+        }
 
-        /*Assign output pointer:*/
+        if (range){
+                local_matrix=xNew<double>(N*range);
+                idxm=xNew<int>(range);  
+                idxn=xNew<int>(N);  
+
+                for (i=0;i<N;i++){
+                        *(idxn+i)=i;
+                }
+                for (i=0;i<range;i++){
+                        *(idxm+i)=lower_row+i;
+                }
+
+                MatGetValues(matrix,range,idxm,N,idxn,local_matrix);     
+        }
+
+        /*Now each node holds its local_matrix containing range rows. 
+         * We send these rows to the matrix on node 0*/
+
+        for (i=1;i<num_procs;i++){
+                if (my_rank==i){ 
+                        buffer[0]=my_rank;
+                        buffer[1]=lower_row;
+                        buffer[2]=range;
+                        ISSM_MPI_Send(buffer,3,ISSM_MPI_INT,0,1,comm);   
+                        if (range)ISSM_MPI_Send(local_matrix,N*range,ISSM_MPI_PDOUBLE,0,1,comm); 
+                }
+                if (my_rank==0){
+                        ISSM_MPI_Recv(buffer,3,ISSM_MPI_INT,i,1,comm,&status); 
+                        if (buffer[2])ISSM_MPI_Recv(outmatrix+(buffer[1]*N),N*buffer[2],ISSM_MPI_PDOUBLE,i,1,comm,&status);
+                }
+        } 
+        if (my_rank==0){ 
+                //Still have the local_matrix on node 0 to take care of.
+                if (range) memcpy(outmatrix,local_matrix,N*range*sizeof(double));
+        } 
+
+        if(broadcast){
+                /*Broadcast:*/
+                ISSM_MPI_Bcast(outmatrix,M*N,ISSM_MPI_PDOUBLE,0,comm);
+        }
+
+        /*Assign output pointer: */
+        xDelete<int>(idxm);
+        xDelete<int>(idxn);
+        xDelete<double>(local_matrix);
         *poutmatrix=outmatrix;
 }

issm/trunk-jpl/src/c/toolkits/petsc/patches/VecToMPISerial.cpp

@@ -82 +82 @@
         if (my_rank==0){ 
                 //Still have the local_vector on node 0 to take care of.
-                if (range)memcpy(gathered_vector+lower_row,local_vector,range*sizeof(double));
+                if (range) memcpy(gathered_vector+lower_row,local_vector,range*sizeof(double));
         }
 

issm/trunk-jpl/src/c/toolkits/petsc/patches/petscpatches.h

@@ -27 +27 @@
 void VecFree(Vec* pvec);
 void KSPFree(KSP* pksp);
-int MatPartition(Mat* poutmatrix,Mat matrixA,double* row_partition_vector,int row_partition_vector_size ,
-                double* col_partition_vector,int col_partition_vector_size);
+int MatPartition(Mat* poutmatrix,Mat matrixA,double* row_partition_vector,int row_partition_vector_size, double* col_partition_vector,int col_partition_vector_size);
 void PetscOptionsDetermineSolverType(int* psolver_type);
 void MatMultPatch(Mat A,Vec X, Vec AX,ISSM_MPI_Comm comm);
-void MatToSerial(double** poutmatrix,Mat matrix,ISSM_MPI_Comm comm);
-void MatToMPISerial(double** poutmatrix,Mat matrix,ISSM_MPI_Comm comm);
+void MatToMPISerial(double** poutmatrix,Mat matrix,ISSM_MPI_Comm comm,bool broadcast=true);
 Vec  SerialToVec(double* vector,int vector_size);
 InsertMode ISSMToPetscInsertMode(InsMode mode);