/[escript]/branches/amg_from_3530/paso/src/Coupler.c
ViewVC logotype

Contents of /branches/amg_from_3530/paso/src/Coupler.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3767 - (show annotations)
Fri Jan 13 01:52:46 2012 UTC (7 years, 5 months ago) by lgao
File MIME type: text/plain
File size: 11721 byte(s)
A version which passed all my small test cases. Ready for the test cases 
in Unittest. 


1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14
15 #include "Coupler.h"
16 #include "esysUtils/error.h"
17
18 /*************************************************************
19 *
20 * allocates a Connector
21 *
22 **************************************************************/
23
24 Paso_Connector* Paso_Connector_alloc(Paso_SharedComponents* send,
25 Paso_SharedComponents* recv)
26 {
27 Paso_Connector*out=NULL;
28 Esys_resetError();
29 out=MEMALLOC(1,Paso_Connector);
30 if ( send->mpi_info != recv->mpi_info ) {
31 Esys_setError(SYSTEM_ERROR,"Paso_Coupler_alloc: send and recv mpi communicator don't match.");
32 return NULL;
33 }
34 if ( send->local_length != recv->local_length ) {
35 Esys_setError(SYSTEM_ERROR,"Paso_Coupler_alloc: local length of send and recv Paso_SharedComponents must match.");
36 return NULL;
37 }
38
39 if (!Esys_checkPtr(out)) {
40 out->send=Paso_SharedComponents_getReference(send);
41 out->recv= Paso_SharedComponents_getReference(recv);
42 out->mpi_info = Esys_MPIInfo_getReference(send->mpi_info);
43 out->reference_counter=1;
44
45 /*
46 { int i;
47 for (i=0; i< out->recv->numNeighbors; ++i)
48 printf("Coupler: %d receive %d data at %d from %d\n",send->mpi_info->rank,out->recv->offsetInShared[i+1]- out->recv->offsetInShared[i],
49 out->recv->offsetInShared[i],out->recv->neighbor[i]);
50 for (i=0; i< out->send->numNeighbors; ++i)
51 printf("Coupler: %d send %d data at %d to %d\n",send->mpi_info->rank,out->send->offsetInShared[i+1]- out->send->offsetInShared[i],
52 out->send->offsetInShared[i],out->send->neighbor[i]);
53 }
54 */
55
56 }
57 if (Esys_noError()) {
58 return out;
59 } else {
60 Paso_Connector_free(out);
61 return NULL;
62 }
63 }
64
65 /* returns a reference to Connector */
66
67 Paso_Connector* Paso_Connector_getReference(Paso_Connector* in) {
68 if (in!=NULL) {
69 ++(in->reference_counter);
70 }
71 return in;
72 }
73
74 /* deallocates a Connector: */
75
76 void Paso_Connector_free(Paso_Connector* in) {
77 if (in!=NULL) {
78 in->reference_counter--;
79 if (in->reference_counter<=0) {
80 Paso_SharedComponents_free(in->send);
81 Paso_SharedComponents_free(in->recv);
82 Esys_MPIInfo_free(in->mpi_info);
83 MEMFREE(in);
84 #ifdef Paso_TRACE
85 printf("Paso_Coupler_dealloc: system matrix pattern as been deallocated.\n");
86 #endif
87 }
88 }
89 }
90
91 Paso_Connector* Paso_Connector_copy(Paso_Connector* in) {
92 return Paso_Connector_unroll(in,1);
93 }
94
95 Paso_Connector* Paso_Connector_unroll(Paso_Connector* in, index_t block_size) {
96 Paso_SharedComponents *new_send_shcomp=NULL, *new_recv_shcomp=NULL;
97 Paso_Connector *out=NULL;
98 if (Esys_noError()) {
99 if (block_size>1) {
100 new_send_shcomp=Paso_SharedComponents_alloc(in->send->local_length,
101 in->send->numNeighbors,
102 in->send->neighbor,
103 in->send->shared,
104 in->send->offsetInShared,
105 block_size,0,in->mpi_info);
106
107 new_recv_shcomp=Paso_SharedComponents_alloc(in->recv->local_length,
108 in->recv->numNeighbors,
109 in->recv->neighbor,
110 in->recv->shared,
111 in->recv->offsetInShared,
112 block_size,0,in->mpi_info);
113 } else {
114 new_send_shcomp=Paso_SharedComponents_getReference(in->send);
115 new_recv_shcomp=Paso_SharedComponents_getReference(in->recv);
116 }
117 if (Esys_noError()) out=Paso_Connector_alloc(new_send_shcomp,new_recv_shcomp);
118 }
119 Paso_SharedComponents_free(new_send_shcomp);
120 Paso_SharedComponents_free(new_recv_shcomp);
121 if (Esys_noError()) {
122 return out;
123 } else {
124 Paso_Connector_free(out);
125 return NULL;
126 }
127 }
128 /*************************************************************
129 *
130 * allocates a Connector
131 *
132 **************************************************************/
133
134 Paso_Coupler* Paso_Coupler_alloc(Paso_Connector* connector, dim_t block_size)
135 {
136 Esys_MPIInfo *mpi_info = connector->mpi_info;
137 Paso_Coupler*out=NULL;
138 Esys_resetError();
139 out=MEMALLOC(1,Paso_Coupler);
140 if (!Esys_checkPtr(out)) {
141 out->data=NULL;
142 out->block_size=block_size;
143 out->connector=Paso_Connector_getReference(connector);
144 out->send_buffer=NULL;
145 out->recv_buffer=NULL;
146 out->mpi_requests=NULL;
147 out->mpi_stati=NULL;
148 out->mpi_info = Esys_MPIInfo_getReference(mpi_info);
149 out->reference_counter=1;
150 out->in_use = FALSE;
151
152 #ifdef ESYS_MPI
153 out->mpi_requests=MEMALLOC(connector->send->numNeighbors+connector->recv->numNeighbors,MPI_Request);
154 out->mpi_stati=MEMALLOC(connector->send->numNeighbors+connector->recv->numNeighbors,MPI_Status);
155 Esys_checkPtr(out->mpi_requests);
156 Esys_checkPtr(out->mpi_stati);
157 #endif
158 if (mpi_info->size>1) {
159 out->send_buffer=MEMALLOC(connector->send->numSharedComponents * block_size,double);
160 out->recv_buffer=MEMALLOC(connector->recv->numSharedComponents * block_size,double);
161 Esys_checkPtr(out->send_buffer);
162 Esys_checkPtr(out->recv_buffer);
163 }
164 }
165 if (Esys_noError()) {
166 return out;
167 } else {
168 Paso_Coupler_free(out);
169 return NULL;
170 }
171 }
172
173 /* returns a reference to Coupler */
174
175 Paso_Coupler* Paso_Coupler_getReference(Paso_Coupler* in) {
176 if (in!=NULL) {
177 ++(in->reference_counter);
178 }
179 return in;
180 }
181
182 /* deallocates a Coupler: */
183
184 void Paso_Coupler_free(Paso_Coupler* in) {
185 if (in!=NULL) {
186 in->reference_counter--;
187 if (in->reference_counter<=0) {
188 Paso_Connector_free(in->connector);
189 MEMFREE(in->send_buffer);
190 MEMFREE(in->recv_buffer);
191 MEMFREE(in->mpi_requests);
192 MEMFREE(in->mpi_stati);
193 Esys_MPIInfo_free(in->mpi_info);
194 MEMFREE(in);
195 }
196 }
197 }
198
199
200 void Paso_Coupler_startCollect(Paso_Coupler* coupler,const double* in)
201 {
202 Esys_MPIInfo *mpi_info = coupler->mpi_info;
203 dim_t block_size=coupler->block_size;
204 size_t block_size_size=block_size*sizeof(double);
205 dim_t i;
206 coupler->data=(double*) in;
207 if ( mpi_info->size>1) {
208 if (coupler->in_use) {
209 Esys_setError(SYSTEM_ERROR,"Paso_Coupler_startCollect: Coupler in use.");
210 }
211 /* start reveiving input */
212 {
213 for (i=0; i< coupler->connector->recv->numNeighbors; ++i) {
214 /*fprintf(stderr, "rank %d COLLECT RECV %d from %d tag %d (i%d) \n", mpi_info->rank,
215 (coupler->connector->recv->offsetInShared[i+1]- coupler->connector->recv->offsetInShared[i])*block_size,
216 coupler->connector->recv->neighbor[i],
217 mpi_info->msg_tag_counter+coupler->connector->recv->neighbor[i], i);*/
218 #ifdef ESYS_MPI
219 MPI_Irecv(&(coupler->recv_buffer[coupler->connector->recv->offsetInShared[i] * block_size]),
220 (coupler->connector->recv->offsetInShared[i+1]- coupler->connector->recv->offsetInShared[i])*block_size,
221 MPI_DOUBLE,
222 coupler->connector->recv->neighbor[i],
223 mpi_info->msg_tag_counter+coupler->connector->recv->neighbor[i],
224 mpi_info->comm,
225 &(coupler->mpi_requests[i]));
226 #endif
227 }
228 }
229 /* collect values into buffer */
230 if (block_size>1) {
231 #pragma omp parallel for private(i)
232 for (i=0; i < coupler->connector->send->numSharedComponents;++i) {
233 memcpy(&(coupler->send_buffer[(block_size)*i]),&(in[ block_size * coupler->connector->send->shared[i]]), block_size_size);
234 }
235 } else {
236 #pragma omp parallel for private(i)
237 for (i=0; i < coupler->connector->send->numSharedComponents;++i) coupler->send_buffer[i]=in[coupler->connector->send->shared[i]];
238 }
239 /* send buffer out */
240 {
241
242 for (i=0; i< coupler->connector->send->numNeighbors; ++i) {
243 /*fprintf(stderr, "rank %d COLLECT SEND %d to %d tag %d (i%d)\n", mpi_info->rank,
244 (coupler->connector->send->offsetInShared[i+1]- coupler->connector->send->offsetInShared[i])*block_size,
245 coupler->connector->send->neighbor[i],
246 mpi_info->msg_tag_counter+mpi_info->rank, i);*/
247 #ifdef ESYS_MPI
248 MPI_Issend(&(coupler->send_buffer[coupler->connector->send->offsetInShared[i] * block_size]),
249 (coupler->connector->send->offsetInShared[i+1]- coupler->connector->send->offsetInShared[i])*block_size,
250 MPI_DOUBLE,
251 coupler->connector->send->neighbor[i],
252 mpi_info->msg_tag_counter+mpi_info->rank,
253 mpi_info->comm,
254 &(coupler->mpi_requests[i+ coupler->connector->recv->numNeighbors]));
255 #endif
256 }
257 }
258 mpi_info->msg_tag_counter+=mpi_info->size;
259 coupler->in_use=TRUE;
260 }
261 }
262
263 double* Paso_Coupler_finishCollect(Paso_Coupler* coupler)
264 {
265 Esys_MPIInfo *mpi_info = coupler->mpi_info;
266 if ( mpi_info->size>1) {
267 if (! coupler->in_use) {
268 Esys_setError(SYSTEM_ERROR,"Paso_Coupler_finishCollect: Communication has not been initiated.");
269 return NULL;
270 }
271 /* wait for receive */
272 #ifdef ESYS_MPI
273 MPI_Waitall(coupler->connector->recv->numNeighbors+coupler->connector->send->numNeighbors,
274 coupler->mpi_requests,
275 coupler->mpi_stati);
276 #endif
277 coupler->in_use=FALSE;
278 }
279
280 return coupler->recv_buffer;
281 }
282
283 void Paso_Coupler_copyAll(const Paso_Coupler* src, Paso_Coupler* target)
284 {
285 dim_t i;
286 #pragma omp parallel
287 {
288 #pragma omp for private(i)
289 for (i =0; i< src->connector->recv->numSharedComponents * src->block_size; ++i) {
290 target->recv_buffer[i]=src->recv_buffer[i];
291 }
292 #pragma omp for private(i)
293 for (i =0; i< Paso_Coupler_getLocalLength(src) * src->block_size; ++i) {
294 target->data[i]=src->data[i];
295 }
296 }
297 }
298
299 /* */
300 void Paso_Coupler_fillOverlap(const dim_t n, double* x, Paso_Coupler *coupler)
301 {
302 double *remote_values = NULL;
303 const dim_t overlap_n = Paso_Coupler_getNumOverlapValues(coupler) ;
304 const dim_t my_n= n - overlap_n;
305 const dim_t block_size = coupler->block_size;
306 const dim_t offset = block_size * my_n;
307 dim_t i;
308
309 Paso_Coupler_startCollect(coupler, x);
310 Paso_Coupler_finishCollect(coupler);
311 remote_values=coupler->recv_buffer;
312
313 #pragma omp parallel for private(i)
314 for (i=0;i<overlap_n * block_size; ++i) {
315 x[offset+i]=remote_values[i];
316 }
317 }
318
319 /* adjusts max values across shared values x */
320 void Paso_Coupler_max(const dim_t n, double* x, Paso_Coupler *coupler)
321 {
322 double *remote_values = NULL;
323 const dim_t overlap_n = Paso_Coupler_getNumOverlapValues(coupler) ;
324 const dim_t my_n= n - overlap_n;
325 dim_t i;
326
327 Paso_Coupler_startCollect(coupler, x);
328 Paso_Coupler_finishCollect(coupler);
329 remote_values=coupler->recv_buffer;
330 #pragma omp parallel for private(i)
331 for (i=0;i<overlap_n; ++i) x[my_n+i]=MAX(x[my_n+i], remote_values[i]);
332 }

  ViewVC Help
Powered by ViewVC 1.1.26