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

Diff of /trunk/paso/src/AMG_Prolongation.c

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

revision 3827 by lgao, Tue Feb 14 11:42:08 2012 UTC revision 3828 by lgao, Wed Feb 15 03:27:58 2012 UTC
# Line 1  Line 1 
 //  
1  /*******************************************************  /*******************************************************
2  *  *
3  * Copyright (c) 2003-2010 by University of Queensland  * Copyright (c) 2003-2010 by University of Queensland
# Line 52  Paso_SystemMatrix* Paso_Preconditioner_A Line 51  Paso_SystemMatrix* Paso_Preconditioner_A
51                                 const dim_t n_C, index_t* counter_C, const index_t interpolation_method)                                 const dim_t n_C, index_t* counter_C, const index_t interpolation_method)
52  {  {
53     Esys_MPIInfo *mpi_info=Esys_MPIInfo_getReference(A_p->mpi_info);     Esys_MPIInfo *mpi_info=Esys_MPIInfo_getReference(A_p->mpi_info);
    Paso_SparseMatrix *main_block=NULL, *couple_block=NULL;  
54     Paso_SystemMatrix *out=NULL;     Paso_SystemMatrix *out=NULL;
55     Paso_SystemMatrixPattern *pattern=NULL;     Paso_SystemMatrixPattern *pattern=NULL;
56     Paso_Distribution *input_dist=NULL, *output_dist=NULL;     Paso_Distribution *input_dist=NULL, *output_dist=NULL;
57     Paso_SharedComponents *send =NULL, *recv=NULL;     Paso_SharedComponents *send =NULL, *recv=NULL;
58     Paso_Connector *col_connector=NULL, *row_connector=NULL;     Paso_Connector *col_connector=NULL;
    Paso_Coupler *coupler=NULL;  
59     Paso_Pattern *main_pattern=NULL, *couple_pattern=NULL;     Paso_Pattern *main_pattern=NULL, *couple_pattern=NULL;
60     const dim_t row_block_size=A_p->row_block_size;     const dim_t row_block_size=A_p->row_block_size;
61     const dim_t col_block_size=A_p->col_block_size;     const dim_t col_block_size=A_p->col_block_size;
62     const dim_t my_n=A_p->mainBlock->numCols;     const dim_t my_n=A_p->mainBlock->numCols;
63     const dim_t overlap_n=A_p->col_coupleBlock->numCols;     const dim_t overlap_n=A_p->col_coupleBlock->numCols;
    const dim_t n = my_n + overlap_n;  
64     const dim_t num_threads=omp_get_max_threads();     const dim_t num_threads=omp_get_max_threads();
65     index_t size=mpi_info->size, rank=mpi_info->rank, *dist=NULL;     index_t size=mpi_info->size, *dist=NULL;
66     index_t *main_p=NULL, *couple_p=NULL, *main_idx=NULL, *couple_idx=NULL;     index_t *main_p=NULL, *couple_p=NULL, *main_idx=NULL, *couple_idx=NULL;
67     index_t *shared=NULL, *offsetInShared=NULL;     index_t *shared=NULL, *offsetInShared=NULL;
68     index_t *recv_shared=NULL, *send_shared=NULL;     index_t *recv_shared=NULL, *send_shared=NULL;
69     index_t sum, j, iptr;;     index_t sum, i, j, k, l, p, q, iptr;
70     dim_t i, my_n_C, k, l, p, q, global_label, num_neighbors;     index_t my_n_C, global_label, num_neighbors;
71       #ifdef ESYS_MPI
72       index_t rank=mpi_info->rank;
73       #endif
74     Esys_MPI_rank *neighbor=NULL;     Esys_MPI_rank *neighbor=NULL;
75     #ifdef ESYS_MPI     #ifdef ESYS_MPI
76       MPI_Request* mpi_requests=NULL;       MPI_Request* mpi_requests=NULL;
# Line 82  Paso_SystemMatrix* Paso_Preconditioner_A Line 81  Paso_SystemMatrix* Paso_Preconditioner_A
81    
82     /* number of C points in current distribution */     /* number of C points in current distribution */
83     my_n_C = 0;     my_n_C = 0;
84       sum=0;
85     if (num_threads>1) {     if (num_threads>1) {
86       #pragma omp parallel private(i,sum)       #pragma omp parallel private(i,sum)
87       {       {
# Line 113  Paso_SystemMatrix* Paso_Preconditioner_A Line 113  Paso_SystemMatrix* Paso_Preconditioner_A
113     output_dist=Paso_Distribution_alloc(mpi_info, dist, 1, 0);     output_dist=Paso_Distribution_alloc(mpi_info, dist, 1, 0);
114     dist = TMPMEMALLOC(size+1, index_t); /* now prepare for col distribution */     dist = TMPMEMALLOC(size+1, index_t); /* now prepare for col distribution */
115     Esys_checkPtr(dist);     Esys_checkPtr(dist);
116       #ifdef ESYS_MPI
117     MPI_Allgather(&my_n_C, 1, MPI_INT, dist, 1, MPI_INT, mpi_info->comm);     MPI_Allgather(&my_n_C, 1, MPI_INT, dist, 1, MPI_INT, mpi_info->comm);
118       #endif
119     global_label=0;     global_label=0;
120     for (i=0; i<size; i++) {     for (i=0; i<size; i++) {
121       k = dist[i];       k = dist[i];
# Line 231  Paso_SystemMatrix* Paso_Preconditioner_A Line 233  Paso_SystemMatrix* Paso_Preconditioner_A
233    
234     for (p=0; p<send->numNeighbors; p++) {     for (p=0; p<send->numNeighbors; p++) {
235       i = send->offsetInShared[p];       i = send->offsetInShared[p];
236         #ifdef ESYS_MPI
237       MPI_Irecv (&(send_shared[i]), send->offsetInShared[p+1]-i, MPI_INT,       MPI_Irecv (&(send_shared[i]), send->offsetInShared[p+1]-i, MPI_INT,
238          send->neighbor[p], mpi_info->msg_tag_counter+send->neighbor[p],          send->neighbor[p], mpi_info->msg_tag_counter+send->neighbor[p],
239          mpi_info->comm, &mpi_requests[p]);          mpi_info->comm, &mpi_requests[p]);
240         #endif
241     }     }
242    
243     num_neighbors = 0;     num_neighbors = 0;
# Line 257  Paso_SystemMatrix* Paso_Preconditioner_A Line 261  Paso_SystemMatrix* Paso_Preconditioner_A
261      num_neighbors++;      num_neighbors++;
262      offsetInShared[num_neighbors] = q;      offsetInShared[num_neighbors] = q;
263       }       }
264         #ifdef ESYS_MPI
265       MPI_Issend(&(recv_shared[recv->offsetInShared[i]]),       MPI_Issend(&(recv_shared[recv->offsetInShared[i]]),
266          k-recv->offsetInShared[i], MPI_INT, recv->neighbor[i],          k-recv->offsetInShared[i], MPI_INT, recv->neighbor[i],
267          mpi_info->msg_tag_counter+rank, mpi_info->comm,          mpi_info->msg_tag_counter+rank, mpi_info->comm,
268          &mpi_requests[i+send->numNeighbors]);          &mpi_requests[i+send->numNeighbors]);
269         #endif
270     }     }
271     recv = Paso_SharedComponents_alloc(my_n_C, num_neighbors, neighbor, shared,     recv = Paso_SharedComponents_alloc(my_n_C, num_neighbors, neighbor, shared,
272                                        offsetInShared, 1, 0, mpi_info);                                        offsetInShared, 1, 0, mpi_info);
273    
274     /* now we can build the sender */     /* now we can build the sender */
275       #ifdef ESYS_MPI
276     MPI_Waitall(recv->numNeighbors+send->numNeighbors, mpi_requests, mpi_stati);     MPI_Waitall(recv->numNeighbors+send->numNeighbors, mpi_requests, mpi_stati);
277       #endif
278     mpi_info->msg_tag_counter += size;     mpi_info->msg_tag_counter += size;
279     TMPMEMFREE(mpi_requests);     TMPMEMFREE(mpi_requests);
280     TMPMEMFREE(mpi_stati);     TMPMEMFREE(mpi_stati);
# Line 368  Paso_SystemMatrix* Paso_Preconditioner_A Line 376  Paso_SystemMatrix* Paso_Preconditioner_A
376  void Paso_Preconditioner_AMG_setDirectProlongation(Paso_SystemMatrix* P,  void Paso_Preconditioner_AMG_setDirectProlongation(Paso_SystemMatrix* P,
377      Paso_SystemMatrix* A, const index_t* offset_S, const dim_t* degree_S,      Paso_SystemMatrix* A, const index_t* offset_S, const dim_t* degree_S,
378          const index_t* S, const index_t *counter_C) {          const index_t* S, const index_t *counter_C) {
    Esys_MPIInfo *mpi_info=A->mpi_info;  
379     Paso_SparseMatrix *main_block=P->mainBlock;     Paso_SparseMatrix *main_block=P->mainBlock;
380     Paso_SparseMatrix *couple_block=P->col_coupleBlock;     Paso_SparseMatrix *couple_block=P->col_coupleBlock;
381     Paso_Pattern *main_pattern=main_block->pattern;     Paso_Pattern *main_pattern=main_block->pattern;
382     Paso_Pattern *couple_pattern=couple_block->pattern;     Paso_Pattern *couple_pattern=couple_block->pattern;
    const dim_t row_block_size=A->row_block_size;  
    const dim_t col_block_size=A->col_block_size;  
383     const dim_t my_n=A->mainBlock->numRows;     const dim_t my_n=A->mainBlock->numRows;
384     const dim_t overlap_n=A->row_coupleBlock->numRows;     index_t range;
    const dim_t n = my_n + overlap_n;  
    const dim_t num_threads=omp_get_max_threads();  
    index_t size=mpi_info->size, rank=mpi_info->rank, range;  
385    
386     dim_t i;     dim_t i;
387     register double alpha, beta, sum_all_neg, sum_all_pos, sum_strong_neg, sum_strong_pos, A_ij, A_ii, rtmp;     register double alpha, beta, sum_all_neg, sum_all_pos, sum_strong_neg, sum_strong_pos, A_ij, A_ii, rtmp;
# Line 510  void Paso_Preconditioner_AMG_setDirectPr Line 512  void Paso_Preconditioner_AMG_setDirectPr
512  void Paso_Preconditioner_AMG_setDirectProlongation_Block(Paso_SystemMatrix* P,  void Paso_Preconditioner_AMG_setDirectProlongation_Block(Paso_SystemMatrix* P,
513          Paso_SystemMatrix* A, const index_t* offset_S, const dim_t* degree_S,          Paso_SystemMatrix* A, const index_t* offset_S, const dim_t* degree_S,
514          const index_t* S, const index_t *counter_C) {          const index_t* S, const index_t *counter_C) {
    Esys_MPIInfo *mpi_info=A->mpi_info;  
515     Paso_SparseMatrix *main_block=P->mainBlock;     Paso_SparseMatrix *main_block=P->mainBlock;
516     Paso_SparseMatrix *couple_block=P->col_coupleBlock;     Paso_SparseMatrix *couple_block=P->col_coupleBlock;
517     Paso_Pattern *main_pattern=main_block->pattern;     Paso_Pattern *main_pattern=main_block->pattern;
518     Paso_Pattern *couple_pattern=couple_block->pattern;     Paso_Pattern *couple_pattern=couple_block->pattern;
519     const dim_t row_block_size=A->row_block_size;     const dim_t row_block_size=A->row_block_size;
    const dim_t col_block_size=A->col_block_size;  
520     const dim_t A_block = A->block_size;     const dim_t A_block = A->block_size;
521     const dim_t my_n=A->mainBlock->numRows;     const dim_t my_n=A->mainBlock->numRows;
522     const dim_t overlap_n=A->row_coupleBlock->numRows;     index_t range;
    const dim_t n = my_n + overlap_n;  
    const dim_t num_threads=omp_get_max_threads();  
    index_t size=mpi_info->size, rank=mpi_info->rank, range;  
523    
524     dim_t i;     dim_t i;
525     double *alpha, *beta, *sum_all_neg, *sum_all_pos, *sum_strong_neg, *sum_strong_pos, *A_ii;     double *alpha, *beta, *sum_all_neg, *sum_all_pos, *sum_strong_neg, *sum_strong_pos, *A_ii;
# Line 709  void Paso_Preconditioner_AMG_setDirectPr Line 706  void Paso_Preconditioner_AMG_setDirectPr
706  void Paso_Preconditioner_AMG_setClassicProlongation(Paso_SystemMatrix* P,  void Paso_Preconditioner_AMG_setClassicProlongation(Paso_SystemMatrix* P,
707      Paso_SystemMatrix* A, const index_t* offset_S, const dim_t* degree_S,      Paso_SystemMatrix* A, const index_t* offset_S, const dim_t* degree_S,
708      const index_t* S, const index_t *counter_C) {      const index_t* S, const index_t *counter_C) {
    Esys_MPIInfo *mpi_info=A->mpi_info;  
709     Paso_SparseMatrix *main_block=P->mainBlock;     Paso_SparseMatrix *main_block=P->mainBlock;
710     Paso_SparseMatrix *couple_block=P->col_coupleBlock;     Paso_SparseMatrix *couple_block=P->col_coupleBlock;
711     Paso_Pattern *main_pattern=main_block->pattern;     Paso_Pattern *main_pattern=main_block->pattern;
712     Paso_Pattern *couple_pattern=couple_block->pattern;     Paso_Pattern *couple_pattern=couple_block->pattern;
    const dim_t row_block_size=A->row_block_size;  
    const dim_t col_block_size=A->col_block_size;  
713     const dim_t my_n=A->mainBlock->numRows;     const dim_t my_n=A->mainBlock->numRows;
714     const dim_t overlap_n=A->row_coupleBlock->numRows;     index_t range, range_j;
    const dim_t n = my_n + overlap_n;  
    const dim_t num_threads=omp_get_max_threads();  
    index_t size=mpi_info->size, rank=mpi_info->rank, range, range_j;  
715    
716     dim_t i, q;     dim_t i, q;
717     double *D_s=NULL;     double *D_s=NULL;
# Line 944  void Paso_Preconditioner_AMG_setClassicP Line 935  void Paso_Preconditioner_AMG_setClassicP
935  void Paso_Preconditioner_AMG_setClassicProlongation_Block(  void Paso_Preconditioner_AMG_setClassicProlongation_Block(
936      Paso_SystemMatrix* P, Paso_SystemMatrix* A, const index_t* offset_S,      Paso_SystemMatrix* P, Paso_SystemMatrix* A, const index_t* offset_S,
937      const dim_t* degree_S, const index_t* S, const index_t *counter_C) {      const dim_t* degree_S, const index_t* S, const index_t *counter_C) {
    Esys_MPIInfo *mpi_info=A->mpi_info;  
938     Paso_SparseMatrix *main_block=P->mainBlock;     Paso_SparseMatrix *main_block=P->mainBlock;
939     Paso_SparseMatrix *couple_block=P->col_coupleBlock;     Paso_SparseMatrix *couple_block=P->col_coupleBlock;
940     Paso_Pattern *main_pattern=main_block->pattern;     Paso_Pattern *main_pattern=main_block->pattern;
941     Paso_Pattern *couple_pattern=couple_block->pattern;     Paso_Pattern *couple_pattern=couple_block->pattern;
942     const dim_t row_block=A->row_block_size;     const dim_t row_block=A->row_block_size;
    const dim_t col_block=A->col_block_size;  
943     const dim_t my_n=A->mainBlock->numRows;     const dim_t my_n=A->mainBlock->numRows;
    const dim_t overlap_n=A->row_coupleBlock->numRows;  
    const dim_t n = my_n + overlap_n;  
944     const dim_t A_block = A->block_size;     const dim_t A_block = A->block_size;
945     const dim_t num_threads=omp_get_max_threads();     index_t range, range_j;
    index_t size=mpi_info->size, rank=mpi_info->rank, range, range_j;  
946    
947     dim_t i, q, ib;     dim_t i, q, ib;
948     double *D_s=NULL;     double *D_s=NULL;
# Line 1112  void Paso_Preconditioner_AMG_setClassicP Line 1098  void Paso_Preconditioner_AMG_setClassicP
1098                     /* first, the row_coupleBlock part */                     /* first, the row_coupleBlock part */
1099                     range_j = A->row_coupleBlock->pattern->ptr[j + 1];                     range_j = A->row_coupleBlock->pattern->ptr[j + 1];
1100                             for (iPtr_j=A->row_coupleBlock->pattern->ptr[j];iPtr_j<range_j; iPtr_j++) {                             for (iPtr_j=A->row_coupleBlock->pattern->ptr[j];iPtr_j<range_j; iPtr_j++) {
 //                              const double* A_jm=&(A->row_coupleBlock->val[iPtr_j*A_block]);  
 //                              const index_t m=A->row_coupleBlock->pattern->index[iPtr_j];  
1101                                      /* is m an interpolation point ? */                                      /* is m an interpolation point ? */
1102                      index_t m, *where_p_m;                      index_t m, *where_p_m;
1103                      double *A_jm;                      double *A_jm;

Legend:
Removed from v.3827  
changed lines
  Added in v.3828

  ViewVC Help
Powered by ViewVC 1.1.26