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

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

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

revision 4828 by caltinay, Tue Apr 1 03:50:23 2014 UTC revision 4829 by caltinay, Thu Apr 3 04:02:53 2014 UTC
# Line 42  Paso_MergedSolver* Paso_MergedSolver_all Line 42  Paso_MergedSolver* Paso_MergedSolver_all
42     const dim_t n_block = A->mainBlock->row_block_size;     const dim_t n_block = A->mainBlock->row_block_size;
43     const dim_t* dist = A->pattern->input_distribution->first_component;     const dim_t* dist = A->pattern->input_distribution->first_component;
44     dim_t i;     dim_t i;
    paso::SparseMatrix *A_temp;  
45    
46     Paso_MergedSolver* out = NULL;     Paso_MergedSolver* out = NULL;
47     A_temp = Paso_MergedSolver_mergeSystemMatrix(A);     paso::SparseMatrix_ptr A_temp(Paso_MergedSolver_mergeSystemMatrix(A));
48     out = new Paso_MergedSolver;     out = new Paso_MergedSolver;
49     if (Esys_noError()) {     if (Esys_noError()) {
     
50         /* merge matrix */         /* merge matrix */
        out->A = NULL;  
51         out->mpi_info=Esys_MPIInfo_getReference(A->mpi_info);         out->mpi_info=Esys_MPIInfo_getReference(A->mpi_info);
52         out->reordering = options->reordering;         out->reordering = options->reordering;
53         out->refinements = options->coarse_matrix_refinements;         out->refinements = options->coarse_matrix_refinements;
# Line 73  Paso_MergedSolver* Paso_MergedSolver_all Line 70  Paso_MergedSolver* Paso_MergedSolver_all
70    
71        if (rank == 0) {        if (rank == 0) {
72           /* solve locally */           /* solve locally */
73           #ifdef MKL  #ifdef MKL
74             out->A = paso::SparseMatrix_unroll(MATRIX_FORMAT_BLK1 + MATRIX_FORMAT_OFFSET1, A_temp);             out->A = A_temp->unroll(MATRIX_FORMAT_BLK1 + MATRIX_FORMAT_OFFSET1);
75             out->A->solver_package = PASO_MKL;             out->A->solver_package = PASO_MKL;
76           #else  #elif defined UMFPACK
77             #ifdef UMFPACK             out->A = A_temp->unroll(MATRIX_FORMAT_BLK1 + MATRIX_FORMAT_CSC);
78               out->A  = paso::SparseMatrix_unroll(MATRIX_FORMAT_BLK1 + MATRIX_FORMAT_CSC, A_temp);             out->A->solver_package = PASO_UMFPACK;
79               out->A->solver_package = PASO_UMFPACK;  #else
80             #else             out->A->solver_p = Paso_Preconditioner_LocalSmoother_alloc(out->A, (options->smoother == PASO_JACOBI), out->verbose);
81               out->A->solver_p = Paso_Preconditioner_LocalSmoother_alloc(out->A, (options->smoother == PASO_JACOBI), out->verbose);             out->A->solver_package = PASO_SMOOTHER;
82               out->A->solver_package = PASO_SMOOTHER;  #endif
            #endif  
          #endif  
83         }         }
84     }     }
85    
    paso::SparseMatrix_free(A_temp);  
