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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3443 - (hide annotations)
Tue Jan 18 05:38:48 2011 UTC (8 years, 8 months ago) by gross
File MIME type: text/plain
File size: 11078 byte(s)
test for AMG fixed. New test set guarantees SPD stiffness matrix. Previously this was not guaranteed and could lead to singular coarse level matrices (which cause MKL based runs to fail)
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     dim_t Paso_Coupler_getNumSharedValues(const Paso_Coupler* in) {
266 gross 3443 return in->connector->send->numSharedComponents * in->block_size;
267 gross 3442 }
268     dim_t Paso_Coupler_getNumOverlapValues(const Paso_Coupler* in) {
269 gross 3443 return in->connector->recv->numSharedComponents * in->block_size;
270 gross 3442 }
271    
272 gross 1562 dim_t Paso_Coupler_getLocalLength(const Paso_Coupler* in) {
273     return in->connector->send->local_length;
274     }
275     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     /* adds a local vector x * a to the vector y including overlap :*/
292     void Paso_Coupler_add(const dim_t my_n, double* x, const double a, const double* y, Paso_Coupler *coupler)
293     {
294     double *remote_values = NULL;
295     const dim_t overlap_n = Paso_Coupler_getNumOverlapValues(coupler) ;
296     const dim_t n= my_n + overlap_n;
297     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     x[i+overlap_n]+=remote_values[i];
321     }
322     } else {
323     #pragma omp parallel for private(i)
324     for (i=0;i<overlap_n; ++i) {
325     x[i+overlap_n]+=a*remote_values[i];
326     }
327     }
328     }
329 gross 3443 }

  ViewVC Help
Powered by ViewVC 1.1.26