/[escript]/temp/finley/src/Mesh_optimizeDOFDistribution.c
ViewVC logotype

Contents of /temp/finley/src/Mesh_optimizeDOFDistribution.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1387 - (show annotations)
Fri Jan 11 07:45:26 2008 UTC (11 years, 8 months ago) by trankine
File MIME type: text/plain
File size: 9514 byte(s)
Restore the trunk that existed before the windows changes were committed to the (now moved to branches) old trunk.
1
2 /* $Id$ */
3
4 /*******************************************************
5 *
6 * Copyright 2003-2007 by ACceSS MNRF
7 * Copyright 2007 by University of Queensland
8 *
9 * http://esscc.uq.edu.au
10 * Primary Business: Queensland, Australia
11 * Licensed under the Open Software License version 3.0
12 * http://www.opensource.org/licenses/osl-3.0.php
13 *
14 *******************************************************/
15
16 /**************************************************************/
17
18 /* Finley: Mesh: optimizes the distribution of DOFs across processors */
19 /* using ParMETIS. On return a new distribution is given and the globalDOF are relabled */
20 /* accordingly but the mesh has not been redesitributed yet */
21
22 /**************************************************************/
23
24 #include "Mesh.h"
25 #include "IndexList.h"
26 #ifdef _OPENMP
27 #include <omp.h>
28 #endif
29
30 /**************************************************************/
31
32 void Finley_Mesh_optimizeDOFDistribution(Finley_Mesh* in,dim_t *distribution) {
33
34 dim_t dim, i,j,k, myNumVertices,p, mpiSize, len, globalNumVertices,*partition_count=NULL, *new_distribution=NULL, *loc_partition_count=NULL;
35 bool_t *setNewDOFId=NULL;
36 index_t myFirstVertex, myLastVertex, firstVertex, lastVertex, *newGlobalDOFID=NULL;
37 size_t mpiSize_size;
38 index_t* partition=NULL;
39 Paso_Pattern *pattern=NULL;
40 Paso_MPI_rank myRank,dest,source,current_rank, rank;
41 Finley_IndexList* index_list=NULL;
42 float *xyz=NULL;
43
44 #ifdef PASO_MPI
45 MPI_Status status;
46 #endif
47
48 if (in==NULL) return;
49 if (in->Nodes == NULL) return;
50
51 myRank=in->MPIInfo->rank;
52 mpiSize=in->MPIInfo->size;
53 mpiSize_size=mpiSize*sizeof(dim_t);
54 dim=in->Nodes->numDim;
55 /* first step is to distribute the elements according to a global X of DOF */
56
57 myFirstVertex=distribution[myRank];
58 myLastVertex=distribution[myRank+1];
59 myNumVertices=myLastVertex-myFirstVertex;
60 globalNumVertices=distribution[mpiSize];
61 len=0;
62 for (p=0;p<mpiSize;++p) len=MAX(len,distribution[p+1]-distribution[p]);
63 partition=TMPMEMALLOC(len,index_t); /* len is used for the sending around of partition later on */
64 xyz=TMPMEMALLOC(myNumVertices*dim,float);
65 partition_count=TMPMEMALLOC(mpiSize+1,dim_t);
66 new_distribution=TMPMEMALLOC(mpiSize+1,dim_t);
67 newGlobalDOFID=TMPMEMALLOC(len,index_t);
68 setNewDOFId=TMPMEMALLOC(in->Nodes->numNodes,bool_t);
69 if (!(Finley_checkPtr(partition) || Finley_checkPtr(xyz) || Finley_checkPtr(partition_count) || Finley_checkPtr(partition_count) || Finley_checkPtr(newGlobalDOFID) || Finley_checkPtr(setNewDOFId))) {
70
71 /* set the coordinates: *?
72 /* it is assumed that at least one node on this processor provides a coordinate */
73 #pragma omp parallel for private(i,j,k)
74 for (i=0;i<in->Nodes->numNodes;++i) {
75 k=in->Nodes->globalDegreesOfFreedom[i]-myFirstVertex;
76 if ((k>=0) && (k<myNumVertices)) {
77 for (j=0;j<dim;++j) xyz[k*dim+j]=(float)(in->Nodes->Coordinates[INDEX2(j,i,dim)]);
78 }
79 }
80
81 index_list=TMPMEMALLOC(myNumVertices,Finley_IndexList);
82 /* ksteube CSR of DOF IDs */
83 /* create the adjacency structure xadj and adjncy */
84 if (! Finley_checkPtr(index_list)) {
85 #pragma omp parallel private(i)
86 {
87 #pragma omp for schedule(static)
88 for(i=0;i<myNumVertices;++i) {
89 index_list[i].extension=NULL;
90 index_list[i].n=0;
91 }
92 /* ksteube build CSR format */
93 /* insert contributions from element matrices into colums index index_list: */
94 Finley_IndexList_insertElementsWithRowRange(index_list, myFirstVertex, myLastVertex,
95 in->Elements,in->Nodes->globalDegreesOfFreedom,
96 in->Nodes->globalDegreesOfFreedom);
97 Finley_IndexList_insertElementsWithRowRange(index_list, myFirstVertex, myLastVertex,
98 in->FaceElements,in->Nodes->globalDegreesOfFreedom,
99 in->Nodes->globalDegreesOfFreedom);
100 Finley_IndexList_insertElementsWithRowRange(index_list, myFirstVertex, myLastVertex,
101 in->ContactElements,in->Nodes->globalDegreesOfFreedom,
102 in->Nodes->globalDegreesOfFreedom);
103 Finley_IndexList_insertElementsWithRowRange(index_list, myFirstVertex, myLastVertex,
104 in->Points,in->Nodes->globalDegreesOfFreedom,
105 in->Nodes->globalDegreesOfFreedom);
106 }
107
108 /* create the local matrix pattern */
109 pattern=Finley_IndexList_createPattern(myNumVertices,index_list,0,globalNumVertices,0);
110
111 /* clean up index list */
112 if (index_list!=NULL) {
113 #pragma omp parallel for private(i)
114 for(i=0;i<myNumVertices;++i) Finley_IndexList_free(index_list[i].extension);
115 }
116
117 if (Finley_noError()) {
118
119 /*
120
121 ParMETIS_V3_PartGeomKway(distribution,
122 pattern->ptr,
123 pattern->index,
124 idxtype *vwgt, +
125 idxtype *adjwgt, +
126 int *wgtflag, +
127 int *numflag, +
128 dim,
129 xyz,
130 int *ncon, +
131 mpiSize,
132 float *tpwgts, +
133 float *ubvec, +
134 int *options, +
135 int *edgecut, +
136 partition,
137 in->MPIInfo->comm);
138 */
139 for (i=0;i<myNumVertices;++i) partition[i]=myRank; /* remove */
140 }
141
142 Paso_Pattern_free(pattern);
143
144 /* create a new distributioin and labeling of the DOF */
145 memset(new_distribution,0,mpiSize_size);
146 #pragma omp parallel private(loc_partition_count)
147 {
148 loc_partition_count=THREAD_MEMALLOC(mpiSize,dim_t);
149 memset(loc_partition_count,0,mpiSize_size);
150 #pragma omp for private(i)
151 for (i=0;i<myNumVertices;++i) loc_partition_count[partition[i]]++ ;
152 #pragma omp critical
153 {
154 for (i=0;i<mpiSize;++i) new_distribution[i]+=loc_partition_count[i];
155 }
156 THREAD_MEMFREE(loc_partition_count);
157 }
158 #ifdef PASO_MPI
159 MPI_Allreduce( new_distribution, partition_count, mpiSize, MPI_INT, MPI_SUM, in->MPIInfo->comm );
160 #else
161 for (i=0;i<mpiSize;++i) partition_count[i]=new_distribution[i];
162 #endif
163 new_distribution[0]=0;
164 for (i=0;i<mpiSize;++i) {
165 new_distribution[i+1]=new_distribution[i]+partition_count[i];
166 partition_count[i]=0;
167 }
168 for (i=0;i<myNumVertices;++i) {
169 rank=partition[i];
170 newGlobalDOFID[i]=new_distribution[rank]+partition_count[rank];
171 partition_count[rank]++;
172 }
173 /* now the overlap needs to be created by sending the partition around*/
174
175 dest=Paso_MPIInfo_mod(mpiSize, myRank + 1);
176 source=Paso_MPIInfo_mod(mpiSize, myRank - 1);
177 current_rank=myRank;
178 #pragma omp parallel for private(i)
179 for (i=0;i<in->Nodes->numNodes;++i) setNewDOFId[i]=TRUE;
180
181 for (p=0; p< mpiSize; ++p) {
182
183 firstVertex=distribution[current_rank];
184 lastVertex=distribution[current_rank+1];
185 #pragma omp parallel for private(i,j,k)
186 for (i=0;i<in->Nodes->numNodes;++i) {
187 k=in->Nodes->globalDegreesOfFreedom[i];
188 if (setNewDOFId[i] && (firstVertex<=k) && (k<lastVertex)) {
189 in->Nodes->globalDegreesOfFreedom[i]=newGlobalDOFID[k-firstVertex];
190 setNewDOFId[i]=FALSE;
191 }
192 }
193
194 if (p<mpiSize-1) { /* the final send can be skipped */
195 #ifdef PASO_MPI
196 MPI_Sendrecv_replace(newGlobalDOFID,len, MPI_INT,
197 dest, in->MPIInfo->msg_tag_counter,
198 source, in->MPIInfo->msg_tag_counter,
199 in->MPIInfo->comm,&status);
200 #endif
201 in->MPIInfo->msg_tag_counter++;
202 current_rank=Paso_MPIInfo_mod(mpiSize, current_rank-1);
203 }
204 }
205 for (i=0;i<mpiSize+1;++i) distribution[i]=new_distribution[i];
206
207
208 }
209 TMPMEMFREE(index_list);
210 }
211 TMPMEMFREE(newGlobalDOFID);
212 TMPMEMFREE(setNewDOFId);
213 TMPMEMFREE(new_distribution);
214 TMPMEMFREE(partition_count);
215 TMPMEMFREE(partition);
216 TMPMEMFREE(xyz);
217 return;
218 }

  ViewVC Help
Powered by ViewVC 1.1.26