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

temp_trunk_copy/finley/src/Mesh_optimizeDOFDistribution.c revision 1384 by phornby, Fri Jan 11 02:29:38 2008 UTC branches/domexper/dudley/src/Mesh_optimizeDOFDistribution.c revision 3231 by jfenwick, Fri Oct 1 01:53:46 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 26  Line 24 
24  #ifdef _OPENMP  #ifdef _OPENMP
25  #include <omp.h>  #include <omp.h>
26  #endif  #endif
27    #ifdef USE_PARMETIS
28    #include "parmetis.h"
29    #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        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  /**************************************************************/  /**************************************************************/
64    
65  void Finley_Mesh_optimizeDOFDistribution(Finley_Mesh* in,dim_t *distribution) {  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        Esys_MPI_rank myRank, dest, source, current_rank, rank;
76        Dudley_IndexList *index_list = NULL;
77        float *xyz = NULL;
78        int c;
79    
80    #ifdef ESYS_MPI
81        MPI_Status status;
82    #endif
83    
84        if (in == NULL)
85        return;
86        if (in->Nodes == NULL)
87        return;
88    
89        myRank = in->MPIInfo->rank;
90        mpiSize = in->MPIInfo->size;
91        mpiSize_size = mpiSize * sizeof(dim_t);
92        dim = in->Nodes->numDim;
93        /* first step is to distribute the elements according to a global X of DOF */
94    
95        myFirstVertex = distribution[myRank];
96        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            if (Dudley_noError())
167            {
168    
169    #ifdef USE_PARMETIS
170    
171            if (mpiSize > 1 && Check_Inputs_For_Parmetis(mpiSize, myRank, distribution, &(in->MPIInfo->comm)) > 0)
172            {
173                int i;
174                int wgtflag = 0;
175                int numflag = 0;    /* pattern->ptr is C style: starting from 0 instead of 1 */
176                int ncon = 1;
177                int edgecut;
178                int options[2];
179                float *tpwgts = TMPMEMALLOC(ncon * mpiSize, float);
180                float *ubvec = TMPMEMALLOC(ncon, float);
181                for (i = 0; i < ncon * mpiSize; i++)
182                tpwgts[i] = 1.0 / (float)mpiSize;
183                for (i = 0; i < ncon; i++)
184                ubvec[i] = 1.05;
185                options[0] = 3;
186                options[1] = 15;
187                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                             &(in->MPIInfo->comm));
189                /* printf("ParMETIS number of edges cut by partitioning per processor: %d\n", edgecut/MAX(in->MPIInfo->size,1)); */
190                TMPMEMFREE(ubvec);
191                TMPMEMFREE(tpwgts);
192            }
193            else
194            {
195                for (i = 0; i < myNumVertices; ++i)
196                partition[i] = 0;   /* CPU 0 owns it */
197            }
198    #else
199            for (i = 0; i < myNumVertices; ++i)
200                partition[i] = myRank;  /* CPU 0 owns it */
201    #endif
202    
203            }
204    
205       dim_t dim, i,j,k, myNumVertices,p, mpiSize, len, globalNumVertices,*partition_count=NULL, *new_distribution=NULL, *loc_partition_count=NULL;          Paso_Pattern_free(pattern);
206       bool_t *setNewDOFId=NULL;  
207       index_t myFirstVertex, myLastVertex, firstVertex, lastVertex, *newGlobalDOFID=NULL;          /* create a new distributioin and labeling of the DOF */
208       size_t mpiSize_size;          memset(new_distribution, 0, mpiSize_size);
209       index_t* partition=NULL;  #pragma omp parallel private(loc_partition_count)
210       Paso_Pattern *pattern=NULL;          {
211       Paso_MPI_rank myRank,dest,source,current_rank, rank;          loc_partition_count = THREAD_MEMALLOC(mpiSize, dim_t);
212       Finley_IndexList* index_list=NULL;          memset(loc_partition_count, 0, mpiSize_size);
213       float *xyz=NULL;  #pragma omp for private(i)
214                for (i = 0; i < myNumVertices; ++i)
215       #ifdef PASO_MPI              loc_partition_count[partition[i]]++;
216       MPI_Status status;  #pragma omp critical
217       #endif          {
218                for (i = 0; i < mpiSize; ++i)
219       if (in==NULL) return;              new_distribution[i] += loc_partition_count[i];
220       if (in->Nodes == NULL) return;          }
221            THREAD_MEMFREE(loc_partition_count);
222       myRank=in->MPIInfo->rank;          }
223       mpiSize=in->MPIInfo->size;  #ifdef ESYS_MPI
224       mpiSize_size=mpiSize*sizeof(dim_t);          /* recvbuf will be the concatenation of each CPU's contribution to new_distribution */
225       dim=in->Nodes->numDim;          MPI_Allgather(new_distribution, mpiSize, MPI_INT, recvbuf, mpiSize, MPI_INT, in->MPIInfo->comm);
226       /* first step is to distribute the elements according to a global X of DOF */  #else
227            for (i = 0; i < mpiSize; ++i)
228       myFirstVertex=distribution[myRank];          recvbuf[i] = new_distribution[i];
229       myLastVertex=distribution[myRank+1];  #endif
230       myNumVertices=myLastVertex-myFirstVertex;          new_distribution[0] = 0;
231       globalNumVertices=distribution[mpiSize];          for (rank = 0; rank < mpiSize; rank++)
232       len=0;          {
233       for (p=0;p<mpiSize;++p) len=MAX(len,distribution[p+1]-distribution[p]);          c = 0;
234       partition=TMPMEMALLOC(len,index_t); /* len is used for the sending around of partition later on */          for (i = 0; i < myRank; ++i)
235       xyz=TMPMEMALLOC(myNumVertices*dim,float);              c += recvbuf[rank + mpiSize * i];
236       partition_count=TMPMEMALLOC(mpiSize+1,dim_t);          for (i = 0; i < myNumVertices; ++i)
237       new_distribution=TMPMEMALLOC(mpiSize+1,dim_t);          {
238       newGlobalDOFID=TMPMEMALLOC(len,index_t);              if (rank == partition[i])
239       setNewDOFId=TMPMEMALLOC(in->Nodes->numNodes,bool_t);              {
240       if (!(Finley_checkPtr(partition) || Finley_checkPtr(xyz) || Finley_checkPtr(partition_count) || Finley_checkPtr(partition_count) || Finley_checkPtr(newGlobalDOFID) || Finley_checkPtr(setNewDOFId))) {              newGlobalDOFID[i] = new_distribution[rank] + c;
241                c++;
242           /* set the coordinates: *?              }
243           /* it is assumed that at least one node on this processor provides a coordinate */          }
244           #pragma omp parallel for private(i,j,k)          for (i = myRank + 1; i < mpiSize; ++i)
245           for (i=0;i<in->Nodes->numNodes;++i) {              c += recvbuf[rank + mpiSize * i];
246               k=in->Nodes->globalDegreesOfFreedom[i]-myFirstVertex;          new_distribution[rank + 1] = new_distribution[rank] + c;
247               if ((k>=0) && (k<myNumVertices)) {          }
248                  for (j=0;j<dim;++j) xyz[k*dim+j]=(float)(in->Nodes->Coordinates[INDEX2(j,i,dim)]);          TMPMEMFREE(recvbuf);
249               }  
250           }          /* now the overlap needs to be created by sending the partition around */
251    
252           index_list=TMPMEMALLOC(myNumVertices,Finley_IndexList);          dest = Esys_MPIInfo_mod(mpiSize, myRank + 1);
253       /* ksteube CSR of DOF IDs */          source = Esys_MPIInfo_mod(mpiSize, myRank - 1);
254           /* create the adjacency structure xadj and adjncy */          current_rank = myRank;
255           if (! Finley_checkPtr(index_list)) {  #pragma omp parallel for private(i)
256              #pragma omp parallel private(i)          for (i = 0; i < in->Nodes->numNodes; ++i)
257              {          setNewDOFId[i] = TRUE;
258                #pragma omp for schedule(static)  
259                for(i=0;i<myNumVertices;++i) {          for (p = 0; p < mpiSize; ++p)
260                     index_list[i].extension=NULL;          {
261                     index_list[i].n=0;  
262                }          firstVertex = distribution[current_rank];
263            /* ksteube build CSR format */          lastVertex = distribution[current_rank + 1];
264                /*  insert contributions from element matrices into colums index index_list: */  #pragma omp parallel for private(i,j,k)
265                Finley_IndexList_insertElementsWithRowRange(index_list, myFirstVertex, myLastVertex,          for (i = 0; i < in->Nodes->numNodes; ++i)
266                                                            in->Elements,in->Nodes->globalDegreesOfFreedom,          {
267                                                            in->Nodes->globalDegreesOfFreedom);              k = in->Nodes->globalDegreesOfFreedom[i];
268                Finley_IndexList_insertElementsWithRowRange(index_list, myFirstVertex, myLastVertex,              if (setNewDOFId[i] && (firstVertex <= k) && (k < lastVertex))
269                                                            in->FaceElements,in->Nodes->globalDegreesOfFreedom,              {
270                                                            in->Nodes->globalDegreesOfFreedom);              in->Nodes->globalDegreesOfFreedom[i] = newGlobalDOFID[k - firstVertex];
271                Finley_IndexList_insertElementsWithRowRange(index_list, myFirstVertex, myLastVertex,              setNewDOFId[i] = FALSE;
272                                                            in->ContactElements,in->Nodes->globalDegreesOfFreedom,              }
273                                                            in->Nodes->globalDegreesOfFreedom);          }
274                Finley_IndexList_insertElementsWithRowRange(index_list, myFirstVertex, myLastVertex,  
275                                                            in->Points,in->Nodes->globalDegreesOfFreedom,          if (p < mpiSize - 1)
276                                                            in->Nodes->globalDegreesOfFreedom);          {       /* the final send can be skipped */
277             }  #ifdef ESYS_MPI
278                          MPI_Sendrecv_replace(newGlobalDOFID, len, MPI_INT,
279             /* create the local matrix pattern */                       dest, in->MPIInfo->msg_tag_counter,
280             pattern=Finley_IndexList_createPattern(myNumVertices,index_list,0,globalNumVertices,0);                       source, in->MPIInfo->msg_tag_counter, in->MPIInfo->comm, &status);
281    #endif
282             /* clean up index list */              in->MPIInfo->msg_tag_counter++;
283             if (index_list!=NULL) {              current_rank = Esys_MPIInfo_mod(mpiSize, current_rank - 1);
284                #pragma omp parallel for private(i)          }
285                for(i=0;i<myNumVertices;++i) Finley_IndexList_free(index_list[i].extension);          }
286             }          for (i = 0; i < mpiSize + 1; ++i)
287            distribution[i] = new_distribution[i];
288             if (Finley_noError()) {      }
289        TMPMEMFREE(index_list);
290  /*      }
291        TMPMEMFREE(newGlobalDOFID);
292          ParMETIS_V3_PartGeomKway(distribution,      TMPMEMFREE(setNewDOFId);
293                                   pattern->ptr,      TMPMEMFREE(new_distribution);
294                                   pattern->index,      TMPMEMFREE(partition_count);
295                                   idxtype *vwgt, +      TMPMEMFREE(partition);
296                                   idxtype *adjwgt, +      TMPMEMFREE(xyz);
297                                   int *wgtflag, +      return;
                                  int *numflag, +  
                                  dim,  
                                  xyz,  
                                  int *ncon, +  
                                  mpiSize,  
                                  float *tpwgts, +  
                                  float *ubvec, +  
                                  int *options, +  
                                  int *edgecut, +  
                                  partition,  
                                  in->MPIInfo->comm);  
 */  
                for (i=0;i<myNumVertices;++i) partition[i]=myRank; /* remove */  
            }  
   
            Paso_Pattern_free(pattern);  
             
            /* create a new distributioin and labeling of the DOF */  
            memset(new_distribution,0,mpiSize_size);  
            #pragma omp parallel private(loc_partition_count)  
            {  
                loc_partition_count=THREAD_MEMALLOC(mpiSize,dim_t);  
                memset(loc_partition_count,0,mpiSize_size);  
                #pragma omp for private(i)  
                for (i=0;i<myNumVertices;++i) loc_partition_count[partition[i]]++ ;  
                #pragma omp critical  
                {  
                   for (i=0;i<mpiSize;++i) new_distribution[i]+=loc_partition_count[i];  
                }  
                THREAD_MEMFREE(loc_partition_count);  
            }  
            #ifdef PASO_MPI  
               MPI_Allreduce( new_distribution, partition_count, mpiSize, MPI_INT, MPI_SUM, in->MPIInfo->comm );  
            #else  
                for (i=0;i<mpiSize;++i) partition_count[i]=new_distribution[i];  
            #endif  
            new_distribution[0]=0;  
            for (i=0;i<mpiSize;++i) {  
                new_distribution[i+1]=new_distribution[i]+partition_count[i];  
                partition_count[i]=0;  
            }  
            for (i=0;i<myNumVertices;++i) {  
               rank=partition[i];  
               newGlobalDOFID[i]=new_distribution[rank]+partition_count[rank];  
               partition_count[rank]++;  
            }  
            /* now the overlap needs to be created by sending the partition around*/  
   
            dest=Paso_MPIInfo_mod(mpiSize, myRank + 1);  
            source=Paso_MPIInfo_mod(mpiSize, myRank - 1);  
            current_rank=myRank;  
            #pragma omp parallel for private(i)  
            for (i=0;i<in->Nodes->numNodes;++i) setNewDOFId[i]=TRUE;  
   
            for (p=0; p< mpiSize; ++p) {  
   
                firstVertex=distribution[current_rank];  
                lastVertex=distribution[current_rank+1];  
                #pragma omp parallel for private(i,j,k)  
                for (i=0;i<in->Nodes->numNodes;++i) {  
                    k=in->Nodes->globalDegreesOfFreedom[i];  
                    if (setNewDOFId[i] && (firstVertex<=k) && (k<lastVertex)) {  
                         in->Nodes->globalDegreesOfFreedom[i]=newGlobalDOFID[k-firstVertex];  
                         setNewDOFId[i]=FALSE;  
                    }  
                }  
   
                if (p<mpiSize-1) {  /* the final send can be skipped */  
                   #ifdef PASO_MPI  
                   MPI_Sendrecv_replace(newGlobalDOFID,len, MPI_INT,  
                                        dest, in->MPIInfo->msg_tag_counter,  
                                        source, in->MPIInfo->msg_tag_counter,  
                                        in->MPIInfo->comm,&status);  
                   #endif  
                   in->MPIInfo->msg_tag_counter++;  
                   current_rank=Paso_MPIInfo_mod(mpiSize, current_rank-1);  
               }  
            }  
            for (i=0;i<mpiSize+1;++i) distribution[i]=new_distribution[i];  
   
             
          }  
          TMPMEMFREE(index_list);  
      }  
      TMPMEMFREE(newGlobalDOFID);  
      TMPMEMFREE(setNewDOFId);  
      TMPMEMFREE(new_distribution);  
      TMPMEMFREE(partition_count);  
      TMPMEMFREE(partition);  
      TMPMEMFREE(xyz);  
      return;  
298  }  }

Legend:
Removed from v.1384  
changed lines
  Added in v.3231

  ViewVC Help
Powered by ViewVC 1.1.26