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

Diff of /trunk/paso/src/AMG_Interpolation.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3763 by lgao, Tue Jan 10 05:10:17 2012 UTC revision 3767 by lgao, Fri Jan 13 01:52:46 2012 UTC
# Line 39  Line 39 
39  ***************************************************************/  ***************************************************************/
40    
41  #define MY_DEBUG 0  #define MY_DEBUG 0
42  #define MY_DEBUG1 1  #define MY_DEBUG1 0
43    
44  /* Extend system matrix B with extra two sparse matrixes:  /* Extend system matrix B with extra two sparse matrixes:
45      B_ext_main and B_ext_couple      B_ext_main and B_ext_couple
# Line 116  if (MY_DEBUG1) fprintf(stderr, "CP1_1\n" Line 116  if (MY_DEBUG1) fprintf(stderr, "CP1_1\n"
116    recv = A->col_coupler->connector->recv;    recv = A->col_coupler->connector->recv;
117    num_neighbors = send->numNeighbors;    num_neighbors = send->numNeighbors;
118    p = send->offsetInShared[num_neighbors];    p = send->offsetInShared[num_neighbors];
119    len = p * num_main_cols * 2;   /*XXX used to be max_num_cols, do I mean num_main_cols?? */    len = p * B->col_distribution->first_component[size];
120    send_buf = TMPMEMALLOC(len * block_size, double);    send_buf = TMPMEMALLOC(len * block_size, double);
121    send_idx = TMPMEMALLOC(len, index_t);    send_idx = TMPMEMALLOC(len, index_t);
122    send_offset = TMPMEMALLOC((p+1)*2, index_t);    send_offset = TMPMEMALLOC((p+1)*2, index_t);
# Line 178  if (MY_DEBUG) { Line 178  if (MY_DEBUG) {
178    for (p=0; p<q; p++) {    for (p=0; p<q; p++) {
179      row = recv->offsetInShared[p];      row = recv->offsetInShared[p];
180      m = recv->offsetInShared[p + 1];      m = recv->offsetInShared[p + 1];
181      MPI_Irecv(&(ptr_ptr[row]), 2 * (m-row), MPI_INT, recv->neighbor[p],      MPI_Irecv(&(ptr_ptr[2*row]), 2 * (m-row), MPI_INT, recv->neighbor[p],
182                  A->mpi_info->msg_tag_counter+recv->neighbor[p],                  A->mpi_info->msg_tag_counter+recv->neighbor[p],
183                  A->mpi_info->comm,                  A->mpi_info->comm,
184                  &(A->col_coupler->mpi_requests[p]));                  &(A->col_coupler->mpi_requests[p]));
185    if(MY_DEBUG1) fprintf(stderr, "rank %d receive %d from %d tag %d (numNeighbors %d)\n", rank, 2*(m-row), recv->neighbor[p], A->mpi_info->msg_tag_counter+recv->neighbor[p], q);
186    }    }
187    
188    /* now prepare the rows to be sent (the degree, the offset and the data) */    /* now prepare the rows to be sent (the degree, the offset and the data) */
# Line 192  if (MY_DEBUG) { Line 193  if (MY_DEBUG) {
193      m_lb = B->col_distribution->first_component[neighbor];      m_lb = B->col_distribution->first_component[neighbor];
194      m_ub = B->col_distribution->first_component[neighbor + 1];      m_ub = B->col_distribution->first_component[neighbor + 1];
195      j_ub = send->offsetInShared[p + 1];      j_ub = send->offsetInShared[p + 1];
196    if (MY_DEBUG1) fprintf(stderr, "rank%d neighbor %d m_lb %d m_ub %d offset %d\n", rank, neighbor, m_lb, m_ub, offset);
197      for (j=send->offsetInShared[p]; j<j_ub; j++) {      for (j=send->offsetInShared[p]; j<j_ub; j++) {
198      row = send->shared[j];      row = send->shared[j];
199      l_m = 0;      l_m = 0;
# Line 205  if (MY_DEBUG) { Line 207  if (MY_DEBUG) {
207        m = global_id[B->col_coupleBlock->pattern->index[k]];        m = global_id[B->col_coupleBlock->pattern->index[k]];
208  //    if (m > q) break;  //    if (m > q) break;
209        if (m > offset) break;        if (m > offset) break;
210    if (MY_DEBUG) fprintf(stderr, "rank%d (1) row %d m %d k %d\n", rank, row, m, B->col_coupleBlock->pattern->index[k]);
211        if (m>= m_lb && m < m_ub) {        if (m>= m_lb && m < m_ub) {
212          /* data to be passed to sparse matrix B_ext_main */          /* data to be passed to sparse matrix B_ext_main */
213          idx_m[l_m] = m - m_lb;          idx_m[l_m] = m - m_lb;
# Line 227  if (MY_DEBUG) fprintf(stderr, "rank%d (1 Line 230  if (MY_DEBUG) fprintf(stderr, "rank%d (1
230      memcpy(&(send_c[l_c*block_size]), &(B->mainBlock->val[block_size*k]), block_size_size * (k_ub-k));      memcpy(&(send_c[l_c*block_size]), &(B->mainBlock->val[block_size*k]), block_size_size * (k_ub-k));
231      for (; k<k_ub; k++) {      for (; k<k_ub; k++) {
232        m = B->mainBlock->pattern->index[k] + offset;        m = B->mainBlock->pattern->index[k] + offset;
233    if (MY_DEBUG) fprintf(stderr, "rank%d (2) row %d m %d k %d\n", rank, row, m, B->mainBlock->pattern->index[k]);
234  if (MY_DEBUG) fprintf(stderr, "rank%d (2) m %d m_lb %d m_ub %d\n", rank, m, m_lb, m_ub);  if (MY_DEBUG) fprintf(stderr, "rank%d (2) m %d m_lb %d m_ub %d\n", rank, m, m_lb, m_ub);
235        idx_c[l_c] = m;        idx_c[l_c] = m;
236        l_c++;        l_c++;
# Line 238  if (MY_DEBUG) fprintf(stderr, "rank%d (2 Line 242  if (MY_DEBUG) fprintf(stderr, "rank%d (2
242      k_ub = B->col_coupleBlock->pattern->ptr[row + 1];      k_ub = B->col_coupleBlock->pattern->ptr[row + 1];
243      for (k=k_lb; k<k_ub; k++) {      for (k=k_lb; k<k_ub; k++) {
244            m = global_id[B->col_coupleBlock->pattern->index[k]];            m = global_id[B->col_coupleBlock->pattern->index[k]];
245        if (m>= m_lb && i < m_ub) {  if (MY_DEBUG) fprintf(stderr, "rank%d (3) row %d m %d k %d\n", rank, row, m, B->col_coupleBlock->pattern->index[k]);
246          if (m>= m_lb && m < m_ub) {
247          /* data to be passed to sparse matrix B_ext_main */          /* data to be passed to sparse matrix B_ext_main */
248          idx_m[l_m] = m - m_lb;          idx_m[l_m] = m - m_lb;
249          memcpy(&(send_m[l_m*block_size]), &(B->col_coupleBlock->val[block_size*k]), block_size_size);          memcpy(&(send_m[l_m*block_size]), &(B->col_coupleBlock->val[block_size*k]), block_size_size);
# Line 252  if (MY_DEBUG) fprintf(stderr, "rank%d (3 Line 257  if (MY_DEBUG) fprintf(stderr, "rank%d (3
257        }        }
258      }      }
259    
260  if (MY_DEBUG) fprintf(stderr, "rank %d p %d i %d j %d l_c %d l_m %d ub %d\n", rank, p, i, j, l_c, l_m, num_couple_cols + num_main_cols);  if (MY_DEBUG && rank == 2) fprintf(stderr, "rank %d p %d i %d j %d l_c %d l_m %d ub %d\n", rank, p, i, j, l_c, l_m, num_couple_cols + num_main_cols);
261      memcpy(&(send_buf[len]), &(send_m[0]), block_size_size*l_m);      memcpy(&(send_buf[len]), send_m, block_size_size*l_m);
262      memcpy(&(send_idx[len]), &(idx_m[0]), l_m * sizeof(index_t));      memcpy(&(send_idx[len]), idx_m, l_m * sizeof(index_t));
263          send_offset[2*i] = l_m;          send_offset[2*i] = l_m;
264      len += l_m;      len += l_m;
265      memcpy(&(send_buf[len]), &(send_c[0]), block_size_size*l_c);      memcpy(&(send_buf[len]), send_c, block_size_size*l_c);
266      memcpy(&(send_idx[len]), &(idx_c[0]), l_c * sizeof(index_t));      memcpy(&(send_idx[len]), idx_c, l_c * sizeof(index_t));
267      send_offset[2*i+1] = l_c;      send_offset[2*i+1] = l_c;
268      len += l_c;      len += l_c;
269        if (MY_DEBUG && rank == 0) {
270    int sum = l_m+l_c, my_i;
271    char *str1, *str2;
272    str1 = TMPMEMALLOC(sum*100+100, char);
273    str2 = TMPMEMALLOC(100, char);
274    sprintf(str1, "rank %d send_idx[%d] offset %d m%d c%d= (", rank, sum, len-sum, l_m, l_c);
275    for (my_i=len-sum; my_i<len; my_i++){
276      sprintf(str2, "%d ", send_idx[my_i]);
277      strcat(str1, str2);
278    }
279    fprintf(stderr, "%s)\n", str1);
280    TMPMEMFREE(str1);
281    TMPMEMFREE(str2);
282    }
283      i++;      i++;
284      }      }
285    
286    if (MY_DEBUG && rank == 0) {
287    int my_i,sum = len;
288    char *str1, *str2;
289    str1 = TMPMEMALLOC(sum*100+100, char);
290    str2 = TMPMEMALLOC(100, char);
291    sprintf(str1, "rank %d send_idx[%d] = (", rank, sum);
292    for (my_i=0; my_i<sum; my_i++){
293      sprintf(str2, "%d ", send_idx[my_i]);
294      strcat(str1, str2);
295    }
296    fprintf(stderr, "%s)\n", str1);
297    TMPMEMFREE(str1);
298    TMPMEMFREE(str2);
299    }
300    
301      /* sending */      /* sending */
302      MPI_Issend(send_offset, 2*i, MPI_INT, send->neighbor[p],      MPI_Issend(send_offset, 2*i, MPI_INT, send->neighbor[p],
303                  A->mpi_info->msg_tag_counter+rank,                  A->mpi_info->msg_tag_counter+rank,
304                  A->mpi_info->comm,                  A->mpi_info->comm,
305                  &(A->col_coupler->mpi_requests[p+recv->numNeighbors]));                  &(A->col_coupler->mpi_requests[p+recv->numNeighbors]));
306    if(MY_DEBUG) fprintf(stderr, "rank %d send %d to %d tag %d (numNeighbors %d)\n", rank, 2*i, send->neighbor[p], A->mpi_info->msg_tag_counter+rank, send->numNeighbors);
307      send_degree[p] = len;      send_degree[p] = len;
308    }    }
309    TMPMEMFREE(send_m);    TMPMEMFREE(send_m);
# Line 289  if (MY_DEBUG) fprintf(stderr, "rank %d p Line 323  if (MY_DEBUG) fprintf(stderr, "rank %d p
323                  A->col_coupler->mpi_stati);                  A->col_coupler->mpi_stati);
324    A->mpi_info->msg_tag_counter += size;    A->mpi_info->msg_tag_counter += size;
325    
326    if (MY_DEBUG) {
327    int sum = recv->offsetInShared[recv->numNeighbors] * 2;
328    char *str1, *str2;
329    str1 = TMPMEMALLOC(sum*100+100, char);
330    str2 = TMPMEMALLOC(100, char);
331    sprintf(str1, "rank %d ptr_ptr[%d] = (", rank, sum);
332    for (i=0; i<sum; i++){
333      sprintf(str2, "%d ", ptr_ptr[i]);
334      strcat(str1, str2);
335    }
336    fprintf(stderr, "%s)\n", str1);
337    TMPMEMFREE(str1);
338    TMPMEMFREE(str2);
339    }
340    
341    j = 0;    j = 0;
342    k = 0;    k = 0;
343    ptr_main[0] = 0;    ptr_main[0] = 0;
# Line 321  if (MY_DEBUG1) fprintf(stderr, "rank %d Line 370  if (MY_DEBUG1) fprintf(stderr, "rank %d
370                  A->mpi_info->msg_tag_counter+recv->neighbor[p],                  A->mpi_info->msg_tag_counter+recv->neighbor[p],
371                  A->mpi_info->comm,                  A->mpi_info->comm,
372                  &(A->col_coupler->mpi_requests[p]));                  &(A->col_coupler->mpi_requests[p]));
373    if(MY_DEBUG) fprintf(stderr, "rank %d (INDEX) recv %d from %d tag %d\n", rank, i, recv->neighbor[p], A->mpi_info->msg_tag_counter+recv->neighbor[p]);
374      } else {      } else {
375  if (MY_DEBUG) fprintf(stderr, "rank%d k %d m %d main(%d %d) couple(%d %d)\n", rank, k, m, ptr_main[m], ptr_main[k], ptr_couple[m], ptr_couple[k]);  if (MY_DEBUG1) fprintf(stderr, "rank%d k %d m %d main(%d %d) couple(%d %d)\n", rank, k, m, ptr_main[m], ptr_main[k], ptr_couple[m], ptr_couple[k]);
376      }      }
377      j += i;      j += i;
378    }    }
379    
380  if (MY_DEBUG1) fprintf(stderr, "rank %d CP1_4_3  %d %d\n", rank, recv->numNeighbors, k_ub);  if (MY_DEBUG) fprintf(stderr, "rank %d CP1_4_3  %d %d\n", rank, recv->numNeighbors, k_ub);
381    j=0;    j=0;
382    k_ub = 0;    k_ub = 0;
383    for (p=0; p<num_neighbors; p++) {    for (p=0; p<num_neighbors; p++) {
384      i = send_degree[p] - j;      i = send_degree[p] - j;
385      if (i > 0){      if (i > 0){
386      k_ub ++;      k_ub ++;
387    if (MY_DEBUG && rank == 0 && send->neighbor[p] == 1) {
388    int sum = i, my_i;
389    char *str1, *str2;
390    str1 = TMPMEMALLOC(sum*100+100, char);
391    str2 = TMPMEMALLOC(100, char);
392    sprintf(str1, "rank %d send_idx[%d] offset %d = (", rank, sum, j);
393    for (my_i=0; my_i<sum; my_i++){
394      sprintf(str2, "%d ", send_idx[j+my_i]);
395      strcat(str1, str2);
396    }
397    fprintf(stderr, "%s)\n", str1);
398    TMPMEMFREE(str1);
399    TMPMEMFREE(str2);
400    }
401          MPI_Issend(&(send_idx[j]), i, MPI_INT, send->neighbor[p],          MPI_Issend(&(send_idx[j]), i, MPI_INT, send->neighbor[p],
402                  A->mpi_info->msg_tag_counter+rank,                  A->mpi_info->msg_tag_counter+rank,
403                  A->mpi_info->comm,                  A->mpi_info->comm,
404                  &(A->col_coupler->mpi_requests[p+recv->numNeighbors]));                  &(A->col_coupler->mpi_requests[p+recv->numNeighbors]));
405    if(MY_DEBUG1) fprintf(stderr, "rank %d (INDEX) send %d to %d tag %d\n", rank, i, send->neighbor[p], A->mpi_info->msg_tag_counter+rank);
406      } else {      } else {
407  if (MY_DEBUG) fprintf(stderr, "rank%d send_degree %d, p %d, j %d\n", rank, send_degree[p], p, j);  if (MY_DEBUG1) fprintf(stderr, "rank%d send_degree %d, p %d, j %d\n", rank, send_degree[p], p, j);
408      }      }
409      j = send_degree[p];      j = send_degree[p];
410    }    }
# Line 350  if (MY_DEBUG1) fprintf(stderr, "rank %d Line 415  if (MY_DEBUG1) fprintf(stderr, "rank %d
415                  A->col_coupler->mpi_stati);                  A->col_coupler->mpi_stati);
416    A->mpi_info->msg_tag_counter += size;    A->mpi_info->msg_tag_counter += size;
417  if (MY_DEBUG1) fprintf(stderr, "rank %d CP1_5 %d %d %d\n", rank, len, ptr_main[len], ptr_couple[len]);  if (MY_DEBUG1) fprintf(stderr, "rank %d CP1_5 %d %d %d\n", rank, len, ptr_main[len], ptr_couple[len]);
418    if (MY_DEBUG && rank == 1) {
419    int my_i, sum1 = recv->offsetInShared[recv->numNeighbors];
420    int sum = ptr_main[sum1] + ptr_couple[sum1];
421    char *str1, *str2;
422    str1 = TMPMEMALLOC(sum*100+100, char);
423    str2 = TMPMEMALLOC(100, char);
424    sprintf(str1, "rank %d ptr_idx[%d] = (", rank, sum);
425    for (my_i=0; my_i<sum; my_i++){
426      sprintf(str2, "%d ", ptr_idx[my_i]);
427      strcat(str1, str2);
428    }
429    fprintf(stderr, "%s)\n", str1);
430    TMPMEMFREE(str1);
431    TMPMEMFREE(str2);
432    }
433    
434    #pragma omp parallel for private(i,j,k,m,p) schedule(static)  
435    //  #pragma omp parallel for private(i,j,k,m,p) schedule(static)
436    for (i=0; i<len; i++) {    for (i=0; i<len; i++) {
437      j = ptr_main[i];      j = ptr_main[i];
438      k = ptr_main[i+1];      k = ptr_main[i+1];
# Line 365  if (MY_DEBUG1) fprintf(stderr, "rank %d Line 446  if (MY_DEBUG1) fprintf(stderr, "rank %d
446      }      }
447    }    }
448    TMPMEMFREE(ptr_idx);    TMPMEMFREE(ptr_idx);
449  if (MY_DEBUG1) fprintf(stderr, "rank %d CP1_5_1\n", rank);  if (MY_DEBUG1) fprintf(stderr, "rank %d CP1_5_1 %d %d\n", rank, num_main_cols, len);
450    if (MY_DEBUG) {
451    int sum = num_main_cols, sum1 = ptr_main[sum];
452    char *str1, *str2;
453    str1 = TMPMEMALLOC(sum1*30+100, char);
454    str2 = TMPMEMALLOC(100, char);
455    sprintf(str1, "rank %d ptr_main[%d] = (", rank, sum);
456    for (i=0; i<sum+1; i++){
457      sprintf(str2, "%d ", ptr_main[i]);
458      strcat(str1, str2);
459    }
460    fprintf(stderr, "%s)\n", str1);
461    sprintf(str1, "rank %d idx_main[%d] = (", rank, sum1);
462    for (i=0; i<sum1; i++){
463      sprintf(str2, "%d ", idx_main[i]);
464      strcat(str1, str2);
465    }
466    fprintf(stderr, "%s)\n", str1);
467    TMPMEMFREE(str1);
468    TMPMEMFREE(str2);
469    }
470    
471    /* allocate pattern and sparsematrix for B_ext_main */    /* allocate pattern and sparsematrix for B_ext_main */
472    pattern_main = Paso_Pattern_alloc(B->col_coupleBlock->pattern->type,    pattern_main = Paso_Pattern_alloc(B->col_coupleBlock->pattern->type,
# Line 672  if (MY_DEBUG1) fprintf(stderr, "rank %d Line 773  if (MY_DEBUG1) fprintf(stderr, "rank %d
773       iptr = 0;       iptr = 0;
774       if (num_Pcouple_cols || sum > 0) {       if (num_Pcouple_cols || sum > 0) {
775      temp = TMPMEMALLOC(num_Pcouple_cols+sum, index_t);      temp = TMPMEMALLOC(num_Pcouple_cols+sum, index_t);
776      for (; iptr<sum; iptr++)      for (; iptr<sum; iptr++){
777        temp[iptr] = P->remote_coupleBlock->pattern->index[iptr];        temp[iptr] = P->remote_coupleBlock->pattern->index[iptr];
778      for (j=0; j<num_Pcouple_cols; j++, iptr++)  if (MY_DEBUG && rank ==1)
779    fprintf(stderr, "rank %d remote_block[%d] = %d\n", rank, iptr, temp[iptr]);
780        }
781        for (j=0; j<num_Pcouple_cols; j++, iptr++){
782        temp[iptr] = P->global_id[j];        temp[iptr] = P->global_id[j];
783    if (MY_DEBUG && rank ==1)
784    fprintf(stderr, "rank %d global_id[%d] = %d\n", rank, j, P->global_id[j]);
785        }
786       }       }
787       if (iptr) {       if (iptr) {
788      #ifdef USE_QSORTG      #ifdef USE_QSORTG
# Line 698  if (MY_DEBUG1) fprintf(stderr, "rank %d Line 805  if (MY_DEBUG1) fprintf(stderr, "rank %d
805      global_id_P = TMPMEMALLOC(num_Pext_cols, index_t);      global_id_P = TMPMEMALLOC(num_Pext_cols, index_t);
806      for (i=0; i<num_Pext_cols; i++)      for (i=0; i<num_Pext_cols; i++)
807        global_id_P[i] = temp[i];        global_id_P[i] = temp[i];
808    if (MY_DEBUG && rank == 1) {
809    int my_i, sum=num_Pext_cols;
810    char *str1, *str2;
811    str1 = TMPMEMALLOC(sum*100+100, char);
812    str2 = TMPMEMALLOC(100, char);
813    sprintf(str1, "rank %d global_id_P[%d] = (", rank, sum);
814    for (my_i=0; my_i<sum; my_i++) {
815      sprintf(str2, "%d ", global_id_P[my_i]);
816      strcat(str1, str2);
817    }
818    fprintf(stderr, "%s)\n", str1);
819    TMPMEMFREE(str1);
820    TMPMEMFREE(str2);
821    }
822    
823       }       }
824       if (num_Pcouple_cols || sum > 0)       if (num_Pcouple_cols || sum > 0)
825      TMPMEMFREE(temp);      TMPMEMFREE(temp);
# Line 1039  if (MY_DEBUG) fprintf(stderr, "rank %d C Line 1161  if (MY_DEBUG) fprintf(stderr, "rank %d C
1161       for (i=0, iptr=0; i<sum; i++) {       for (i=0, iptr=0; i<sum; i++) {
1162      if (RAP_ext_idx[i] < offset || RAP_ext_idx[i] >= j1_ub)  /* XXX */ {      if (RAP_ext_idx[i] < offset || RAP_ext_idx[i] >= j1_ub)  /* XXX */ {
1163        temp[iptr++] = RAP_ext_idx[i];          /* XXX */        temp[iptr++] = RAP_ext_idx[i];          /* XXX */
1164  if (MY_DEBUG && rank ==0) fprintf(stderr, "RAP_ext_idx[%d]=%d\n", i, RAP_ext_idx[i]);  if (MY_DEBUG && rank ==1) fprintf(stderr, "RAP_ext_idx[%d]=%d\n", i, RAP_ext_idx[i]);
1165  }  }
1166       }       }
1167       for (j=0; j<num_Pext_cols; j++, iptr++)       for (j=0; j<num_Pext_cols; j++, iptr++){
1168      temp[iptr] = global_id_P[j];      temp[iptr] = global_id_P[j];
1169         }
1170    
1171       if (iptr) {       if (iptr) {
1172      #ifdef USE_QSORTG      #ifdef USE_QSORTG
# Line 1546  if (MY_DEBUG){ Line 1669  if (MY_DEBUG){
1669  if (rank == 0) fprintf(stderr, "A[%d]=%f A[%d]=%f A[%d]=%f A[%d]=%f \n", RAP_main_idx[0], RAP_main_val[0], RAP_main_idx[1], RAP_main_val[1], RAP_main_idx[2], RAP_main_val[2], RAP_main_idx[3], RAP_main_val[3]);  if (rank == 0) fprintf(stderr, "A[%d]=%f A[%d]=%f A[%d]=%f A[%d]=%f \n", RAP_main_idx[0], RAP_main_val[0], RAP_main_idx[1], RAP_main_val[1], RAP_main_idx[2], RAP_main_val[2], RAP_main_idx[3], RAP_main_val[3]);
1670  }  }
1671  if (MY_DEBUG1) fprintf(stderr, "rank %d CP10\n", rank);  if (MY_DEBUG1) fprintf(stderr, "rank %d CP10\n", rank);
1672    if (MY_DEBUG1 && rank == 1) {
1673    int sum = num_RAPext_cols, my_i;
1674    char *str1, *str2;
1675    str1 = TMPMEMALLOC(sum*100+100, char);
1676    str2 = TMPMEMALLOC(100, char);
1677    sprintf(str1, "rank %d global_id_RAP[%d] = (", rank, sum);
1678    for (my_i=0; my_i<sum; my_i++){
1679      sprintf(str2, "%d ", global_id_RAP[my_i]);
1680      strcat(str1, str2);
1681    }
1682    fprintf(stderr, "%s)\n", str1);
1683    sprintf(str1, "rank %d distribution = (", rank);
1684    for (my_i=0; my_i<size; my_i++){
1685      sprintf(str2, "%d ", P->pattern->input_distribution->first_component[my_i+1]);
1686      strcat(str1, str2);
1687    }
1688    fprintf(stderr, "%s)\n", str1);
1689    TMPMEMFREE(str1);
1690    TMPMEMFREE(str2);
1691    }
1692    
1693    
1694     /* Check whether there are empty columns in RAP_couple */     /* Check whether there are empty columns in RAP_couple */
1695     #pragma omp parallel for schedule(static) private(i)     #pragma omp parallel for schedule(static) private(i)
# Line 1599  if (MY_DEBUG) fprintf(stderr, "rank %d C Line 1743  if (MY_DEBUG) fprintf(stderr, "rank %d C
1743     num_neighbors = 0;     num_neighbors = 0;
1744     offsetInShared[0] = 0;     offsetInShared[0] = 0;
1745     for (i=0, j=0, k=dist[j+1]; i<num_RAPext_cols; i++) {     for (i=0, j=0, k=dist[j+1]; i<num_RAPext_cols; i++) {
1746  if (MY_DEBUG && rank == 0) fprintf(stderr, "i%d j%d k%d %d %d\n", i, j, k, global_id_RAP[i], recv_len[j]);  if (MY_DEBUG && rank == 1) fprintf(stderr, "i%d j%d k%d %d %d\n", i, j, k, global_id_RAP[i], recv_len[j]);
1747       shared[i] = i + num_Pmain_cols;       shared[i] = i + num_Pmain_cols;
1748       if (k <= global_id_RAP[i]) {       if (k <= global_id_RAP[i]) {
1749      if (recv_len[j] > 0) {      if (recv_len[j] > 0) {
# Line 1709  if (MY_DEBUG) fprintf(stderr, "rank %d C Line 1853  if (MY_DEBUG) fprintf(stderr, "rank %d C
1853        for (p=0; p<recv->numNeighbors; p++) {        for (p=0; p<recv->numNeighbors; p++) {
1854          if (i_c < recv->offsetInShared[p+1]) {          if (i_c < recv->offsetInShared[p+1]) {
1855            k = recv->neighbor[p];            k = recv->neighbor[p];
1856    if (MY_DEBUG1 && rank == 1 && k == 1)
1857    fprintf(stderr, "******* i_r %d i_c %d p %d ub %d numNeighbors %d\n", i_r, i_c, p, recv->offsetInShared[p+1], recv->numNeighbors);
1858            if (send_ptr[k][i_r] == 0) sum++;            if (send_ptr[k][i_r] == 0) sum++;
1859            send_ptr[k][i_r] ++;            send_ptr[k][i_r] ++;
1860            send_idx[k][len[k]] = global_id_RAP[i_c];            send_idx[k][len[k]] = global_id_RAP[i_c];
# Line 1733  if (MY_DEBUG) fprintf(stderr, "rank %d C Line 1879  if (MY_DEBUG) fprintf(stderr, "rank %d C
1879     offsetInShared[0] = 0;     offsetInShared[0] = 0;
1880  if (MY_DEBUG) fprintf(stderr, "rank %d CP14_5 %d\n", rank, size);  if (MY_DEBUG) fprintf(stderr, "rank %d CP14_5 %d\n", rank, size);
1881     for (p=0, k=0; p<size; p++) {     for (p=0, k=0; p<size; p++) {
1882  if (MY_DEBUG) fprintf(stderr, "rank %d CP14_6 %d\n", rank, p);  if (MY_DEBUG && rank == 1) fprintf(stderr, "rank %d CP14_6 %d %d %d\n", rank, p, k, num_neighbors);
1883       for (i=0; i<num_Pmain_cols; i++) {       for (i=0; i<num_Pmain_cols; i++) {
1884      if (send_ptr[p][i] > 0) {      if (send_ptr[p][i] > 0) {
1885        shared[k] = i;        shared[k] = i;
# Line 1794  if (MY_DEBUG1) fprintf(stderr, "rank %d Line 1940  if (MY_DEBUG1) fprintf(stderr, "rank %d
1940       MPI_Irecv(&(row_couple_ptr[j]), i-j, MPI_INT, neighbor[p],       MPI_Irecv(&(row_couple_ptr[j]), i-j, MPI_INT, neighbor[p],
1941          mpi_info->msg_tag_counter+neighbor[p],          mpi_info->msg_tag_counter+neighbor[p],
1942          mpi_info->comm, &mpi_requests[p]);          mpi_info->comm, &mpi_requests[p]);
1943    if (MY_DEBUG) fprintf(stderr, "rank %d RECV %d from %d tag %d\n", rank, i-j, neighbor[p], mpi_info->msg_tag_counter+neighbor[p]);
1944     }     }
1945     send = row_connector->send;     send = row_connector->send;
1946     for (p=0; p<send->numNeighbors; p++) {     for (p=0; p<send->numNeighbors; p++) {
# Line 1801  if (MY_DEBUG1) fprintf(stderr, "rank %d Line 1948  if (MY_DEBUG1) fprintf(stderr, "rank %d
1948          MPI_INT, send->neighbor[p],          MPI_INT, send->neighbor[p],
1949          mpi_info->msg_tag_counter+rank,          mpi_info->msg_tag_counter+rank,
1950          mpi_info->comm, &mpi_requests[p+num_neighbors]);          mpi_info->comm, &mpi_requests[p+num_neighbors]);
1951    if (MY_DEBUG) fprintf(stderr, "rank %d SEND %d TO %d tag %d\n", rank, send_len[send->neighbor[p]], send->neighbor[p], mpi_info->msg_tag_counter+rank);
1952     }     }
1953     MPI_Waitall(num_neighbors + send->numNeighbors,     MPI_Waitall(num_neighbors + send->numNeighbors,
1954      mpi_requests, mpi_stati);      mpi_requests, mpi_stati);
1955     mpi_info->msg_tag_counter += size;     mpi_info->msg_tag_counter += size;
1956     TMPMEMFREE(send_len);     TMPMEMFREE(send_len);
1957  if (MY_DEBUG) fprintf(stderr, "rank %d CP17\n", rank);  if (MY_DEBUG1) fprintf(stderr, "rank %d CP17\n", rank);
1958    
1959     sum = 0;     sum = 0;
1960     for (i=0; i<num_RAPext_rows; i++) {     for (i=0; i<num_RAPext_rows; i++) {
# Line 1854  if (MY_DEBUG) fprintf(stderr, "rank %d C Line 2002  if (MY_DEBUG) fprintf(stderr, "rank %d C
2002     TMPMEMFREE(mpi_requests);     TMPMEMFREE(mpi_requests);
2003     TMPMEMFREE(mpi_stati);     TMPMEMFREE(mpi_stati);
2004     Esys_MPIInfo_free(mpi_info);     Esys_MPIInfo_free(mpi_info);
2005  if (MY_DEBUG) fprintf(stderr, "rank %d CP18\n", rank);  if (MY_DEBUG1) fprintf(stderr, "rank %d CP18\n", rank);
2006    
2007     /* Now, we can create pattern for mainBlock and coupleBlock */     /* Now, we can create pattern for mainBlock and coupleBlock */
2008     main_pattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, num_Pmain_cols,     main_pattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, num_Pmain_cols,

Legend:
Removed from v.3763  
changed lines
  Added in v.3767

  ViewVC Help
Powered by ViewVC 1.1.26