/[escript]/trunk/paso/src/Coupler.c
ViewVC logotype

Annotation of /trunk/paso/src/Coupler.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3458 - (hide annotations)
Mon Jan 31 07:06:42 2011 UTC (8 years, 7 months ago) by gross
File MIME type: text/plain
File size: 11536 byte(s)
more on AMG under MPI
1 ksteube 1313
2     /*******************************************************
3 ksteube 1811 *
4 jfenwick 2881 * Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1811 * 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 ksteube 1313
14 ksteube 1811
15 gross 1552 #include "Coupler.h"
16 jfenwick 3259 #include "esysUtils/error.h"
17 ksteube 1313
18 gross 1552 /*************************************************************
19     *
20     * allocates a Connector
21     *
22     **************************************************************/
23 ksteube 1313
24 gross 1552 Paso_Connector* Paso_Connector_alloc(Paso_SharedComponents* send,
25     Paso_SharedComponents* recv)
26     {
27     Paso_Connector*out=NULL;
28 jfenwick 3259 Esys_resetError();
29 gross 1552 out=MEMALLOC(1,Paso_Connector);
30     if ( send->mpi_info != recv->mpi_info ) {
31 jfenwick 3259 Esys_setError(SYSTEM_ERROR,"Paso_Coupler_alloc: send and recv mpi communicator don't match.");
32 gross 1552 return NULL;
33     }
34 gross 1562 if ( send->local_length != recv->local_length ) {
35 jfenwick 3259 Esys_setError(SYSTEM_ERROR,"Paso_Coupler_alloc: local length of send and recv Paso_SharedComponents must match.");
36 gross 1562 return NULL;
37     }
38    
39 jfenwick 3259 if (!Esys_checkPtr(out)) {
40 gross 1552 out->send=Paso_SharedComponents_getReference(send);
41     out->recv= Paso_SharedComponents_getReference(recv);
42 jfenwick 3259 out->mpi_info = Esys_MPIInfo_getReference(send->mpi_info);
43 gross 1552 out->reference_counter=1;
44 gross 2385
45 gross 1758 /*
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 gross 2385
56 gross 1552 }
57 jfenwick 3259 if (Esys_noError()) {
58 gross 1552 return out;
59     } else {
60     Paso_Connector_free(out);
61     return NULL;
62     }
63     }
64 ksteube 1313
65 gross 1552 /* returns a reference to Connector */
66 ksteube 1313
67 gross 1552 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 ksteube 1313
76 gross 1552 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 jfenwick 3259 Esys_MPIInfo_free(in->mpi_info);
83 gross 1563 MEMFREE(in);
84 gross 1552 #ifdef Paso_TRACE
85     printf("Paso_Coupler_dealloc: system matrix pattern as been deallocated.\n");
86     #endif
87     }
88     }
89     }
90 ksteube 1313
91 gross 1552 Paso_Connector* Paso_Connector_copy(Paso_Connector* in) {
92     return Paso_Connector_unroll(in,1);
93     }
94 ksteube 1313
95 gross 1552 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 jfenwick 3259 if (Esys_noError()) {
99 gross 1552 if (block_size>1) {
100 gross 1562 new_send_shcomp=Paso_SharedComponents_alloc(in->send->local_length,
101     in->send->numNeighbors,
102 gross 1552 in->send->neighbor,
103     in->send->shared,
104     in->send->offsetInShared,
105     block_size,0,in->mpi_info);
106 ksteube 1313
107 gross 1562 new_recv_shcomp=Paso_SharedComponents_alloc(in->recv->local_length,
108     in->recv->numNeighbors,
109 gross 1552 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 jfenwick 3259 if (Esys_noError()) out=Paso_Connector_alloc(new_send_shcomp,new_recv_shcomp);
118 gross 1552 }
119     Paso_SharedComponents_free(new_send_shcomp);
120     Paso_SharedComponents_free(new_recv_shcomp);
121 jfenwick 3259 if (Esys_noError()) {
122 gross 1552 return out;
123     } else {
124     Paso_Connector_free(out);
125     return NULL;
126     }
127     }
128     /*************************************************************
129     *
130     * allocates a Connector
131     *
132     **************************************************************/
133 ksteube 1313
134 gross 1552 Paso_Coupler* Paso_Coupler_alloc(Paso_Connector* connector, dim_t block_size)
135 ksteube 1313 {
136 jfenwick 3259 Esys_MPIInfo *mpi_info = connector->mpi_info;
137 ksteube 1313 Paso_Coupler*out=NULL;
138 jfenwick 3259 Esys_resetError();
139 ksteube 1313 out=MEMALLOC(1,Paso_Coupler);
140 jfenwick 3259 if (!Esys_checkPtr(out)) {
141 gross 1562 out->data=NULL;
142 gross 1552 out->block_size=block_size;
143     out->connector=Paso_Connector_getReference(connector);
144 ksteube 1313 out->send_buffer=NULL;
145     out->recv_buffer=NULL;
146     out->mpi_requests=NULL;
147     out->mpi_stati=NULL;
148 jfenwick 3259 out->mpi_info = Esys_MPIInfo_getReference(mpi_info);
149 ksteube 1313 out->reference_counter=1;
150 gross 3458 out->in_use = FALSE;
151 ksteube 1313
152 jfenwick 3259 #ifdef ESYS_MPI
153 gross 1552 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 jfenwick 3259 Esys_checkPtr(out->mpi_requests);
156     Esys_checkPtr(out->mpi_stati);
157 ksteube 1313 #endif
158 gross 1552 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 jfenwick 3259 Esys_checkPtr(out->send_buffer);
162     Esys_checkPtr(out->recv_buffer);
163 gross 1552 }
164 ksteube 1313 }
165 jfenwick 3259 if (Esys_noError()) {
166 ksteube 1313 return out;
167     } else {
168     Paso_Coupler_free(out);
169     return NULL;
170     }
171     }
172    
173 gross 1552 /* returns a reference to Coupler */
174 ksteube 1313
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 gross 1552 Paso_Connector_free(in->connector);
189 ksteube 1313 MEMFREE(in->send_buffer);
190     MEMFREE(in->recv_buffer);
191     MEMFREE(in->mpi_requests);
192     MEMFREE(in->mpi_stati);
193 jfenwick 3259 Esys_MPIInfo_free(in->mpi_info);
194 ksteube 1313 MEMFREE(in);
195     }
196     }
197     }
198    
199    
200 gross 1639 void Paso_Coupler_startCollect(Paso_Coupler* coupler,const double* in)
201 ksteube 1313 {
202 jfenwick 3259 Esys_MPIInfo *mpi_info = coupler->mpi_info;
203 ksteube 1313 dim_t block_size=coupler->block_size;
204     size_t block_size_size=block_size*sizeof(double);
205 phornby 1628 dim_t i;
206 gross 1639 coupler->data=(double*) in;
207 ksteube 1313 if ( mpi_info->size>1) {
208 gross 3458 if (coupler->in_use) {
209     Esys_setError(SYSTEM_ERROR,"Paso_Coupler_startCollect: Coupler in use.");
210     }
211 ksteube 1313 /* start reveiving input */
212     {
213 gross 1552 for (i=0; i< coupler->connector->recv->numNeighbors; ++i) {
214 jfenwick 3259 #ifdef ESYS_MPI
215 gross 1552 MPI_Irecv(&(coupler->recv_buffer[coupler->connector->recv->offsetInShared[i] * block_size]),
216     (coupler->connector->recv->offsetInShared[i+1]- coupler->connector->recv->offsetInShared[i])*block_size,
217 ksteube 1313 MPI_DOUBLE,
218 gross 1552 coupler->connector->recv->neighbor[i],
219 gross 1553 mpi_info->msg_tag_counter+coupler->connector->recv->neighbor[i],
220 ksteube 1313 mpi_info->comm,
221     &(coupler->mpi_requests[i]));
222     #endif
223    
224     }
225     }
226     /* collect values into buffer */
227 gross 1758 if (block_size>1) {
228     #pragma omp parallel for private(i)
229     for (i=0; i < coupler->connector->send->numSharedComponents;++i) {
230     memcpy(&(coupler->send_buffer[(block_size)*i]),&(in[ block_size * coupler->connector->send->shared[i]]), block_size_size);
231     }
232     } else {
233     #pragma omp parallel for private(i)
234     for (i=0; i < coupler->connector->send->numSharedComponents;++i) coupler->send_buffer[i]=in[coupler->connector->send->shared[i]];
235 ksteube 1313 }
236     /* send buffer out */
237     {
238 gross 1552 for (i=0; i< coupler->connector->send->numNeighbors; ++i) {
239 jfenwick 3259 #ifdef ESYS_MPI
240 gross 1553 MPI_Issend(&(coupler->send_buffer[coupler->connector->send->offsetInShared[i] * block_size]),
241 gross 1552 (coupler->connector->send->offsetInShared[i+1]- coupler->connector->send->offsetInShared[i])*block_size,
242 ksteube 1313 MPI_DOUBLE,
243 gross 1552 coupler->connector->send->neighbor[i],
244 ksteube 1313 mpi_info->msg_tag_counter+mpi_info->rank,
245     mpi_info->comm,
246 gross 1553 &(coupler->mpi_requests[i+ coupler->connector->recv->numNeighbors]));
247 ksteube 1313 #endif
248     }
249     }
250 gross 1413 mpi_info->msg_tag_counter+=mpi_info->size;
251 gross 3458 coupler->in_use=TRUE;
252 ksteube 1313 }
253     }
254    
255     double* Paso_Coupler_finishCollect(Paso_Coupler* coupler)
256     {
257 jfenwick 3259 Esys_MPIInfo *mpi_info = coupler->mpi_info;
258 ksteube 1313 if ( mpi_info->size>1) {
259 gross 3458 if (! coupler->in_use) {
260     Esys_setError(SYSTEM_ERROR,"Paso_Coupler_finishCollect: Communication has not been initiated.");
261     return NULL;
262     }
263 ksteube 1313 /* wait for receive */
264 gross 3458 #ifdef ESYS_MPI
265 gross 1552 MPI_Waitall(coupler->connector->recv->numNeighbors+coupler->connector->send->numNeighbors,
266 ksteube 1313 coupler->mpi_requests,
267     coupler->mpi_stati);
268 gross 3458 #endif
269     coupler->in_use=FALSE;
270 ksteube 1313 }
271 gross 1758
272 ksteube 1313 return coupler->recv_buffer;
273     }
274 gross 3442
275 gross 1562 void Paso_Coupler_copyAll(const Paso_Coupler* src, Paso_Coupler* target)
276     {
277     dim_t i;
278     #pragma omp parallel
279     {
280     #pragma omp for private(i)
281     for (i =0; i< src->connector->recv->numSharedComponents * src->block_size; ++i) {
282     target->recv_buffer[i]=src->recv_buffer[i];
283     }
284     #pragma omp for private(i)
285     for (i =0; i< Paso_Coupler_getLocalLength(src) * src->block_size; ++i) {
286     target->data[i]=src->data[i];
287     }
288     }
289     }
290 gross 3442
291 gross 3445 /* aadds a local vector x * a to the vector y including overlap :*/
292     void Paso_Coupler_add(const dim_t n, double* x, const double a, const double* y, Paso_Coupler *coupler)
293 gross 3442 {
294     double *remote_values = NULL;
295     const dim_t overlap_n = Paso_Coupler_getNumOverlapValues(coupler) ;
296 gross 3445 const dim_t my_n= n - overlap_n;
297 gross 3442 dim_t i;
298    
299     if (ABS(a) > 0 ) {
300     Paso_Coupler_startCollect(coupler, x);
301    
302     if ( a == 1.) {
303     #pragma omp parallel for private(i)
304     for (i=0; i<my_n; i++) {
305     x[i]+=y[i];
306     }
307     } else {
308     #pragma omp parallel for private(i)
309     for (i=0; i<my_n; i++) {
310     x[i]+=a*y[i];
311     }
312     }
313    
314     Paso_Coupler_finishCollect(coupler);
315     remote_values=coupler->recv_buffer;
316    
317     if ( a == 1.) {
318     #pragma omp parallel for private(i)
319     for (i=0;i<overlap_n; ++i) {
320 gross 3445 x[my_n+i]+=remote_values[i];
321 gross 3442 }
322     } else {
323     #pragma omp parallel for private(i)
324     for (i=0;i<overlap_n; ++i) {
325 gross 3445 x[my_n+i]+=a*remote_values[i];
326 gross 3442 }
327     }
328     }
329 gross 3443 }
330 gross 3445
331     /* adjusts max values across shared values x */
332     void Paso_Coupler_max(const dim_t n, double* x, Paso_Coupler *coupler)
333     {
334     double *remote_values = NULL;
335     const dim_t overlap_n = Paso_Coupler_getNumOverlapValues(coupler) ;
336     const dim_t my_n= n - overlap_n;
337     dim_t i;
338    
339     Paso_Coupler_startCollect(coupler, x);
340     Paso_Coupler_finishCollect(coupler);
341     remote_values=coupler->recv_buffer;
342     #pragma omp parallel for private(i)
343     for (i=0;i<overlap_n; ++i) x[my_n+i]=MAX(x[my_n+i], remote_values[i]);
344     }

  ViewVC Help
Powered by ViewVC 1.1.26