86     if ( Esys_noError()) {     if ( Esys_noError()) {
87        return out;        return out;
88     } else {     } else {
# Line 97  Paso_MergedSolver* Paso_MergedSolver_all Line 91  Paso_MergedSolver* Paso_MergedSolver_all
91     }     }
92  }  }
93    
94  void Paso_MergedSolver_free(Paso_MergedSolver* in) {  void Paso_MergedSolver_free(Paso_MergedSolver* in)
95       if (in!=NULL) {  {
96           paso::SparseMatrix_free(in->A);      if (in!=NULL) {
97      delete[] in->x;          delete[] in->x;
98      delete[] in->b;          delete[] in->b;
99      delete[] in->counts;          delete[] in->counts;
100      delete[] in->offset;          delete[] in->offset;
101      delete in;          delete in;
102       }      }
103  }  }
104    
105  /* Merge the system matrix which is distributed on ranks into a complete  /* Merge the system matrix which is distributed on ranks into a complete
106     matrix on rank 0, then solve this matrix on rank 0 only */     matrix on rank 0, then solve this matrix on rank 0 only */
107  paso::SparseMatrix* Paso_MergedSolver_mergeSystemMatrix(Paso_SystemMatrix* A) {  paso::SparseMatrix_ptr Paso_MergedSolver_mergeSystemMatrix(Paso_SystemMatrix* A)
108    {
109    index_t i, iptr, j, n, remote_n, global_n, len, offset, tag;    index_t i, iptr, j, n, remote_n, global_n, len, offset, tag;
110    index_t row_block_size, col_block_size, block_size;    index_t row_block_size, col_block_size, block_size;
111    index_t size=A->mpi_info->size;    index_t size=A->mpi_info->size;
# Line 118  paso::SparseMatrix* Paso_MergedSolver_me Line 113  paso::SparseMatrix* Paso_MergedSolver_me
113    index_t *ptr=NULL, *idx=NULL, *ptr_global=NULL, *idx_global=NULL;    index_t *ptr=NULL, *idx=NULL, *ptr_global=NULL, *idx_global=NULL;
114    index_t *temp_n=NULL, *temp_len=NULL;    index_t *temp_n=NULL, *temp_len=NULL;
115    double  *val=NULL;    double  *val=NULL;
116    paso::SparseMatrix *out=NULL;    paso::SparseMatrix_ptr out;
117    #ifdef ESYS_MPI    #ifdef ESYS_MPI
118      MPI_Request* mpi_requests=NULL;      MPI_Request* mpi_requests=NULL;
119      MPI_Status* mpi_stati=NULL;      MPI_Status* mpi_stati=NULL;
# Line 131  paso::SparseMatrix* Paso_MergedSolver_me Line 126  paso::SparseMatrix* Paso_MergedSolver_me
126      ptr = new index_t[n];      ptr = new index_t[n];
127      #pragma omp parallel for private(i)      #pragma omp parallel for private(i)
128      for (i=0; i<n; i++) ptr[i] = i;      for (i=0; i<n; i++) ptr[i] = i;
129      out = paso::SparseMatrix_getSubmatrix(A->mainBlock, n, n, ptr, ptr);      out = A->mainBlock->getSubmatrix(n, n, ptr, ptr);
130      delete[] ptr;      delete[] ptr;
131      return out;      return out;
132    }    }
# Line 169  paso::SparseMatrix* Paso_MergedSolver_me Line 164  paso::SparseMatrix* Paso_MergedSolver_me
164      /* Second, receive ptr values from other ranks */      /* Second, receive ptr values from other ranks */
165      for (i=1; i<size; i++) {      for (i=1; i<size; i++) {
166        remote_n = A->row_distribution->first_component[i+1] -        remote_n = A->row_distribution->first_component[i+1] -
167           A->row_distribution->first_component[i];                   A->row_distribution->first_component[i];
168        #ifdef ESYS_MPI        #ifdef ESYS_MPI
169        MPI_Irecv(&(ptr_global[iptr]), remote_n, MPI_INT, i,        MPI_Irecv(&(ptr_global[iptr]), remote_n, MPI_INT, i,
170              A->mpi_info->msg_tag_counter+i,                          A->mpi_info->msg_tag_counter+i,
171              A->mpi_info->comm,                          A->mpi_info->comm,
172              &mpi_requests[i]);                          &mpi_requests[i]);
173        #endif        #endif
174        temp_n[i] = remote_n;        temp_n[i] = remote_n;
175        iptr += remote_n;        iptr += remote_n;
# Line 189  paso::SparseMatrix* Paso_MergedSolver_me Line 184  paso::SparseMatrix* Paso_MergedSolver_me
184      offset = -1;      offset = -1;
185      for (i=0; i<size; i++) {      for (i=0; i<size; i++) {
186        if (temp_n[i] > 0) {        if (temp_n[i] > 0) {
187      offset += temp_n[i];          offset += temp_n[i];
188      len += ptr_global[offset];          len += ptr_global[offset];
189      temp_len[i] = ptr_global[offset];          temp_len[i] = ptr_global[offset];
190        }else        }else
191      temp_len[i] = 0;          temp_len[i] = 0;
192      }      }
193    
194      idx_global = new index_t[len];      idx_global = new index_t[len];
# Line 203  paso::SparseMatrix* Paso_MergedSolver_me Line 198  paso::SparseMatrix* Paso_MergedSolver_me
198        len = temp_len[i];        len = temp_len[i];
199        #ifdef ESYS_MPI        #ifdef ESYS_MPI
200        MPI_Irecv(&(idx_global[iptr]), len, MPI_INT, i,        MPI_Irecv(&(idx_global[iptr]), len, MPI_INT, i,
201              A->mpi_info->msg_tag_counter+i,                          A->mpi_info->msg_tag_counter+i,
202              A->mpi_info->comm,                          A->mpi_info->comm,
203              &mpi_requests[i]);                          &mpi_requests[i]);
204        #endif        #endif
205        remote_n = temp_n[i];        remote_n = temp_n[i];
206        for (j=0; j<remote_n; j++) {        for (j=0; j<remote_n; j++) {
207      ptr_global[j+offset] = ptr_global[j+offset] + iptr;          ptr_global[j+offset] = ptr_global[j+offset] + iptr;
208        }        }
209        offset += remote_n;        offset += remote_n;
210        iptr += len;        iptr += len;
# Line 227  paso::SparseMatrix* Paso_MergedSolver_me Line 222  paso::SparseMatrix* Paso_MergedSolver_me
222      /* Then generate the sparse matrix */      /* Then generate the sparse matrix */
223      paso::Pattern_ptr pattern(new paso::Pattern(A->mainBlock->pattern->type,      paso::Pattern_ptr pattern(new paso::Pattern(A->mainBlock->pattern->type,
224                  global_n, global_n, ptr_global, idx_global));                  global_n, global_n, ptr_global, idx_global));
225      out = paso::SparseMatrix_alloc(A->mainBlock->type, pattern,      out.reset(new paso::SparseMatrix(A->mainBlock->type, pattern,
226              row_block_size, col_block_size, FALSE);                          row_block_size, col_block_size, false));
227    
228      /* Finally, receive and copy the value */      /* Finally, receive and copy the value */
229      iptr = temp_len[0] * block_size;      iptr = temp_len[0] * block_size;
# Line 255  paso::SparseMatrix* Paso_MergedSolver_me Line 250  paso::SparseMatrix* Paso_MergedSolver_me
250      tag = A->mpi_info->msg_tag_counter+rank;      tag = A->mpi_info->msg_tag_counter+rank;
251      #ifdef ESYS_MPI      #ifdef ESYS_MPI
252      MPI_Issend(&(ptr[1]), n, MPI_INT, 0, tag, A->mpi_info->comm,      MPI_Issend(&(ptr[1]), n, MPI_INT, 0, tag, A->mpi_info->comm,
253              &mpi_requests[0]);                          &mpi_requests[0]);
254      #endif      #endif
255    
256      /* Next, send out the local idx */      /* Next, send out the local idx */
# Line 263  paso::SparseMatrix* Paso_MergedSolver_me Line 258  paso::SparseMatrix* Paso_MergedSolver_me
258      tag += size;      tag += size;
259      #ifdef ESYS_MPI      #ifdef ESYS_MPI
260      MPI_Issend(idx, len, MPI_INT, 0, tag, A->mpi_info->comm,      MPI_Issend(idx, len, MPI_INT, 0, tag, A->mpi_info->comm,
261              &mpi_requests[1]);                          &mpi_requests[1]);
262      #endif      #endif
263    
264      /* At last, send out the local val */      /* At last, send out the local val */
# Line 271  paso::SparseMatrix* Paso_MergedSolver_me Line 266  paso::SparseMatrix* Paso_MergedSolver_me
266      tag += size;      tag += size;
267      #ifdef ESYS_MPI      #ifdef ESYS_MPI
268      MPI_Issend(val, len, MPI_DOUBLE, 0, tag, A->mpi_info->comm,      MPI_Issend(val, len, MPI_DOUBLE, 0, tag, A->mpi_info->comm,
269              &mpi_requests[2]);                          &mpi_requests[2]);
270    
271      MPI_Waitall(3, mpi_requests, mpi_stati);      MPI_Waitall(3, mpi_requests, mpi_stati);
272      #endif      #endif
# Line 279  paso::SparseMatrix* Paso_MergedSolver_me Line 274  paso::SparseMatrix* Paso_MergedSolver_me
274      delete[] ptr;      delete[] ptr;
275      delete[] idx;      delete[] idx;
276      delete[] val;      delete[] val;
277        out.reset();
     out = NULL;  
278    }    }
279    
280    delete[] mpi_requests;    delete[] mpi_requests;

Legend:
Removed from v.4828  
changed lines
  Added in v.4829

  ViewVC Help
Powered by ViewVC 1.1.26