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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3445 - (hide annotations)
Wed Jan 19 06:02:15 2011 UTC (8 years, 8 months ago) by gross
File MIME type: text/plain
File size: 11199 byte(s)
more work towards MPI AMG
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    
151 jfenwick 3259 #ifdef ESYS_MPI
152 gross 1552 out->mpi_requests=MEMALLOC(connector->send->numNeighbors+connector->recv->numNeighbors,MPI_Request);
153     out->mpi_stati=MEMALLOC(connector->send->numNeighbors+connector->recv->numNeighbors,MPI_Status);
154 jfenwick 3259 Esys_checkPtr(out->mpi_requests);
155     Esys_checkPtr(out->mpi_stati);
156 ksteube 1313 #endif
157 gross 1552 if (mpi_info->size>1) {
158     out->send_buffer=MEMALLOC(connector->send->numSharedComponents * block_size,double);
159     out->recv_buffer=MEMALLOC(connector->recv->numSharedComponents * block_size,double);
160 jfenwick 3259 Esys_checkPtr(out->send_buffer);
161     Esys_checkPtr(out->recv_buffer);
162 gross 1552 }
163 ksteube 1313 }
164 jfenwick 3259 if (Esys_noError()) {
165 ksteube 1313 return out;
166     } else {
167     Paso_Coupler_free(out);
168     return NULL;
169     }
170     }
171    
172 gross 1552 /* returns a reference to Coupler */
173 ksteube 1313
174     Paso_Coupler* Paso_Coupler_getReference(Paso_Coupler* in) {
175     if (in!=NULL) {
176     ++(in->reference_counter);
177     }
178     return in;
179     }
180    
181     /* deallocates a Coupler: */
182    
183     void Paso_Coupler_free(Paso_Coupler* in) {
184     if (in!=NULL) {
185     in->reference_counter--;
186     if (in->reference_counter<=0) {
187 gross 1552 Paso_Connector_free(in->connector);
188 ksteube 1313 MEMFREE(in->send_buffer);
189     MEMFREE(in->recv_buffer);
190     MEMFREE(in->mpi_requests);
191     MEMFREE(in->mpi_stati);
192 jfenwick 3259 Esys_MPIInfo_free(in->mpi_info);
193 ksteube 1313 MEMFREE(in);
194     }
195     }
196     }
197    
198    
199 gross 1639 void Paso_Coupler_startCollect(Paso_Coupler* coupler,const double* in)
200 ksteube 1313 {
201 jfenwick 3259 Esys_MPIInfo *mpi_info = coupler->mpi_info;
202 ksteube 1313 dim_t block_size=coupler->block_size;
203     size_t block_size_size=block_size*sizeof(double);
204 phornby 1628 dim_t i;
205 gross 1639 coupler->data=(double*) in;
206 ksteube 1313 if ( mpi_info->size>1) {
207     /* start reveiving input */
208     {
209 gross 1552 for (i=0; i< coupler->connector->recv->numNeighbors; ++i) {
210 jfenwick 3259 #ifdef ESYS_MPI
211 gross 1552 MPI_Irecv(&(coupler->recv_buffer[coupler->connector->recv->offsetInShared[i] * block_size]),
212     (coupler->connector->recv->offsetInShared[i+1]- coupler->connector->recv->offsetInShared[i])*block_size,
213 ksteube 1313 MPI_DOUBLE,
214 gross 1552 coupler->connector->recv->neighbor[i],
215 gross 1553 mpi_info->msg_tag_counter+coupler->connector->recv->neighbor[i],
216 ksteube 1313 mpi_info->comm,
217     &(coupler->mpi_requests[i]));
218     #endif
219    
220     }
221     }
222     /* collect values into buffer */
223 gross 1758 if (block_size>1) {
224     #pragma omp parallel for private(i)
225     for (i=0; i < coupler->connector->send->numSharedComponents;++i) {
226     memcpy(&(coupler->send_buffer[(block_size)*i]),&(in[ block_size * coupler->connector->send->shared[i]]), block_size_size);
227     }
228     } else {
229     #pragma omp parallel for private(i)
230     for (i=0; i < coupler->connector->send->numSharedComponents;++i) coupler->send_buffer[i]=in[coupler->connector->send->shared[i]];
231 ksteube 1313 }
232     /* send buffer out */
233     {
234 gross 1552 for (i=0; i< coupler->connector->send->numNeighbors; ++i) {
235 jfenwick 3259 #ifdef ESYS_MPI
236 gross 1553 MPI_Issend(&(coupler->send_buffer[coupler->connector->send->offsetInShared[i] * block_size]),
237 gross 1552 (coupler->connector->send->offsetInShared[i+1]- coupler->connector->send->offsetInShared[i])*block_size,
238 ksteube 1313 MPI_DOUBLE,
239 gross 1552 coupler->connector->send->neighbor[i],
240 ksteube 1313 mpi_info->msg_tag_counter+mpi_info->rank,
241     mpi_info->comm,
242 gross 1553 &(coupler->mpi_requests[i+ coupler->connector->recv->numNeighbors]));
243 ksteube 1313 #endif
244     }
245     }
246 gross 1413 mpi_info->msg_tag_counter+=mpi_info->size;
247 ksteube 1313 }
248     }
249    
250     double* Paso_Coupler_finishCollect(Paso_Coupler* coupler)
251     {
252 jfenwick 3259 Esys_MPIInfo *mpi_info = coupler->mpi_info;
253 ksteube 1313 if ( mpi_info->size>1) {
254     /* wait for receive */
255 jfenwick 3259 #ifdef ESYS_MPI
256 gross 1552 MPI_Waitall(coupler->connector->recv->numNeighbors+coupler->connector->send->numNeighbors,
257 ksteube 1313 coupler->mpi_requests,
258     coupler->mpi_stati);
259     #endif
260     }
261 gross 1758
262 ksteube 1313 return coupler->recv_buffer;
263     }
264 gross 3442
265 gross 1562 void Paso_Coupler_copyAll(const Paso_Coupler* src, Paso_Coupler* target)
266     {
267     dim_t i;
268     #pragma omp parallel
269     {
270     #pragma omp for private(i)
271     for (i =0; i< src->connector->recv->numSharedComponents * src->block_size; ++i) {
272     target->recv_buffer[i]=src->recv_buffer[i];
273     }
274     #pragma omp for private(i)
275     for (i =0; i< Paso_Coupler_getLocalLength(src) * src->block_size; ++i) {
276     target->data[i]=src->data[i];
277     }
278     }
279     }
280 gross 3442
281 gross 3445 /* aadds a local vector x * a to the vector y including overlap :*/
282     void Paso_Coupler_add(const dim_t n, double* x, const double a, const double* y, Paso_Coupler *coupler)
283 gross 3442 {
284     double *remote_values = NULL;
285     const dim_t overlap_n = Paso_Coupler_getNumOverlapValues(coupler) ;
286 gross 3445 const dim_t my_n= n - overlap_n;
287 gross 3442 dim_t i;
288    
289     if (ABS(a) > 0 ) {
290     Paso_Coupler_startCollect(coupler, x);
291    
292     if ( a == 1.) {
293     #pragma omp parallel for private(i)
294     for (i=0; i<my_n; i++) {
295     x[i]+=y[i];
296     }
297     } else {
298     #pragma omp parallel for private(i)
299     for (i=0; i<my_n; i++) {
300     x[i]+=a*y[i];
301     }
302     }
303    
304     Paso_Coupler_finishCollect(coupler);
305     remote_values=coupler->recv_buffer;
306    
307     if ( a == 1.) {
308     #pragma omp parallel for private(i)
309     for (i=0;i<overlap_n; ++i) {
310 gross 3445 x[my_n+i]+=remote_values[i];
311 gross 3442 }
312     } else {
313     #pragma omp parallel for private(i)
314     for (i=0;i<overlap_n; ++i) {
315 gross 3445 x[my_n+i]+=a*remote_values[i];
316 gross 3442 }
317     }
318     }
319 gross 3443 }
320 gross 3445
321     /* adjusts max values across shared values x */
322     void Paso_Coupler_max(const dim_t n, double* x, Paso_Coupler *coupler)
323     {
324     double *remote_values = NULL;
325     const dim_t overlap_n = Paso_Coupler_getNumOverlapValues(coupler) ;
326     const dim_t my_n= n - overlap_n;
327     dim_t i;
328    
329     Paso_Coupler_startCollect(coupler, x);
330     Paso_Coupler_finishCollect(coupler);
331     remote_values=coupler->recv_buffer;
332     #pragma omp parallel for private(i)
333     for (i=0;i<overlap_n; ++i) x[my_n+i]=MAX(x[my_n+i], remote_values[i]);
334     }

  ViewVC Help
Powered by ViewVC 1.1.26