/[escript]/trunk/paso/src/SystemMatrix_copyRemoteCoupleBlock.cpp
ViewVC logotype

Contents of /trunk/paso/src/SystemMatrix_copyRemoteCoupleBlock.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4829 - (show annotations)
Thu Apr 3 04:02:53 2014 UTC (5 years, 5 months ago) by caltinay
File size: 10458 byte(s)
checkpointing some SparseMatrix cleanup.

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2014 by University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17
18 /****************************************************************************/
19 /* Paso: SystemMatrix */
20 /* */
21 /* Copy mainBlock and col_coupleBlock in other ranks */
22 /* into remote_coupleBlock */
23 /* */
24 /* WARNING: function uses mpi_requests of the coupler attached to A. */
25 /* */
26 /****************************************************************************/
27
28 /* Copyrights by ACcESS Australia 2003 */
29 /* Author: Lin Gao, l.gao@uq.edu.au */
30
31 /****************************************************************************/
32
33 #include "SystemMatrix.h"
34
35 void Paso_SystemMatrix_copyRemoteCoupleBlock(Paso_SystemMatrix* A, bool recreatePattern)
36 {
37 paso::Coupler_ptr coupler;
38 paso::SharedComponents_ptr send, recv;
39 double *cols=NULL, *send_buf=NULL;
40 index_t *global_id=NULL, *cols_array=NULL, *ptr_ptr=NULL, *ptr_idx=NULL;
41 index_t *send_idx=NULL, *send_offset=NULL, *recv_buf=NULL, *recv_offset=NULL;
42 index_t i, j, k, l, m, n, p, q, j_ub, k_lb, k_ub, l_lb, l_ub, i0;
43 index_t offset, len, overlapped_n, base, block_size, block_size_size;
44 index_t num_main_cols, num_couple_cols, num_neighbors, row, neighbor;
45 index_t *recv_degree=NULL, *send_degree=NULL;
46 dim_t rank=A->mpi_info->rank, mpi_size=A->mpi_info->size;
47
48 if (mpi_size == 1)
49 return;
50
51 if (recreatePattern)
52 A->remote_coupleBlock.reset();
53
54 if (A->remote_coupleBlock)
55 return;
56
57 /* sending/receiving unknown's global ID */
58 num_main_cols = A->mainBlock->numCols;
59 cols = new double[num_main_cols];
60 offset = A->col_distribution->first_component[rank];
61 #pragma omp parallel for private(i) schedule(static)
62 for (i=0; i<num_main_cols; ++i) cols[i] = offset + i;
63 if (A->global_id == NULL) {
64 coupler.reset(new paso::Coupler(A->col_coupler->connector, 1));
65 coupler->startCollect(cols);
66 }
67
68 recv_buf = new index_t[mpi_size];
69 recv_degree = new index_t[mpi_size];
70 recv_offset = new index_t[mpi_size+1];
71 #pragma omp parallel for private(i) schedule(static)
72 for (i=0; i<mpi_size; i++){
73 recv_buf[i] = 0;
74 recv_degree[i] = 1;
75 recv_offset[i] = i;
76 }
77
78 num_couple_cols = A->col_coupleBlock->numCols;
79 overlapped_n = A->row_coupleBlock->numRows;
80 send = A->row_coupler->connector->send;
81 recv = A->row_coupler->connector->recv;
82 num_neighbors = send->numNeighbors;
83 block_size = A->block_size;
84 block_size_size = block_size * sizeof(double);
85
86 /* waiting for receiving unknown's global ID */
87 if (A->global_id == NULL) {
88 coupler->finishCollect();
89 global_id = new index_t[num_couple_cols+1];
90 #pragma omp parallel for private(i) schedule(static)
91 for (i=0; i<num_couple_cols; ++i)
92 global_id[i] = coupler->recv_buffer[i];
93 A->global_id = global_id;
94 coupler.reset();
95 } else
96 global_id = A->global_id;
97
98 /* distribute the number of cols in current col_coupleBlock to all ranks */
99 #ifdef ESYS_MPI
100 MPI_Allgatherv(&num_couple_cols, 1, MPI_INT, recv_buf, recv_degree, recv_offset, MPI_INT, A->mpi_info->comm);
101 #endif
102
103 /* distribute global_ids of cols to be considered to all ranks*/
104 len = 0;
105 for (i=0; i<mpi_size; i++){
106 recv_degree[i] = recv_buf[i];
107 recv_offset[i] = len;
108 len += recv_buf[i];
109 }
110 recv_offset[mpi_size] = len;
111 cols_array = new index_t[len];
112
113 #ifdef ESYS_MPI
114 MPI_Allgatherv(global_id, num_couple_cols, MPI_INT, cols_array, recv_degree, recv_offset, MPI_INT, A->mpi_info->comm);
115 #endif
116
117 /* first, prepare the ptr_ptr to be received */
118 ptr_ptr = new index_t[overlapped_n+1];
119 for (p=0; p<recv->numNeighbors; p++) {
120 row = recv->offsetInShared[p];
121 i = recv->offsetInShared[p+1];
122 #ifdef ESYS_MPI
123 MPI_Irecv(&(ptr_ptr[row]), i-row, MPI_INT, recv->neighbor[p],
124 A->mpi_info->msg_tag_counter+recv->neighbor[p],
125 A->mpi_info->comm,
126 &(A->row_coupler->mpi_requests[p]));
127 #endif
128 }
129
130 /* now prepare the rows to be sent (the degree, the offset and the data) */
131 p = send->offsetInShared[num_neighbors];
132 len = 0;
133 for (i=0; i<num_neighbors; i++) {
134 /* #cols per row X #rows */
135 len += recv_buf[send->neighbor[i]] *
136 (send->offsetInShared[i+1] - send->offsetInShared[i]);
137 }
138 send_buf = new double[len*block_size];
139 send_idx = new index_t[len];
140 send_offset = new index_t[p+1];
141 send_degree = new index_t[num_neighbors];
142
143 len = 0;
144 base = 0;
145 i0 = 0;
146 for (p=0; p<num_neighbors; p++) {
147 i = i0;
148 neighbor = send->neighbor[p];
149 l_ub = recv_offset[neighbor+1];
150 l_lb = recv_offset[neighbor];
151 j_ub = send->offsetInShared[p + 1];
152 for (j=send->offsetInShared[p]; j<j_ub; j++) {
153 row = send->shared[j];
154
155 /* check col_coupleBlock for data to be passed @row */
156 l = l_lb;
157 k_ub = A->col_coupleBlock->pattern->ptr[row+1];
158 k = A->col_coupleBlock->pattern->ptr[row];
159 q = A->mainBlock->pattern->index[A->mainBlock->pattern->ptr[row]] + offset;
160 while (k<k_ub && l<l_ub) {
161 m = global_id[A->col_coupleBlock->pattern->index[k]];
162 if (m > q) break;
163 n = cols_array[l];
164 if (m == n) {
165 send_idx[len] = l - l_lb;
166 memcpy(&(send_buf[len*block_size]), &(A->col_coupleBlock->val[block_size*k]), block_size_size);
167 len++;
168 l++;
169 k++;
170 } else if (m < n) {
171 k++;
172 } else {
173 l++;
174 }
175 }
176 k_lb = k;
177
178 /* check mainBlock for data to be passed @row */
179 k_ub = A->mainBlock->pattern->ptr[row+1];
180 k=A->mainBlock->pattern->ptr[row];
181 while (k<k_ub && l<l_ub) {
182 m = A->mainBlock->pattern->index[k] + offset;
183 n = cols_array[l];
184 if (m == n) {
185 send_idx[len] = l - l_lb;
186 memcpy(&(send_buf[len*block_size]), &(A->mainBlock->val[block_size*k]), block_size_size);
187 len++;
188 l++;
189 k++;
190 } else if (m < n) {
191 k++;
192 } else {
193 l++;
194 }
195 }
196
197 /* check col_coupleBlock for data to be passed @row */
198 k_ub = A->col_coupleBlock->pattern->ptr[row+1];
199 k=k_lb;
200 while (k<k_ub && l<l_ub) {
201 m = global_id[A->col_coupleBlock->pattern->index[k]];
202 n = cols_array[l];
203 if (m == n) {
204 send_idx[len] = l - l_lb;
205 memcpy(&(send_buf[len*block_size]), &(A->col_coupleBlock->val[block_size*k]), block_size_size);
206 len++;
207 l++;
208 k++;
209 } else if (m < n) {
210 k++;
211 } else {
212 l++;
213 }
214 }
215
216 send_offset[i] = len - base;
217 base = len;
218 i++;
219 }
220
221 /* sending */
222 #ifdef ESYS_MPI
223 MPI_Issend(&(send_offset[i0]), i-i0, MPI_INT, send->neighbor[p],
224 A->mpi_info->msg_tag_counter+rank,
225 A->mpi_info->comm,
226 &(A->row_coupler->mpi_requests[p+recv->numNeighbors]));
227 #endif
228 send_degree[p] = len;
229 i0 = i;
230 }
231
232 #ifdef ESYS_MPI
233 MPI_Waitall(A->row_coupler->connector->send->numNeighbors+A->row_coupler->connector->recv->numNeighbors,
234 A->row_coupler->mpi_requests,
235 A->row_coupler->mpi_stati);
236 #endif
237 ESYS_MPI_INC_COUNTER(*(A->mpi_info), mpi_size)
238
239 len = 0;
240 for (i=0; i<overlapped_n; i++) {
241 p = ptr_ptr[i];
242 ptr_ptr[i] = len;
243 len += p;
244 }
245 ptr_ptr[overlapped_n] = len;
246 ptr_idx = new index_t[len];
247
248 /* send/receive index array */
249 j=0;
250 for (p=0; p<recv->numNeighbors; p++) {
251 i = ptr_ptr[recv->offsetInShared[p+1]] - ptr_ptr[recv->offsetInShared[p]];
252 #ifdef ESYS_MPI
253 if (i > 0)
254 MPI_Irecv(&(ptr_idx[j]), i, MPI_INT, recv->neighbor[p],
255 A->mpi_info->msg_tag_counter+recv->neighbor[p],
256 A->mpi_info->comm,
257 &(A->row_coupler->mpi_requests[p]));
258 #endif
259 j += i;
260 }
261
262 j=0;
263 for (p=0; p<num_neighbors; p++) {
264 i = send_degree[p] - j;
265 #ifdef ESYS_MPI
266 if (i > 0)
267 MPI_Issend(&(send_idx[j]), i, MPI_INT, send->neighbor[p],
268 A->mpi_info->msg_tag_counter+rank,
269 A->mpi_info->comm,
270 &(A->row_coupler->mpi_requests[p+recv->numNeighbors]));
271 #endif
272 j = send_degree[p];
273 }
274
275 #ifdef ESYS_MPI
276 MPI_Waitall(A->row_coupler->connector->send->numNeighbors+A->row_coupler->connector->recv->numNeighbors,
277 A->row_coupler->mpi_requests,
278 A->row_coupler->mpi_stati);
279 #endif
280 ESYS_MPI_INC_COUNTER(*(A->mpi_info), mpi_size)
281
282 /* allocate pattern and sparsematrix for remote_coupleBlock */
283 paso::Pattern_ptr pattern(new paso::Pattern(A->row_coupleBlock->pattern->type,
284 overlapped_n, num_couple_cols, ptr_ptr, ptr_idx));
285 A->remote_coupleBlock.reset(new paso::SparseMatrix(A->row_coupleBlock->type,
286 pattern, A->row_block_size, A->col_block_size, false));
287
288 /* send/receive value array */
289 j=0;
290 for (p=0; p<recv->numNeighbors; p++) {
291 i = ptr_ptr[recv->offsetInShared[p+1]] - ptr_ptr[recv->offsetInShared[p]];
292 #ifdef ESYS_MPI
293 if (i > 0)
294 MPI_Irecv(&(A->remote_coupleBlock->val[j]), i * block_size,
295 MPI_DOUBLE, recv->neighbor[p],
296 A->mpi_info->msg_tag_counter+recv->neighbor[p],
297 A->mpi_info->comm,
298 &(A->row_coupler->mpi_requests[p]));
299 #endif
300 j += (i * block_size);
301 }
302
303 j=0;
304 for (p=0; p<num_neighbors; p++) {
305 i = send_degree[p] - j;
306 #ifdef ESYS_MPI
307 if (i > 0)
308 MPI_Issend(&(send_buf[j*block_size]), i*block_size, MPI_DOUBLE, send->neighbor[p],
309 A->mpi_info->msg_tag_counter+rank,
310 A->mpi_info->comm,
311 &(A->row_coupler->mpi_requests[p+recv->numNeighbors]));
312 #endif
313 j = send_degree[p];
314 }
315
316 #ifdef ESYS_MPI
317 MPI_Waitall(A->row_coupler->connector->send->numNeighbors+A->row_coupler->connector->recv->numNeighbors,
318 A->row_coupler->mpi_requests,
319 A->row_coupler->mpi_stati);
320 #endif
321 ESYS_MPI_INC_COUNTER(*(A->mpi_info), mpi_size)
322
323 /* release all temp memory allocation */
324 delete[] cols;
325 delete[] cols_array;
326 delete[] recv_offset;
327 delete[] recv_degree;
328 delete[] recv_buf;
329 delete[] send_buf;
330 delete[] send_offset;
331 delete[] send_degree;
332 delete[] send_idx;
333 }
334

