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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3482 - (show annotations)
Wed Mar 23 04:06:52 2011 UTC (8 years, 5 months ago) by gross
File MIME type: text/plain
File size: 14038 byte(s)
some work on AMG MPI
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 /**************************************************************
16
17 Paso: SystemMatrix: copies the col_coupleBlock into
18 row_coupleBlock.
19
20 WARNING: function uses mpi_reqests of the coupler attached to A.
21
22 No reordering on the colums received is performed.
23 In practice this means that components in
24 A->row_coupleBlock->pattern->index
25 and
26 A->row_coupler->connector->recv->shared
27 are ordered by increasing value.
28
29 Notice: that send and receive A->row_coupler->connectors
30 are swapping roles.
31
32 **************************************************************/
33
34 /* Author: Lutz Gross, l.gross@uq.edu.au */
35
36 /**************************************************************/
37
38 #include "SystemMatrix.h"
39 #include "esysUtils/error.h"
40
41 /**************************************************************/
42
43 void Paso_SystemMatrix_copyColCoupleBlock(Paso_SystemMatrix *A)
44 {
45 /* const dim_t n=Paso_SystemMatrix_getNumRows(A); */
46 dim_t p;
47 index_t z0, iPtr, rPtr;
48 const dim_t block_size=A->block_size;
49 const size_t block_size_size=block_size*sizeof(double);
50 double * send_buffer = NULL;
51
52 if ( A->mpi_info->size > 1 ) {
53 if ( A->row_coupleBlock == NULL ) {
54 Esys_setError(VALUE_ERROR, "SystemMatrix_copyColCoupleBlock: creation of row_coupleBlock pattern not supported yet.");
55 return;
56 }
57 if ( A->row_coupler->in_use ) {
58 Esys_setError(SYSTEM_ERROR,"Paso_Coupler_finishCollect: Communication has not been initiated.");
59 return;
60 }
61
62 /* start receiving */
63 for (p=0; p<A->row_coupler->connector->recv->numNeighbors; p++) {
64 #ifdef ESYS_MPI
65 const index_t irow1= A->row_coupler->connector->recv->offsetInShared[p];
66 const index_t irow2= A->row_coupler->connector->recv->offsetInShared[p+1];
67 const index_t a = A->row_coupleBlock->pattern->ptr[irow1];
68 const index_t b = A->row_coupleBlock->pattern->ptr[irow2];
69 printf(" %d %d : %d %d : %d %d\n",A->row_coupler->connector->recv->offsetInShared[p], A->row_coupler->connector->recv->offsetInShared[p+1],a,b,irow1,irow2);
70 printf("reveive from %d len = %d\n",A->row_coupler->connector->recv->neighbor[p], (b-a) * block_size);
71
72 MPI_Irecv(&(A->row_coupleBlock->val[a * block_size]), (b-a) * block_size, MPI_DOUBLE,
73 A->row_coupler->connector->recv->neighbor[p],
74 A->mpi_info->msg_tag_counter+A->row_coupler->connector->recv->neighbor[p],
75 A->mpi_info->comm,
76 &(A->row_coupler->mpi_requests[p]) );
77
78 #endif
79 }
80 /* start sending */
81 z0=0;
82 send_buffer = TMPMEMALLOC(A->col_coupleBlock->len, double);
83
84 for (p=0; p<A->row_coupler->connector->send->numNeighbors; p++) {
85 /* j_min, j_max defines the range of columns to be sent to processor p*/
86 const index_t j_min = A->col_coupler->connector->recv->offsetInShared[p];
87 const index_t j_max = A->col_coupler->connector->recv->offsetInShared[p+1];
88 index_t z = z0;
89
90 /* run over the rows to be connected to processor p */
91 for (rPtr= A->row_coupler->connector->send->offsetInShared[p]; rPtr< A->row_coupler->connector->send->offsetInShared[p+1]; ++rPtr) {
92 const index_t row=A->row_coupler->connector->send->shared[rPtr];
93
94 /* collect the entries in the col. couple block refering to columns on processor p */
95 for (iPtr =A->col_coupleBlock->pattern->ptr[row]; iPtr<A->col_coupleBlock->pattern->ptr[row+1]; ++iPtr) {
96 register index_t j =A->col_coupleBlock->pattern->index[iPtr];
97 if ( (j_min <= j) && (j < j_max) ) {
98 memcpy(&(send_buffer[z]),&(A->col_coupleBlock->val[ block_size * iPtr]), block_size_size);
99 z+=block_size;
100 }
101 }
102
103 }
104 #ifdef ESYS_MPI
105 printf("send to %d len = %d\n",A->row_coupler->connector->send->neighbor[p], z-z0);
106 MPI_Issend(&(send_buffer[z0]),z-z0, MPI_DOUBLE,
107 A->row_coupler->connector->send->neighbor[p],
108 A->mpi_info->msg_tag_counter+A->mpi_info->rank,
109 A->mpi_info->comm,
110 &(A->row_coupler->mpi_requests[p+A->row_coupler->connector->recv->numNeighbors]));
111
112 #endif
113 }
114
115
116 /* wait til everything is done */
117 #ifdef ESYS_MPI
118 MPI_Waitall(A->row_coupler->connector->send->numNeighbors+A->row_coupler->connector->recv->numNeighbors,
119 A->row_coupler->mpi_requests,
120 A->row_coupler->mpi_stati);
121 #endif
122 A->mpi_info->msg_tag_counter+=A->mpi_info->size;
123 TMPMEMFREE(send_buffer);
124 }
125 return;
126 }
127
128 #ifdef AAAA
129 void Paso_SystemMatrix_copyRemoteCoupleBlock(Paso_SystemMatrix *A, const bool_t recreatePattern)
130 {
131 Paso_Pattern* couple_pattern=Null;
132 /* const dim_t n=Paso_SystemMatrix_getNumRows(A); */
133 dim_t p;
134 index_t z0, iPtr, rPtr;
135 const dim_t block_size=A->block_size;
136 const size_t block_size_size=block_size*sizeof(double);
137 double * send_buffer = NULL;
138
139 /* for each processor p we need to identify the rows which have connection to processor p.
140 this list is available at
141
142 A ->row_coupler->connector->send->shared[rPtr...] rptr=A->row_coupler->connector->send->offsetInShared[p]*
143
144 These rows correspond to certain block of rows in the A->row_coupleBlock and A->remote_coupleBlock on processor p
145 where (on p) the first row is given by A->row_coupler->connector->recv->offsetInShared[q] wher q is my processor.
146
147 for each row conencted to p we need to ideinfy the list of columns in mainBlock which are used by processor p
148 this list is given by
149
150 A ->col_coupler->connector->send->shared[rPtr...] rptr=A->col_coupler->connector->send->offsetInShared[p]
151
152 where the order in the list corresponds to the order in which the columns are ordered on processor p.
153
154
155 the algorithm:
156
157 for each processor p
158
159 a) mark the the columns used by processor p in mask by new label on processor starting from 0
160 b) collect the columnsin mainBlock to be send to processor p:
161 - for each row the list of columns (defines pattern)
162 - the corresponding entries in mainBlock
163 c) send length of pattern to processor p
164 d) send pattern to processor p
165 e) send entries to processor p
166
167 receive all pattern length
168 allocate pattern
169 receive pattern
170 merge patterns ?????
171 allocate pattern and matrix ?????
172
173
174 */
175
176
177 for these rows we need to find the
178 ==========================
179 if ( A->mpi_info->size > 1 ) {
180
181 if (recreatePattern) Paso_SparseMatrix_free(in->remote_coupleBlock);
182
183
184 if (A->remote_coupleBlock == NULL) {
185
186
187
188
189 A->remote_coupleBlock=Paso_SparseMatrix_alloc(Paso_SparseMatrixType,
190 couple_pattern,
191 A->row_block_size, A->col_block_size,
192 const bool_t patternIsUnrolled)
193
194 Paso_Pattern_free(couple_pattern);
195 }
196 /* and now the matrix : */
197 if (! Esys_noError()) {
198 couple_pattern=Paso_Pattern_getReference(remote_coupleBlock->pattern);
199 /* start receiving */
200 for (p=0; p<A->row_coupler->connector->recv->numNeighbors; p++) {
201 #ifdef ESYS_MPI
202 const index_t irow1= A->row_coupler->connector->recv->offsetInShared[p];
203 const index_t irow2= A->row_coupler->connector->recv->offsetInShared[p+1];
204 const index_t a = A->row_coupleBlock->pattern->ptr[irow1];
205 const index_t b = A->row_coupleBlock->pattern->ptr[irow2];
206 printf(" %d %d : %d %d : %d %d\n",A->row_coupler->connector->recv->offsetInShared[p], A->row_coupler->connector->recv->offsetInShared[p+1],a,b,irow1,irow2);
207 printf("reveive from %d len = %d\n",A->row_coupler->connector->recv->neighbor[p], (b-a) * block_size);
208
209 MPI_Irecv(&(A->row_coupleBlock->val[a * block_size]), (b-a) * block_size, MPI_DOUBLE,
210 A->row_coupler->connector->recv->neighbor[p],
211 A->mpi_info->msg_tag_counter+A->row_coupler->connector->recv->neighbor[p],
212 A->mpi_info->comm,
213 &(A->row_coupler->mpi_requests[p]) );
214
215 #endif
216 }
217 /* start sending */
218 z0=0;
219 send_buffer = TMPMEMALLOC(A->col_coupleBlock->len, double);
220
221 for (p=0; p<A->row_coupler->connector->send->numNeighbors; p++) {
222 /* j_min, j_max defines the range of columns to be sent to processor p*/
223 const index_t j_min = A->col_coupler->connector->recv->offsetInShared[p];
224 const index_t j_max = A->col_coupler->connector->recv->offsetInShared[p+1];
225 index_t z = z0;
226
227 /* run over the rows to be connected to processor p */
228 for (rPtr= A->row_coupler->connector->send->offsetInShared[p]; rPtr< A->row_coupler->connector->send->offsetInShared[p+1]; ++rPtr) {
229 const index_t row=A->row_coupler->connector->send->shared[rPtr];
230
231 /* collect the entries in the col. couple block refering to columns on processor p */
232 for (iPtr =A->col_coupleBlock->pattern->ptr[row]; iPtr<A->col_coupleBlock->pattern->ptr[row+1]; ++iPtr) {
233 register index_t j =A->col_coupleBlock->pattern->index[iPtr];
234 if ( (j_min <= j) && (j < j_max) ) {
235 memcpy(&(send_buffer[z]),&(A->col_coupleBlock->val[ block_size * iPtr]), block_size_size);
236 z+=block_size;
237 }
238 }
239 }
240 #ifdef ESYS_MPI
241 printf("send to %d len = %d\n",A->row_coupler->connector->send->neighbor[p], z-z0);
242 MPI_Issend(&(send_buffer[z0]),z-z0, MPI_DOUBLE,
243 A->row_coupler->connector->send->neighbor[p],
244 A->mpi_info->msg_tag_counter+A->mpi_info->rank,
245 A->mpi_info->comm,
246 &(A->row_coupler->mpi_requests[p+A->row_coupler->connector->recv->numNeighbors]));
247
248 #endif
249 }
250
251
252 /* wait til everything is done */
253 #ifdef ESYS_MPI
254 MPI_Waitall(A->row_coupler->connector->send->numNeighbors+A->row_coupler->connector->recv->numNeighbors,
255 A->row_coupler->mpi_requests,
256 A->row_coupler->mpi_stati);
257 #endif
258 A->mpi_info->msg_tag_counter+=A->mpi_info->size;
259 TMPMEMFREE(send_buffer);
260
261
262 Paso_Pattern_free(couple_pattern);
263 }
264
265 }
266
267 if ( A->row_coupleBlock == NULL ) {
268 Esys_setError(VALUE_ERROR, "SystemMatrix_copyColCoupleBlock: creation of row_coupleBlock pattern not supported yet.");
269 return;
270 }
271
272
273 if ( A->row_coupleBlock == NULL ) {
274 Esys_setError(VALUE_ERROR, "SystemMatrix_copyColCoupleBlock: creation of row_coupleBlock pattern not supported yet.");
275 return;
276 }
277
278 if ( A->row_coupler->in_use ) {
279 Esys_setError(SYSTEM_ERROR,"Paso_Coupler_finishCollect: Communication has not been initiated.");
280 return;
281 }
282
283
284 /* start receiving */
285 for (p=0; p<A->row_coupler->connector->recv->numNeighbors; p++) {
286 #ifdef ESYS_MPI
287 const index_t irow1= A->row_coupler->connector->recv->offsetInShared[p];
288 const index_t irow2= A->row_coupler->connector->recv->offsetInShared[p+1];
289 const index_t a = A->row_coupleBlock->pattern->ptr[irow1];
290 const index_t b = A->row_coupleBlock->pattern->ptr[irow2];
291 printf(" %d %d : %d %d : %d %d\n",A->row_coupler->connector->recv->offsetInShared[p], A->row_coupler->connector->recv->offsetInShared[p+1],a,b,irow1,irow2);
292 printf("reveive from %d len = %d\n",A->row_coupler->connector->recv->neighbor[p], (b-a) * block_size);
293
294 MPI_Irecv(&(A->row_coupleBlock->val[a * block_size]), (b-a) * block_size, MPI_DOUBLE,
295 A->row_coupler->connector->recv->neighbor[p],
296 A->mpi_info->msg_tag_counter+A->row_coupler->connector->recv->neighbor[p],
297 A->mpi_info->comm,
298 &(A->row_coupler->mpi_requests[p]) );
299
300 #endif
301 }
302 /* start sending */
303 z0=0;
304 send_buffer = TMPMEMALLOC(A->col_coupleBlock->len, double);
305
306 for (p=0; p<A->row_coupler->connector->send->numNeighbors; p++) {
307 const index_t j_min = A->col_coupler->connector->recv->offsetInShared[p];
308 const index_t j_max = A->col_coupler->connector->recv->offsetInShared[p+1];
309 index_t z = z0;
310
311 for (rPtr= A->row_coupler->connector->send->offsetInShared[p]; rPtr< A->row_coupler->connector->send->offsetInShared[p+1]; ++rPtr) {
312
313 const index_t row=A->row_coupler->connector->send->shared[rPtr];
314
315
316 for (iPtr =A->col_coupleBlock->pattern->ptr[row]; iPtr<A->col_coupleBlock->pattern->ptr[row+1]; ++iPtr) {
317 register index_t j =A->col_coupleBlock->pattern->index[iPtr];
318 if ( (j_min <= j) && (j < j_max) ) {
319 memcpy(&(send_buffer[z]),&(A->col_coupleBlock->val[ block_size * iPtr]), block_size_size);
320 z+=block_size;
321 }
322 }
323 }
324 #ifdef ESYS_MPI
325 printf("send to %d len = %d\n",A->row_coupler->connector->send->neighbor[p], z-z0);
326 MPI_Issend(&(send_buffer[z0]),z-z0, MPI_DOUBLE,
327 A->row_coupler->connector->send->neighbor[p],
328 A->mpi_info->msg_tag_counter+A->mpi_info->rank,
329 A->mpi_info->comm,
330 &(A->row_coupler->mpi_requests[p+A->row_coupler->connector->recv->numNeighbors]));
331
332 #endif
333 }
334
335
336 /* wait til everything is done */
337 #ifdef ESYS_MPI
338 MPI_Waitall(A->row_coupler->connector->send->numNeighbors+A->row_coupler->connector->recv->numNeighbors,
339 A->row_coupler->mpi_requests,
340 A->row_coupler->mpi_stati);
341 #endif
342 A->mpi_info->msg_tag_counter+=A->mpi_info->size;
343 TMPMEMFREE(send_buffer);
344 }
345 return;
346 }
347 #endif

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26