/[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

trunk/finley/src/Mesh_optimizeDOFDistribution.c revision 1763 by gross, Mon Sep 8 02:47:55 2008 UTC branches/domexper/dudley/src/Mesh_optimizeDOFDistribution.c revision 3224 by jfenwick, Wed Sep 29 05:19:37 2010 UTC
# Line 1  Line 1 
1    
 /* $Id$ */  
   
2  /*******************************************************  /*******************************************************
3   *  *
4   *           Copyright 2003-2007 by ACceSS MNRF  * Copyright (c) 2003-2010 by University of Queensland
5   *       Copyright 2007 by University of Queensland  * Earth Systems Science Computational Center (ESSCC)
6   *  * http://www.uq.edu.au/esscc
7   *                http://esscc.uq.edu.au  *
8   *        Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
9   *  Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
10   *     http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
11   *  *
12   *******************************************************/  *******************************************************/
13    
14  /**************************************************************/  /**************************************************************/
15    
16  /*   Finley: Mesh: optimizes the distribution of DOFs across processors */  /*   Dudley: Mesh: optimizes the distribution of DOFs across processors */
17  /*   using ParMETIS. On return a new distribution is given and the globalDOF are relabled */  /*   using ParMETIS. On return a new distribution is given and the globalDOF are relabled */
18  /*   accordingly but the mesh has not been redesitributed yet                             */  /*   accordingly but the mesh has not been redesitributed yet                             */
19    
# Line 30  Line 28 
28  #include "parmetis.h"  #include "parmetis.h"
29  #endif  #endif
30    
31  /**************************************************************/  /**************************************************************
32       Check whether there is any node which has no vertex. In case
33       such node exists, we don't use parmetis since parmetis requires
34       that every node has at least 1 vertex (at line 129 of file
35       "xyzpart.c" in parmetis 3.1.1, variable "nvtxs" would be 0 if
36       any node has no vertex).
37     **************************************************************/
38    #ifdef USE_PARMETIS
39    int Check_Inputs_For_Parmetis(dim_t mpiSize, dim_t rank, dim_t * distribution, MPI_Comm * comm)
40    {
41        dim_t i, len;
42        int ret_val = 1;
43    
44  void Finley_Mesh_optimizeDOFDistribution(Finley_Mesh* in,dim_t *distribution) {      if (rank == 0)
45        {
46        for (i = 0; i < mpiSize; i++)
47        {
48            len = distribution[i + 1] - distribution[i];
49            if (len == 0)
50            {
51            ret_val = 0;
52            break;
53            }
54        }
55        }
56        MPI_Bcast(&ret_val, 1, MPI_INTEGER, 0, *comm);
57        if (ret_val == 0)
58        printf("INFO: Parmetis is not used since some nodes have no vertex!\n");
59        return ret_val;
60    }
61    #endif
62    
63       dim_t dim, i,j,k, myNumVertices,p, mpiSize, len, globalNumVertices,*partition_count=NULL, *new_distribution=NULL, *loc_partition_count=NULL;  /**************************************************************/
      bool_t *setNewDOFId=NULL;  
      index_t myFirstVertex, myLastVertex, firstVertex, lastVertex, *newGlobalDOFID=NULL;  
      size_t mpiSize_size;  
      index_t* partition=NULL;  
      Paso_Pattern *pattern=NULL;  
      Paso_MPI_rank myRank,dest,source,current_rank, rank;  
      Finley_IndexList* index_list=NULL;  
      float *xyz=NULL;  
      int c;  
       
      #ifdef PASO_MPI  
      MPI_Status status;  
      #endif  
   
      if (in==NULL) return;  
      if (in->Nodes == NULL) return;  
   
      myRank=in->MPIInfo->rank;  
      mpiSize=in->MPIInfo->size;  
      mpiSize_size=mpiSize*sizeof(dim_t);  
      dim=in->Nodes->numDim;  
      /* first step is to distribute the elements according to a global X of DOF */  
   
      myFirstVertex=distribution[myRank];  
      myLastVertex=distribution[myRank+1];  
      myNumVertices=myLastVertex-myFirstVertex;  
      globalNumVertices=distribution[mpiSize];  
      len=0;  
      for (p=0;p<mpiSize;++p) len=MAX(len,distribution[p+1]-distribution[p]);  
      partition=TMPMEMALLOC(len,index_t); /* len is used for the sending around of partition later on */  
      xyz=TMPMEMALLOC(myNumVertices*dim,float);  
      partition_count=TMPMEMALLOC(mpiSize+1,dim_t);  
      new_distribution=TMPMEMALLOC(mpiSize+1,dim_t);  
      newGlobalDOFID=TMPMEMALLOC(len,index_t);  
      setNewDOFId=TMPMEMALLOC(in->Nodes->numNodes,bool_t);  
      if (!(Finley_checkPtr(partition) || Finley_checkPtr(xyz) || Finley_checkPtr(partition_count) || Finley_checkPtr(partition_count) || Finley_checkPtr(newGlobalDOFID) || Finley_checkPtr(setNewDOFId))) {  
          dim_t *recvbuf=TMPMEMALLOC(mpiSize*mpiSize,dim_t);  
   
          /* set the coordinates: */  
          /* it is assumed that at least one node on this processor provides a coordinate */  
          #pragma omp parallel for private(i,j,k)  
          for (i=0;i<in->Nodes->numNodes;++i) {  
              k=in->Nodes->globalDegreesOfFreedom[i]-myFirstVertex;  
              if ((k>=0) && (k<myNumVertices)) {  
                 for (j=0;j<dim;++j) xyz[k*dim+j]=(float)(in->Nodes->Coordinates[INDEX2(j,i,dim)]);  
              }  
          }  
   
          index_list=TMPMEMALLOC(myNumVertices,Finley_IndexList);  
      /* ksteube CSR of DOF IDs */  
          /* create the adjacency structure xadj and adjncy */  
          if (! Finley_checkPtr(index_list)) {  
             #pragma omp parallel private(i)  
             {  
               #pragma omp for schedule(static)  
               for(i=0;i<myNumVertices;++i) {  
                    index_list[i].extension=NULL;  
                    index_list[i].n=0;  
               }  
           /* ksteube build CSR format */  
               /*  insert contributions from element matrices into colums index index_list: */  
               Finley_IndexList_insertElementsWithRowRangeNoMainDiagonal(index_list, myFirstVertex, myLastVertex,  
                                                                         in->Elements,in->Nodes->globalDegreesOfFreedom,  
                                                                         in->Nodes->globalDegreesOfFreedom);  
               Finley_IndexList_insertElementsWithRowRangeNoMainDiagonal(index_list, myFirstVertex, myLastVertex,  
                                                                         in->FaceElements,in->Nodes->globalDegreesOfFreedom,  
                                                                         in->Nodes->globalDegreesOfFreedom);  
               Finley_IndexList_insertElementsWithRowRangeNoMainDiagonal(index_list, myFirstVertex, myLastVertex,  
                                                                         in->ContactElements,in->Nodes->globalDegreesOfFreedom,  
                                                                         in->Nodes->globalDegreesOfFreedom);  
               Finley_IndexList_insertElementsWithRowRangeNoMainDiagonal(index_list, myFirstVertex, myLastVertex,  
                                                                         in->Points,in->Nodes->globalDegreesOfFreedom,  
                                                                         in->Nodes->globalDegreesOfFreedom);  
            }  
             
            /* create the local matrix pattern */  
            pattern=Finley_IndexList_createPattern(0,myNumVertices,index_list,0,globalNumVertices,0);  
   
            /* clean up index list */  
            if (index_list!=NULL) {  
               #pragma omp parallel for private(i)  
               for(i=0;i<myNumVertices;++i) Finley_IndexList_free(index_list[i].extension);  
            }  
64    
65             if (Finley_noError()) {  void Dudley_Mesh_optimizeDOFDistribution(Dudley_Mesh * in, dim_t * distribution)
66    {
67    
68  #ifdef USE_PARMETIS      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            if (in->MPIInfo->size>1) {  #ifdef PASO_MPI
81           int i;      MPI_Status status;
82           int wgtflag = 0;  #endif
          int numflag = 0;   /* pattern->ptr is C style: starting from 0 instead of 1 */  
          int ncon = 1;  
          int edgecut;  
          int options[2];  
          float *tpwgts = TMPMEMALLOC(ncon*mpiSize,float);  
          float *ubvec = TMPMEMALLOC(ncon,float);  
          for (i=0; i<ncon*mpiSize; i++) tpwgts[i] = 1.0/(float)mpiSize;  
          for (i=0; i<ncon; i++) ubvec[i] = 1.05;  
          options[0] = 3;  
          options[1] = 15;  
83    
84  {      if (in == NULL)
85  int k=distribution[in->MPIInfo->rank+1]-distribution[in->MPIInfo->rank];      return;
86  int min_i,max_i;      if (in->Nodes == NULL)
87  printf("INPUT PARMETIS: %d\n",k);      return;
88  for(i=0;i<in->MPIInfo->size+1;++i) printf("%d ",distribution[i]);  
89  printf("\n");      myRank = in->MPIInfo->rank;
90  min_i=pattern->index[0];      mpiSize = in->MPIInfo->size;
91  max_i=pattern->index[0];      mpiSize_size = mpiSize * sizeof(dim_t);
92  for(i=0;i<pattern->ptr[k];++i) {      dim = in->Nodes->numDim;
93  min_i=MIN(min_i,pattern->index[i]);      /* first step is to distribute the elements according to a global X of DOF */
94  max_i=MAX(max_i,pattern->index[i]);  
95  }      myFirstVertex = distribution[myRank];
96  printf("index range = %d : %d\n",min_i,max_i);      myLastVertex = distribution[myRank + 1];
97        myNumVertices = myLastVertex - myFirstVertex;
98        globalNumVertices = distribution[mpiSize];
99        len = 0;
100        for (p = 0; p < mpiSize; ++p)
101        len = MAX(len, distribution[p + 1] - distribution[p]);
102        partition = TMPMEMALLOC(len, index_t);  /* len is used for the sending around of partition later on */
103        xyz = TMPMEMALLOC(myNumVertices * dim, float);
104        partition_count = TMPMEMALLOC(mpiSize + 1, dim_t);
105        new_distribution = TMPMEMALLOC(mpiSize + 1, dim_t);
106        newGlobalDOFID = TMPMEMALLOC(len, index_t);
107        setNewDOFId = TMPMEMALLOC(in->Nodes->numNodes, bool_t);
108        if (!
109        (Dudley_checkPtr(partition) || Dudley_checkPtr(xyz) || Dudley_checkPtr(partition_count)
110         || Dudley_checkPtr(partition_count) || Dudley_checkPtr(newGlobalDOFID) || Dudley_checkPtr(setNewDOFId)))
111        {
112        dim_t *recvbuf = TMPMEMALLOC(mpiSize * mpiSize, dim_t);
113    
114        /* set the coordinates: */
115        /* it is assumed that at least one node on this processor provides a coordinate */
116    #pragma omp parallel for private(i,j,k)
117        for (i = 0; i < in->Nodes->numNodes; ++i)
118        {
119            k = in->Nodes->globalDegreesOfFreedom[i] - myFirstVertex;
120            if ((k >= 0) && (k < myNumVertices))
121            {
122            for (j = 0; j < dim; ++j)
123                xyz[k * dim + j] = (float)(in->Nodes->Coordinates[INDEX2(j, i, dim)]);
124            }
125        }
126    
127        index_list = TMPMEMALLOC(myNumVertices, Dudley_IndexList);
128        /* ksteube CSR of DOF IDs */
129        /* create the adjacency structure xadj and adjncy */
130        if (!Dudley_checkPtr(index_list))
131        {
132    #pragma omp parallel private(i)
133            {
134    #pragma omp for schedule(static)
135            for (i = 0; i < myNumVertices; ++i)
136            {
137                index_list[i].extension = NULL;
138                index_list[i].n = 0;
139            }
140            /* ksteube build CSR format */
141            /*  insert contributions from element matrices into colums index index_list: */
142            Dudley_IndexList_insertElementsWithRowRangeNoMainDiagonal(index_list, myFirstVertex, myLastVertex,
143                                          in->Elements,
144                                          in->Nodes->globalDegreesOfFreedom,
145                                          in->Nodes->globalDegreesOfFreedom);
146            Dudley_IndexList_insertElementsWithRowRangeNoMainDiagonal(index_list, myFirstVertex, myLastVertex,
147                                          in->FaceElements,
148                                          in->Nodes->globalDegreesOfFreedom,
149                                          in->Nodes->globalDegreesOfFreedom);
150            Dudley_IndexList_insertElementsWithRowRangeNoMainDiagonal(index_list, myFirstVertex, myLastVertex,
151                                          in->Points, in->Nodes->globalDegreesOfFreedom,
152                                          in->Nodes->globalDegreesOfFreedom);
153            }
154    
155            /* create the local matrix pattern */
156            pattern = Dudley_IndexList_createPattern(0, myNumVertices, index_list, 0, globalNumVertices, 0);
157    
158            /* clean up index list */
159            if (index_list != NULL)
160            {
161    #pragma omp parallel for private(i)
162            for (i = 0; i < myNumVertices; ++i)
163                Dudley_IndexList_free(index_list[i].extension);
164            }
165    
166  for(i=0;i<k+1;++i) printf("%d ",pattern->ptr[i]);          if (Dudley_noError())
167  printf("\n");          {
168  /*  
169  for(i=0;i<pattern->ptr[k];++i) printf("%d ",pattern->index[i]);  #ifdef USE_PARMETIS
 printf("\n");  
 */  
 }  
170    
171           ParMETIS_V3_PartGeomKway(distribution,          if (mpiSize > 1 && Check_Inputs_For_Parmetis(mpiSize, myRank, distribution, &(in->MPIInfo->comm)) > 0)
172                                   pattern->ptr,          {
173                                   pattern->index,              int i;
174                                   NULL,              int wgtflag = 0;
175                                   NULL,              int numflag = 0;    /* pattern->ptr is C style: starting from 0 instead of 1 */
176                                   &wgtflag,              int ncon = 1;
177                                   &numflag,              int edgecut;
178                                   &dim,              int options[2];
179                                   xyz,              float *tpwgts = TMPMEMALLOC(ncon * mpiSize, float);
180                                   &ncon,              float *ubvec = TMPMEMALLOC(ncon, float);
181                                   &mpiSize,              for (i = 0; i < ncon * mpiSize; i++)
182                                   tpwgts,              tpwgts[i] = 1.0 / (float)mpiSize;
183                                   ubvec,              for (i = 0; i < ncon; i++)
184                                   options,              ubvec[i] = 1.05;
185                                   &edgecut,              options[0] = 3;
186                                   partition,             /* new CPU ownership of elements */              options[1] = 15;
187                                   &(in->MPIInfo->comm));              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           printf("ParMETIS number of edges cut by partitioning per processor: %d\n", edgecut/MAX(in->MPIInfo->size,1));                           &(in->MPIInfo->comm));
189                   TMPMEMFREE(ubvec);              /* printf("ParMETIS number of edges cut by partitioning per processor: %d\n", edgecut/MAX(in->MPIInfo->size,1)); */
190                   TMPMEMFREE(tpwgts);              TMPMEMFREE(ubvec);
191            } else {              TMPMEMFREE(tpwgts);
192                   for (i=0;i<myNumVertices;++i) partition[i]=0;      /* CPU 0 owns it */          }
193            }          else
194            {
195                for (i = 0; i < myNumVertices; ++i)
196                partition[i] = 0;   /* CPU 0 owns it */
197            }
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.1763  
changed lines
  Added in v.3224

  ViewVC Help
Powered by ViewVC 1.1.26