Properties

Name Value
svn:mergeinfo /branches/amg_from_3530/paso/src/SystemMatrix_copyRemoteCoupleBlock.cpp:3531-3826 /branches/lapack2681/paso/src/SystemMatrix_copyRemoteCoupleBlock.cpp:2682-2741 /branches/pasowrap/paso/src/SystemMatrix_copyRemoteCoupleBlock.cpp:3661-3674 /branches/py3_attempt2/paso/src/SystemMatrix_copyRemoteCoupleBlock.cpp:3871-3891 /branches/restext/paso/src/SystemMatrix_copyRemoteCoupleBlock.cpp:2610-2624 /branches/ripleygmg_from_3668/paso/src/SystemMatrix_copyRemoteCoupleBlock.cpp:3669-3791 /branches/stage3.0/paso/src/SystemMatrix_copyRemoteCoupleBlock.cpp:2569-2590 /branches/symbolic_from_3470/paso/src/SystemMatrix_copyRemoteCoupleBlock.cpp:3471-3974 /branches/symbolic_from_3470/ripley/test/python/paso/src/SystemMatrix_copyRemoteCoupleBlock.cpp:3517-3974 /release/3.0/paso/src/SystemMatrix_copyRemoteCoupleBlock.cpp:2591-2601 /trunk/paso/src/SystemMatrix_copyRemoteCoupleBlock.cpp:4257-4344 /trunk/ripley/test/python/paso/src/SystemMatrix_copyRemoteCoupleBlock.cpp:3480-3515

  ViewVC Help
Powered by ViewVC 1.1.26