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

Diff of /trunk/paso/src/Smoother.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 116  void Paso_Preconditioner_Smoother_solve( Line 116  void Paso_Preconditioner_Smoother_solve(
116                      const dim_t sweeps, const bool_t x_is_initial)                      const dim_t sweeps, const bool_t x_is_initial)
117  {  {
118     const dim_t n= (A_p->mainBlock->numRows) * (A_p->mainBlock->row_block_size);     const dim_t n= (A_p->mainBlock->numRows) * (A_p->mainBlock->row_block_size);
 //   const dim_t n= (A_p->mainBlock->numRows) * (A_p->mainBlock->row_block_size) + (A_p->col_coupleBlock->numCols) * (A_p->col_coupleBlock->col_block_size);  
119        
120     double *b_new = smoother->localSmoother->buffer;     double *b_new = smoother->localSmoother->buffer;
121     dim_t nsweeps=sweeps;     dim_t nsweeps=sweeps;
# Line 234  void Paso_Preconditioner_LocalSmoother_S Line 233  void Paso_Preconditioner_LocalSmoother_S
233        x[0]*=diag[0];        x[0]*=diag[0];
234        for (i = 1; i < n; ++i) {        for (i = 1; i < n; ++i) {
235       mm=ptr_main[i];       mm=ptr_main[i];
 //   mm=A_p->pattern->ptr[i+1];  
236       /* x_i=x_i-a_ik*x_k (with k<i) - a_ik*x_k (with k>i) */       /* x_i=x_i-a_ik*x_k (with k<i) - a_ik*x_k (with k>i) */
237       rtmp=x[i];       rtmp=x[i];
238       for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<mm; ++iptr_ik) {       for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<mm; ++iptr_ik) {
239          k=A_p->pattern->index[iptr_ik];          k=A_p->pattern->index[iptr_ik];
240          if (k != i)          rtmp-=A_p->val[iptr_ik]*x[k];
           rtmp-=A_p->val[iptr_ik]*x[k];  
241       }       }
242       x[i]=rtmp*diag[i];       x[i]=rtmp*diag[i];
243        }        }
# Line 248  void Paso_Preconditioner_LocalSmoother_S Line 245  void Paso_Preconditioner_LocalSmoother_S
245        Paso_BlockOps_MViP_2(&diag[0], &x[0]);        Paso_BlockOps_MViP_2(&diag[0], &x[0]);
246        for (i = 1; i < n; ++i) {        for (i = 1; i < n; ++i) {
247       mm=ptr_main[i];       mm=ptr_main[i];
 //   mm=A_p->pattern->ptr[i+1];  
248       for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<mm; ++iptr_ik) {       for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<mm; ++iptr_ik) {
249          k=A_p->pattern->index[iptr_ik];                                    k=A_p->pattern->index[iptr_ik];                          
250          if (k != i)          Paso_BlockOps_SMV_2(&x[2*i], &A_p->val[4*iptr_ik], &x[2*k]);
           Paso_BlockOps_SMV_2(&x[2*i], &A_p->val[4*iptr_ik], &x[2*k]);  
251       }       }
252       Paso_BlockOps_MViP_2(&diag[4*i], &x[2*i]);       Paso_BlockOps_MViP_2(&diag[4*i], &x[2*i]);
253        }        }
# Line 260  void Paso_Preconditioner_LocalSmoother_S Line 255  void Paso_Preconditioner_LocalSmoother_S
255        Paso_BlockOps_MViP_3(&diag[0], &x[0]);        Paso_BlockOps_MViP_3(&diag[0], &x[0]);
256        for (i = 1; i < n; ++i) {        for (i = 1; i < n; ++i) {
257       mm=ptr_main[i];       mm=ptr_main[i];
 //   mm=A_p->pattern->ptr[i+1];  
258       for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<mm; ++iptr_ik) {       for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<mm; ++iptr_ik) {
259          k=A_p->pattern->index[iptr_ik];          k=A_p->pattern->index[iptr_ik];
260          if (k != i)          Paso_BlockOps_SMV_3(&x[3*i], &A_p->val[9*iptr_ik], &x[3*k]);
           Paso_BlockOps_SMV_3(&x[3*i], &A_p->val[9*iptr_ik], &x[3*k]);  
261       }       }
262       Paso_BlockOps_MViP_3(&diag[9*i], &x[3*i]);       Paso_BlockOps_MViP_3(&diag[9*i], &x[3*i]);
263        }        }
# Line 272  void Paso_Preconditioner_LocalSmoother_S Line 265  void Paso_Preconditioner_LocalSmoother_S
265        Paso_BlockOps_solve_N(n_block, &x[0], &diag[0], &pivot[0], &failed);        Paso_BlockOps_solve_N(n_block, &x[0], &diag[0], &pivot[0], &failed);
266        for (i = 1; i < n; ++i) {        for (i = 1; i < n; ++i) {
267       mm=ptr_main[i];       mm=ptr_main[i];
 //   mm=A_p->pattern->ptr[i+1];  
268       for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<mm; ++iptr_ik) {       for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<mm; ++iptr_ik) {
269          k=A_p->pattern->index[iptr_ik];          k=A_p->pattern->index[iptr_ik];
270          if (k != i)          Paso_BlockOps_SMV_N(n_block, &x[n_block*i], &A_p->val[block_len*iptr_ik], &x[n_block*k]);
           Paso_BlockOps_SMV_N(n_block, &x[n_block*i], &A_p->val[block_len*iptr_ik], &x[n_block*k]);  
271       }       }
272       Paso_BlockOps_solve_N(n_block, &x[n_block*i], &diag[block_len*i], &pivot[n_block*i], &failed);       Paso_BlockOps_solve_N(n_block, &x[n_block*i], &diag[block_len*i], &pivot[n_block*i], &failed);
273        }        }
274     }     }
 //fprintf(stderr, "f x[0]=%g\n", x[0]);  
