/[escript]/trunk/dudley/src/Mesh_optimizeDOFDistribution.c
ViewVC logotype

Diff of /trunk/dudley/src/Mesh_optimizeDOFDistribution.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3126 by jfenwick, Wed Sep 1 00:37:53 2010 UTC revision 3224 by jfenwick, Wed Sep 29 05:19:37 2010 UTC
# Line 11  Line 11 
11  *  *
12  *******************************************************/  *******************************************************/
13    
   
14  /**************************************************************/  /**************************************************************/
15    
16  /*   Dudley: Mesh: optimizes the distribution of DOFs across processors */  /*   Dudley: Mesh: optimizes the distribution of DOFs across processors */
# Line 37  Line 36 
36     any node has no vertex).     any node has no vertex).
37   **************************************************************/   **************************************************************/
38  #ifdef USE_PARMETIS  #ifdef USE_PARMETIS
39  int Check_Inputs_For_Parmetis(dim_t mpiSize, dim_t rank, dim_t *distribution, MPI_Comm *comm)  int Check_Inputs_For_Parmetis(dim_t mpiSize, dim_t rank, dim_t * distribution, MPI_Comm * comm)
40  {  {
41    dim_t i, len;      dim_t i, len;
42    int ret_val = 1;      int ret_val = 1;
43    
44    if (rank == 0){      if (rank == 0)
45      for (i=0; i<mpiSize; i++){      {
46        len = distribution[i+1] - distribution[i];      for (i = 0; i < mpiSize; i++)
47        if (len == 0){      {
48          ret_val = 0;          len = distribution[i + 1] - distribution[i];
49          break;          if (len == 0)
50        }          {
51            ret_val = 0;
52            break;
53            }
54        }
55      }      }
56    }      MPI_Bcast(&ret_val, 1, MPI_INTEGER, 0, *comm);
57    MPI_Bcast(&ret_val, 1, MPI_INTEGER, 0, *comm);      if (ret_val == 0)
58    if (ret_val == 0)      printf("INFO: Parmetis is not used since some nodes have no vertex!\n");
59      printf("INFO: Parmetis is not used since some nodes have no vertex!\n");      return ret_val;
   return ret_val;  
60  }  }
61  #endif  #endif
62    
63    /**************************************************************/
64    
65    void Dudley_Mesh_optimizeDOFDistribution(Dudley_Mesh * in, dim_t * distribution)
66    {
67    
68  /**************************************************************/      dim_t dim, i, j, k, myNumVertices, p, mpiSize, len, globalNumVertices, *partition_count = NULL, *new_distribution =
69        NULL, *loc_partition_count = NULL;
70        bool_t *setNewDOFId = NULL;
71        index_t myFirstVertex, myLastVertex, firstVertex, lastVertex, *newGlobalDOFID = NULL;
72        size_t mpiSize_size;
73        index_t *partition = NULL;
74        Paso_Pattern *pattern = NULL;
75        Paso_MPI_rank myRank, dest, source, current_rank, rank;
76        Dudley_IndexList *index_list = NULL;
77        float *xyz = NULL;
78        int c;
79    
80  void Dudley_Mesh_optimizeDOFDistribution(Dudley_Mesh* in,dim_t *distribution) {  #ifdef PASO_MPI
81        MPI_Status status;
82    #endif
83    
84       dim_t dim, i,j,k, myNumVertices,p, mpiSize, len, globalNumVertices,*partition_count=NULL, *new_distribution=NULL, *loc_partition_count=NULL;      if (in == NULL)
85       bool_t *setNewDOFId=NULL;      return;
86       index_t myFirstVertex, myLastVertex, firstVertex, lastVertex, *newGlobalDOFID=NULL;      if (in->Nodes == NULL)
87       size_t mpiSize_size;      return;
88       index_t* partition=NULL;  
89       Paso_Pattern *pattern=NULL;      myRank = in->MPIInfo->rank;
90       Paso_MPI_rank myRank,dest,source,current_rank, rank;      mpiSize = in->MPIInfo->size;
91       Dudley_IndexList* index_list=NULL;      mpiSize_size = mpiSize * sizeof(dim_t);
92       float *xyz=NULL;      dim = in->Nodes->numDim;
93       int c;      /* first step is to distribute the elements according to a global X of DOF */
94        
95       #ifdef PASO_MPI      myFirstVertex = distribution[myRank];
96       MPI_Status status;      myLastVertex = distribution[myRank + 1];
97       #endif      myNumVertices = myLastVertex - myFirstVertex;
98        globalNumVertices = distribution[mpiSize];
99       if (in==NULL) return;      len = 0;
100       if (in->Nodes == NULL) return;      for (p = 0; p < mpiSize; ++p)
101        len = MAX(len, distribution[p + 1] - distribution[p]);
102       myRank=in->MPIInfo->rank;      partition = TMPMEMALLOC(len, index_t);  /* len is used for the sending around of partition later on */
103       mpiSize=in->MPIInfo->size;      xyz = TMPMEMALLOC(myNumVertices * dim, float);
104       mpiSize_size=mpiSize*sizeof(dim_t);      partition_count = TMPMEMALLOC(mpiSize + 1, dim_t);
105       dim=in->Nodes->numDim;      new_distribution = TMPMEMALLOC(mpiSize + 1, dim_t);
106       /* first step is to distribute the elements according to a global X of DOF */      newGlobalDOFID = TMPMEMALLOC(len, index_t);
107        setNewDOFId = TMPMEMALLOC(in->Nodes->numNodes, bool_t);
108       myFirstVertex=distribution[myRank];      if (!
109       myLastVertex=distribution[myRank+1];      (Dudley_checkPtr(partition) || Dudley_checkPtr(xyz) || Dudley_checkPtr(partition_count)
110       myNumVertices=myLastVertex-myFirstVertex;       || Dudley_checkPtr(partition_count) || Dudley_checkPtr(newGlobalDOFID) || Dudley_checkPtr(setNewDOFId)))
111       globalNumVertices=distribution[mpiSize];      {
112       len=0;      dim_t *recvbuf = TMPMEMALLOC(mpiSize * mpiSize, dim_t);
113       for (p=0;p<mpiSize;++p) len=MAX(len,distribution[p+1]-distribution[p]);  
114       partition=TMPMEMALLOC(len,index_t); /* len is used for the sending around of partition later on */      /* set the coordinates: */
115       xyz=TMPMEMALLOC(myNumVertices*dim,float);      /* it is assumed that at least one node on this processor provides a coordinate */
116       partition_count=TMPMEMALLOC(mpiSize+1,dim_t);  #pragma omp parallel for private(i,j,k)
117       new_distribution=TMPMEMALLOC(mpiSize+1,dim_t);      for (i = 0; i < in->Nodes->numNodes; ++i)
118       newGlobalDOFID=TMPMEMALLOC(len,index_t);      {
119       setNewDOFId=TMPMEMALLOC(in->Nodes->numNodes,bool_t);          k = in->Nodes->globalDegreesOfFreedom[i] - myFirstVertex;
120       if (!(Dudley_checkPtr(partition) || Dudley_checkPtr(xyz) || Dudley_checkPtr(partition_count) || Dudley_checkPtr(partition_count) || Dudley_checkPtr(newGlobalDOFID) || Dudley_checkPtr(setNewDOFId))) {          if ((k >= 0) && (k < myNumVertices))
121           dim_t *recvbuf=TMPMEMALLOC(mpiSize*mpiSize,dim_t);          {
122            for (j = 0; j < dim; ++j)
123           /* set the coordinates: */              xyz[k * dim + j] = (float)(in->Nodes->Coordinates[INDEX2(j, i, dim)]);
124           /* it is assumed that at least one node on this processor provides a coordinate */          }
125           #pragma omp parallel for private(i,j,k)      }
126           for (i=0;i<in->Nodes->numNodes;++i) {  
127               k=in->Nodes->globalDegreesOfFreedom[i]-myFirstVertex;      index_list = TMPMEMALLOC(myNumVertices, Dudley_IndexList);
128               if ((k>=0) && (k<myNumVertices)) {      /* ksteube CSR of DOF IDs */
129                  for (j=0;j<dim;++j) xyz[k*dim+j]=(float)(in->Nodes->Coordinates[INDEX2(j,i,dim)]);      /* create the adjacency structure xadj and adjncy */
130               }      if (!Dudley_checkPtr(index_list))
131           }      {
132    #pragma omp parallel private(i)
133           index_list=TMPMEMALLOC(myNumVertices,Dudley_IndexList);          {
134       /* ksteube CSR of DOF IDs */  #pragma omp for schedule(static)
135           /* create the adjacency structure xadj and adjncy */          for (i = 0; i < myNumVertices; ++i)
136           if (! Dudley_checkPtr(index_list)) {          {
137              #pragma omp parallel private(i)              index_list[i].extension = NULL;
138              {              index_list[i].n = 0;
139                #pragma omp for schedule(static)          }
140                for(i=0;i<myNumVertices;++i) {          /* ksteube build CSR format */
141                     index_list[i].extension=NULL;          /*  insert contributions from element matrices into colums index index_list: */
142                     index_list[i].n=0;          Dudley_IndexList_insertElementsWithRowRangeNoMainDiagonal(index_list, myFirstVertex, myLastVertex,
143                }                                        in->Elements,
144            /* ksteube build CSR format */                                        in->Nodes->globalDegreesOfFreedom,
145                /*  insert contributions from element matrices into colums index index_list: */                                        in->Nodes->globalDegreesOfFreedom);
146                Dudley_IndexList_insertElementsWithRowRangeNoMainDiagonal(index_list, myFirstVertex, myLastVertex,          Dudley_IndexList_insertElementsWithRowRangeNoMainDiagonal(index_list, myFirstVertex, myLastVertex,
147                                                                          in->Elements,in->Nodes->globalDegreesOfFreedom,                                        in->FaceElements,
148                                                                          in->Nodes->globalDegreesOfFreedom);                                        in->Nodes->globalDegreesOfFreedom,
149                Dudley_IndexList_insertElementsWithRowRangeNoMainDiagonal(index_list, myFirstVertex, myLastVertex,                                        in->Nodes->globalDegreesOfFreedom);
150                                                                          in->FaceElements,in->Nodes->globalDegreesOfFreedom,          Dudley_IndexList_insertElementsWithRowRangeNoMainDiagonal(index_list, myFirstVertex, myLastVertex,
151                                                                          in->Nodes->globalDegreesOfFreedom);                                        in->Points, in->Nodes->globalDegreesOfFreedom,
152                Dudley_IndexList_insertElementsWithRowRangeNoMainDiagonal(index_list, myFirstVertex, myLastVertex,                                        in->Nodes->globalDegreesOfFreedom);
153                                                                          in->Points,in->Nodes->globalDegreesOfFreedom,          }
154                                                                          in->Nodes->globalDegreesOfFreedom);  
155             }          /* create the local matrix pattern */
156                      pattern = Dudley_IndexList_createPattern(0, myNumVertices, index_list, 0, globalNumVertices, 0);
157             /* create the local matrix pattern */  
158             pattern=Dudley_IndexList_createPattern(0,myNumVertices,index_list,0,globalNumVertices,0);          /* clean up index list */
159            if (index_list != NULL)
160             /* clean up index list */          {
161             if (index_list!=NULL) {  #pragma omp parallel for private(i)
162                #pragma omp parallel for private(i)          for (i = 0; i < myNumVertices; ++i)
163                for(i=0;i<myNumVertices;++i) Dudley_IndexList_free(index_list[i].extension);              Dudley_IndexList_free(index_list[i].extension);
164             }          }
165    
166             if (Dudley_noError()) {          if (Dudley_noError())
167            {
168    
169  #ifdef USE_PARMETIS  #ifdef USE_PARMETIS
170    
171            if (mpiSize>1 &&          if (mpiSize > 1 && Check_Inputs_For_Parmetis(mpiSize, myRank, distribution, &(in->MPIInfo->comm)) > 0)
172            Check_Inputs_For_Parmetis(mpiSize, myRank, distribution, &(in->MPIInfo->comm))>0 ) {          {
173           int i;              int i;
174           int wgtflag = 0;              int wgtflag = 0;
175           int numflag = 0;   /* pattern->ptr is C style: starting from 0 instead of 1 */              int numflag = 0;    /* pattern->ptr is C style: starting from 0 instead of 1 */
176           int ncon = 1;              int ncon = 1;
177           int edgecut;              int edgecut;
178           int options[2];              int options[2];
179           float *tpwgts = TMPMEMALLOC(ncon*mpiSize,float);              float *tpwgts = TMPMEMALLOC(ncon * mpiSize, float);
180           float *ubvec = TMPMEMALLOC(ncon,float);              float *ubvec = TMPMEMALLOC(ncon, float);
181           for (i=0; i<ncon*mpiSize; i++) tpwgts[i] = 1.0/(float)mpiSize;              for (i = 0; i < ncon * mpiSize; i++)
182           for (i=0; i<ncon; i++) ubvec[i] = 1.05;              tpwgts[i] = 1.0 / (float)mpiSize;
183           options[0] = 3;              for (i = 0; i < ncon; i++)
184           options[1] = 15;              ubvec[i] = 1.05;
185               ParMETIS_V3_PartGeomKway(distribution,              options[0] = 3;
186                                            pattern->ptr,              options[1] = 15;
187                                            pattern->index,              ParMETIS_V3_PartGeomKway(distribution, pattern->ptr, pattern->index, NULL, NULL, &wgtflag, &numflag, &dim, xyz, &ncon, &mpiSize, tpwgts, ubvec, options, &edgecut, partition,   /* new CPU ownership of elements */
188                                            NULL,                           &(in->MPIInfo->comm));
189                                            NULL,              /* printf("ParMETIS number of edges cut by partitioning per processor: %d\n", edgecut/MAX(in->MPIInfo->size,1)); */
190                                            &wgtflag,              TMPMEMFREE(ubvec);
191                                            &numflag,              TMPMEMFREE(tpwgts);
192                                            &dim,          }
193                                            xyz,          else
194                                            &ncon,          {
195                                            &mpiSize,              for (i = 0; i < myNumVertices; ++i)
196                                            tpwgts,              partition[i] = 0;   /* CPU 0 owns it */
197                                            ubvec,          }
                                           options,  
                                           &edgecut,  
                                           partition,                /* new CPU ownership of elements */  
                                           &(in->MPIInfo->comm));  
          /* printf("ParMETIS number of edges cut by partitioning per processor: %d\n", edgecut/MAX(in->MPIInfo->size,1)); */  
                  TMPMEMFREE(ubvec);  
                  TMPMEMFREE(tpwgts);  
           } else {  
                  for (i=0;i<myNumVertices;++i) partition[i]=0;      /* CPU 0 owns it */  
           }  
198  #else  #else
199                for (i=0;i<myNumVertices;++i) partition[i]=myRank;    /* CPU 0 owns it */          for (i = 0; i < myNumVertices; ++i)
200                partition[i] = myRank;  /* CPU 0 owns it */
201  #endif  #endif
202    
203             }          }
204    
205             Paso_Pattern_free(pattern);          Paso_Pattern_free(pattern);
206              
207             /* create a new distributioin and labeling of the DOF */          /* create a new distributioin and labeling of the DOF */
208             memset(new_distribution,0,mpiSize_size);          memset(new_distribution, 0, mpiSize_size);
209             #pragma omp parallel private(loc_partition_count)  #pragma omp parallel private(loc_partition_count)
210             {          {
211                 loc_partition_count=THREAD_MEMALLOC(mpiSize,dim_t);          loc_partition_count = THREAD_MEMALLOC(mpiSize, dim_t);
212                 memset(loc_partition_count,0,mpiSize_size);          memset(loc_partition_count, 0, mpiSize_size);
213                 #pragma omp for private(i)  #pragma omp for private(i)
214                 for (i=0;i<myNumVertices;++i) loc_partition_count[partition[i]]++ ;          for (i = 0; i < myNumVertices; ++i)
215                 #pragma omp critical              loc_partition_count[partition[i]]++;
216                 {  #pragma omp critical
217                    for (i=0;i<mpiSize;++i) new_distribution[i]+=loc_partition_count[i];          {
218                 }              for (i = 0; i < mpiSize; ++i)
219                 THREAD_MEMFREE(loc_partition_count);              new_distribution[i] += loc_partition_count[i];
220             }          }
221             #ifdef PASO_MPI          THREAD_MEMFREE(loc_partition_count);
222            /* recvbuf will be the concatenation of each CPU's contribution to new_distribution */          }
223            MPI_Allgather(new_distribution, mpiSize, MPI_INT, recvbuf, mpiSize, MPI_INT, in->MPIInfo->comm);  #ifdef PASO_MPI
224             #else          /* recvbuf will be the concatenation of each CPU's contribution to new_distribution */
225                 for (i=0;i<mpiSize;++i) recvbuf[i]=new_distribution[i];          MPI_Allgather(new_distribution, mpiSize, MPI_INT, recvbuf, mpiSize, MPI_INT, in->MPIInfo->comm);
226             #endif  #else
227             new_distribution[0]=0;          for (i = 0; i < mpiSize; ++i)
228         for (rank=0; rank<mpiSize;rank++) {          recvbuf[i] = new_distribution[i];
229            c=0;  #endif
230                for (i=0;i<myRank;++i) c+=recvbuf[rank+mpiSize*i];          new_distribution[0] = 0;
231                for (i=0;i<myNumVertices;++i) {          for (rank = 0; rank < mpiSize; rank++)
232                   if (rank==partition[i]) {          {
233                      newGlobalDOFID[i]=new_distribution[rank]+c;          c = 0;
234                      c++;          for (i = 0; i < myRank; ++i)
235               }              c += recvbuf[rank + mpiSize * i];
236                }          for (i = 0; i < myNumVertices; ++i)
237                for (i=myRank+1;i<mpiSize;++i) c+=recvbuf[rank+mpiSize*i];          {
238                new_distribution[rank+1]=new_distribution[rank]+c;              if (rank == partition[i])
239             }              {
240             TMPMEMFREE(recvbuf);              newGlobalDOFID[i] = new_distribution[rank] + c;
241                c++;
242             /* now the overlap needs to be created by sending the partition around*/              }
243            }
244             dest=Paso_MPIInfo_mod(mpiSize, myRank + 1);          for (i = myRank + 1; i < mpiSize; ++i)
245             source=Paso_MPIInfo_mod(mpiSize, myRank - 1);              c += recvbuf[rank + mpiSize * i];
246             current_rank=myRank;          new_distribution[rank + 1] = new_distribution[rank] + c;
247             #pragma omp parallel for private(i)          }
248             for (i=0;i<in->Nodes->numNodes;++i) setNewDOFId[i]=TRUE;          TMPMEMFREE(recvbuf);
249    
250             for (p=0; p< mpiSize; ++p) {          /* now the overlap needs to be created by sending the partition around */
251    
252                 firstVertex=distribution[current_rank];          dest = Paso_MPIInfo_mod(mpiSize, myRank + 1);
253                 lastVertex=distribution[current_rank+1];          source = Paso_MPIInfo_mod(mpiSize, myRank - 1);
254                 #pragma omp parallel for private(i,j,k)          current_rank = myRank;
255                 for (i=0;i<in->Nodes->numNodes;++i) {  #pragma omp parallel for private(i)
256                     k=in->Nodes->globalDegreesOfFreedom[i];          for (i = 0; i < in->Nodes->numNodes; ++i)
257                     if (setNewDOFId[i] && (firstVertex<=k) && (k<lastVertex)) {          setNewDOFId[i] = TRUE;
258                          in->Nodes->globalDegreesOfFreedom[i]=newGlobalDOFID[k-firstVertex];  
259                          setNewDOFId[i]=FALSE;          for (p = 0; p < mpiSize; ++p)
260                     }          {
261                 }  
262            firstVertex = distribution[current_rank];
263                 if (p<mpiSize-1) {  /* the final send can be skipped */          lastVertex = distribution[current_rank + 1];
264                    #ifdef PASO_MPI  #pragma omp parallel for private(i,j,k)
265                    MPI_Sendrecv_replace(newGlobalDOFID,len, MPI_INT,          for (i = 0; i < in->Nodes->numNodes; ++i)
266                                         dest, in->MPIInfo->msg_tag_counter,          {
267                                         source, in->MPIInfo->msg_tag_counter,              k = in->Nodes->globalDegreesOfFreedom[i];
268                                         in->MPIInfo->comm,&status);              if (setNewDOFId[i] && (firstVertex <= k) && (k < lastVertex))
269                    #endif              {
270                    in->MPIInfo->msg_tag_counter++;              in->Nodes->globalDegreesOfFreedom[i] = newGlobalDOFID[k - firstVertex];
271                    current_rank=Paso_MPIInfo_mod(mpiSize, current_rank-1);              setNewDOFId[i] = FALSE;
272                }              }
273             }          }
274             for (i=0;i<mpiSize+1;++i) distribution[i]=new_distribution[i];  
275           }          if (p < mpiSize - 1)
276           TMPMEMFREE(index_list);          {       /* the final send can be skipped */
277       }  #ifdef PASO_MPI
278       TMPMEMFREE(newGlobalDOFID);              MPI_Sendrecv_replace(newGlobalDOFID, len, MPI_INT,
279       TMPMEMFREE(setNewDOFId);                       dest, in->MPIInfo->msg_tag_counter,
280       TMPMEMFREE(new_distribution);                       source, in->MPIInfo->msg_tag_counter, in->MPIInfo->comm, &status);
281       TMPMEMFREE(partition_count);  #endif
282       TMPMEMFREE(partition);              in->MPIInfo->msg_tag_counter++;
283       TMPMEMFREE(xyz);              current_rank = Paso_MPIInfo_mod(mpiSize, current_rank - 1);
284       return;          }
285            }
286            for (i = 0; i < mpiSize + 1; ++i)
287            distribution[i] = new_distribution[i];
288        }
289        TMPMEMFREE(index_list);
290        }
291        TMPMEMFREE(newGlobalDOFID);
292        TMPMEMFREE(setNewDOFId);
293        TMPMEMFREE(new_distribution);
294        TMPMEMFREE(partition_count);
295        TMPMEMFREE(partition);
296        TMPMEMFREE(xyz);
297        return;
298  }  }

Legend:
Removed from v.3126  
changed lines
  Added in v.3224

  ViewVC Help
Powered by ViewVC 1.1.26