/[escript]/trunk-mpi-branch/finley/src/Mesh_optimizeDOFDistribution.c
ViewVC logotype

Contents of /trunk-mpi-branch/finley/src/Mesh_optimizeDOFDistribution.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1301 - (show annotations)
Mon Sep 17 04:13:48 2007 UTC (12 years, 8 months ago) by ksteube
File MIME type: text/plain
File size: 9764 byte(s)
Lumping was not working because a loop over colors was missing.
Also resolved a few more compiler warnings that turned up.

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

  ViewVC Help
Powered by ViewVC 1.1.26