275    
276     /* backward sweeps */     /* backward sweeps */
277     if (n_block==1) {     if (n_block==1) {
# Line 289  void Paso_Preconditioner_LocalSmoother_S Line 279  void Paso_Preconditioner_LocalSmoother_S
279          rtmp=x[i];          rtmp=x[i];
280          mm=ptr_main[i];          mm=ptr_main[i];
281          for (iptr_ik=mm+1; iptr_ik < A_p->pattern->ptr[i+1]; ++iptr_ik) {          for (iptr_ik=mm+1; iptr_ik < A_p->pattern->ptr[i+1]; ++iptr_ik) {
 //      for (iptr_ik=A_p->pattern->ptr[i]; iptr_ik < A_p->pattern->ptr[i+1]; ++iptr_ik) {  
282             k=A_p->pattern->index[iptr_ik];             k=A_p->pattern->index[iptr_ik];
283             if (k != i)             rtmp-=A_p->val[iptr_ik]*x[k];
              rtmp-=A_p->val[iptr_ik]*x[k];  
284          }          }
285          x[i]=diag[i]*rtmp;          x[i]=diag[i]*rtmp;
286        }        }
# Line 301  void Paso_Preconditioner_LocalSmoother_S Line 289  void Paso_Preconditioner_LocalSmoother_S
289        for (i = n-2; i > -1; --i) {        for (i = n-2; i > -1; --i) {
290          mm=ptr_main[i];          mm=ptr_main[i];
291              for (iptr_ik=mm+1; iptr_ik < A_p->pattern->ptr[i+1]; ++iptr_ik) {              for (iptr_ik=mm+1; iptr_ik < A_p->pattern->ptr[i+1]; ++iptr_ik) {
 //      for (iptr_ik=A_p->pattern->ptr[i]; iptr_ik < A_p->pattern->ptr[i+1]; ++iptr_ik) {  
   
292             k=A_p->pattern->index[iptr_ik];             k=A_p->pattern->index[iptr_ik];
293             if (k != i)             Paso_BlockOps_SMV_2(&x[2*i], &A_p->val[4*iptr_ik], &x[2*k]);
              Paso_BlockOps_SMV_2(&x[2*i], &A_p->val[4*iptr_ik], &x[2*k]);  
294          }          }
295          Paso_BlockOps_MViP_2(&diag[i*4], &x[2*i]);          Paso_BlockOps_MViP_2(&diag[i*4], &x[2*i]);
296                    
# Line 314  void Paso_Preconditioner_LocalSmoother_S Line 299  void Paso_Preconditioner_LocalSmoother_S
299        for (i = n-2; i > -1; --i) {        for (i = n-2; i > -1; --i) {
300          mm=ptr_main[i];          mm=ptr_main[i];
301              for (iptr_ik=mm+1; iptr_ik < A_p->pattern->ptr[i+1]; ++iptr_ik) {              for (iptr_ik=mm+1; iptr_ik < A_p->pattern->ptr[i+1]; ++iptr_ik) {
 //      for (iptr_ik=A_p->pattern->ptr[i]; iptr_ik < A_p->pattern->ptr[i+1]; ++iptr_ik) {  
302             k=A_p->pattern->index[iptr_ik];                 k=A_p->pattern->index[iptr_ik];    
303             if (k != i)             Paso_BlockOps_SMV_3(&x[3*i], &A_p->val[9*iptr_ik], &x[3*k]);
              Paso_BlockOps_SMV_3(&x[3*i], &A_p->val[9*iptr_ik], &x[3*k]);  
304          }          }
305          Paso_BlockOps_MViP_3(&diag[i*9], &x[3*i]);          Paso_BlockOps_MViP_3(&diag[i*9], &x[3*i]);
306        }        }
# Line 326  void Paso_Preconditioner_LocalSmoother_S Line 309  void Paso_Preconditioner_LocalSmoother_S
309        for (i = n-2; i > -1; --i) {        for (i = n-2; i > -1; --i) {
310       mm=ptr_main[i];       mm=ptr_main[i];
311           for (iptr_ik=mm+1; iptr_ik < A_p->pattern->ptr[i+1]; ++iptr_ik) {           for (iptr_ik=mm+1; iptr_ik < A_p->pattern->ptr[i+1]; ++iptr_ik) {
 //   for (iptr_ik=A_p->pattern->ptr[i]; iptr_ik < A_p->pattern->ptr[i+1]; ++iptr_ik) {  
312          k=A_p->pattern->index[iptr_ik];              k=A_p->pattern->index[iptr_ik];    
313          if (k != i)          Paso_BlockOps_SMV_N(n_block, &y[0], &A_p->val[block_len*iptr_ik], &x[n_block*k]);
           Paso_BlockOps_SMV_N(n_block, &y[0], &A_p->val[block_len*iptr_ik], &x[n_block*k]);  
314       }       }
315       Paso_BlockOps_Cpy_N(n_block ,&x[n_block*i], &y[0]);       Paso_BlockOps_Cpy_N(n_block ,&x[n_block*i], &y[0]);
316       Paso_BlockOps_solve_N(n_block, &x[n_block*i], &diag[i*block_len], &pivot[i*n_block], &failed);       Paso_BlockOps_solve_N(n_block, &x[n_block*i], &diag[i*block_len], &pivot[i*n_block], &failed);
317        }          }  
318        TMPMEMFREE(y);        TMPMEMFREE(y);
 //fprintf(stderr, "b x[0]=%g\n", x[0]);  
319     }     }
320        
321     if (failed > 0) {     if (failed > 0) {

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

  ViewVC Help
Powered by ViewVC 1.1.26