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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3443 - (show 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
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
151 #ifdef ESYS_MPI
152 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 Esys_checkPtr(out->mpi_requests);
155 Esys_checkPtr(out->mpi_stati);
156 #endif
157 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 Esys_checkPtr(out->send_buffer);
161 Esys_checkPtr(out->recv_buffer);
162 }
163 }
164 if (Esys_noError()) {
165 return out;
166 } else {
167 Paso_Coupler_free(out);
168 return NULL;
169 }
170 }
171
172 /* returns a reference to Coupler */
173
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 Paso_Connector_free(in->connector);
188 MEMFREE(in->send_buffer);
189 MEMFREE(in->recv_buffer);
190 MEMFREE(in->mpi_requests);
191 MEMFREE(in->mpi_stati);
192 Esys_MPIInfo_free(in->mpi_info);
193 MEMFREE(in);
194 }
195 }
196 }
197
198
199 void Paso_Coupler_startCollect(Paso_Coupler* coupler,const double* in)
200 {
201 Esys_MPIInfo *mpi_info = coupler->mpi_info;
202 dim_t block_size=coupler->block_size;
203 size_t block_size_size=block_size*sizeof(double);
204 dim_t i;
205 coupler->data=(double*) in;
206 if ( mpi_info->size>1) {
207 /* start reveiving input */
208 {
209 for (i=0; i< coupler->connector->recv->numNeighbors; ++i) {
210 #ifdef ESYS_MPI
211 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 MPI_DOUBLE,
214 coupler->connector->recv->neighbor[i],
215 mpi_info->msg_tag_counter+coupler->connector->recv->neighbor[i],
216 mpi_info->comm,
217 &(coupler->mpi_requests[i]));
218 #endif
219
220 }
221 }
222 /* collect values into buffer */
223 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 }
232 /* send buffer out */
233 {
234 for (i=0; i< coupler->connector->send->numNeighbors; ++i) {
235 #ifdef ESYS_MPI
236 MPI_Issend(&(coupler->send_buffer[coupler->connector->send->offsetInShared[i] * block_size]),
237 (coupler->connector->send->offsetInShared[i+1]- coupler->connector->send->offsetInShared[i])*block_size,
238 MPI_DOUBLE,
239 coupler->connector->send->neighbor[i],
240 mpi_info->msg_tag_counter+mpi_info->rank,
241 mpi_info->comm,
242 &(coupler->mpi_requests[i+ coupler->connector->recv->numNeighbors]));
243 #endif
244 }
245 }
246 mpi_info->msg_tag_counter+=mpi_info->size;
247 }
248 }
249
250 double* Paso_Coupler_finishCollect(Paso_Coupler* coupler)
251 {
252 Esys_MPIInfo *mpi_info = coupler->mpi_info;
253 if ( mpi_info->size>1) {
254 /* wait for receive */
255 #ifdef ESYS_MPI
256 MPI_Waitall(coupler->connector->recv->numNeighbors+coupler->connector->send->numNeighbors,
257 coupler->mpi_requests,
258 coupler->mpi_stati);
259 #endif
260 }
261
262 return coupler->recv_buffer;
263 }
264
265 dim_t Paso_Coupler_getNumSharedValues(const Paso_Coupler* in) {
266 return in->connector->send->numSharedComponents * in->block_size;
267 }
268 dim_t Paso_Coupler_getNumOverlapValues(const Paso_Coupler* in) {
269 return in->connector->recv->numSharedComponents * in->block_size;
270 }
271
272 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
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 }

  ViewVC Help
Powered by ViewVC 1.1.26