/[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 3379 by gross, Wed Nov 24 04:48:49 2010 UTC revision 3402 by gross, Tue Dec 7 07:36:12 2010 UTC
# Line 166  void Paso_Preconditioner_LocalSmoother_S Line 166  void Paso_Preconditioner_LocalSmoother_S
166  {  {
167     const dim_t nt=omp_get_max_threads();     const dim_t nt=omp_get_max_threads();
168     if (smoother->Jacobi) {     if (smoother->Jacobi) {
169        Paso_BlockOps_allMV(A->row_block_size,A->numRows,smoother->diag,smoother->pivot,x);        Paso_BlockOps_solveAll(A->row_block_size,A->numRows,smoother->diag,smoother->pivot,x);
170     } else {     } else {
171        if (nt < 2) {        if (nt < 2) {
172       Paso_Preconditioner_LocalSmoother_Sweep_sequential(A,smoother,x);       Paso_Preconditioner_LocalSmoother_Sweep_sequential(A,smoother,x);
# Line 183  void Paso_Preconditioner_LocalSmoother_S Line 183  void Paso_Preconditioner_LocalSmoother_S
183     const dim_t n=A_p->numRows;     const dim_t n=A_p->numRows;
184     const dim_t n_block=A_p->row_block_size;     const dim_t n_block=A_p->row_block_size;
185     const double *diag = smoother->diag;     const double *diag = smoother->diag;
186     /* const index_t* pivot = smoother->pivot;     const index_t* pivot = smoother->pivot;
187     const dim_t block_size=A_p->block_size;  use for block size >3*/     const dim_t block_len=A_p->block_size;
188        
189     register dim_t i,k;     register dim_t i,k;
190     register index_t iptr_ik, mm;     register index_t iptr_ik, mm;
191       register double rtmp;
192       int failed = 0;
193        
194     const index_t* ptr_main = Paso_SparseMatrix_borrowMainDiagonalPointer(A_p);     const index_t* ptr_main = Paso_SparseMatrix_borrowMainDiagonalPointer(A_p);
195     /* forward substitution */     /* forward substitution */
# Line 197  void Paso_Preconditioner_LocalSmoother_S Line 199  void Paso_Preconditioner_LocalSmoother_S
199        for (i = 1; i < n; ++i) {        for (i = 1; i < n; ++i) {
200       mm=ptr_main[i];       mm=ptr_main[i];
201       /* x_i=x_i-a_ik*x_k  (with k<i) */       /* x_i=x_i-a_ik*x_k  (with k<i) */
202         rtmp=x[i];
203       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) {
204          k=A_p->pattern->index[iptr_ik];          k=A_p->pattern->index[iptr_ik];
205          x[i]-=A_p->val[iptr_ik]*x[k];          rtmp-=A_p->val[iptr_ik]*x[k];
206       }       }
207       x[i]*=diag[i];       x[i]=rtmp*diag[i];
208        }        }
209     } else if (n_block==2) {     } else if (n_block==2) {
210        Paso_BlockOps_MV_2(&x[0], &diag[0], &x[0]);        Paso_BlockOps_MViP_2(&diag[0], &x[0]);
211        for (i = 1; i < n; ++i) {        for (i = 1; i < n; ++i) {
212       mm=ptr_main[i];       mm=ptr_main[i];
213       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) {
214          k=A_p->pattern->index[iptr_ik];                                    k=A_p->pattern->index[iptr_ik];                          
215          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]);
216       }       }
217       Paso_BlockOps_MV_2(&x[2*i], &diag[4*i], &x[2*i]);       Paso_BlockOps_MViP_2(&diag[4*i], &x[2*i]);
218        }        }
219     } else if (n_block==3) {     } else if (n_block==3) {
220        Paso_BlockOps_MV_3(&x[0], &diag[0], &x[0]);        Paso_BlockOps_MViP_3(&diag[0], &x[0]);
221        for (i = 1; i < n; ++i) {        for (i = 1; i < n; ++i) {
222       mm=ptr_main[i];       mm=ptr_main[i];
223       /* x_i=x_i-a_ik*x_k */       /* x_i=x_i-a_ik*x_k */
# Line 222  void Paso_Preconditioner_LocalSmoother_S Line 225  void Paso_Preconditioner_LocalSmoother_S
225          k=A_p->pattern->index[iptr_ik];          k=A_p->pattern->index[iptr_ik];
226          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]);
227       }       }
228       Paso_BlockOps_MV_3(&x[3*i], &diag[9*i], &x[3*i]);       Paso_BlockOps_MViP_3(&diag[9*i], &x[3*i]);
229        }        }
230     } /* add block size >3 */     } else {
231          Paso_BlockOps_solve_N(n_block, &x[0], &diag[0], &pivot[0], &failed);
232          for (i = 1; i < n; ++i) {
233         mm=ptr_main[i];
234         /* x_i=x_i-a_ik*x_k */
235         for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<mm; ++iptr_ik) {
236            k=A_p->pattern->index[iptr_ik];
237            Paso_BlockOps_SMV_N(n_block, &x[n_block*i], &A_p->val[block_len*iptr_ik], &x[n_block*k]);
238         }
239         Paso_BlockOps_solve_N(n_block, &x[n_block*i], &diag[block_len*i], &pivot[n_block*i], &failed)
240          }
241       }
242     /* backward substitution */     /* backward substitution */
243        
244     if (n_block==1) {     if (n_block==1) {
245        for (i = n-2; i > -1; --i) {                for (i = n-2; i > -1; --i) {        
246          mm=ptr_main[i];          mm=ptr_main[i];
247          Paso_BlockOps_MV_1(&x[i], &A_p->val[mm], &x[i]);          rtmp=x[i]*A_p->val[mm];
248          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) {
249             k=A_p->pattern->index[iptr_ik];               k=A_p->pattern->index[iptr_ik];
250             Paso_BlockOps_SMV_1(&x[i], &A_p->val[iptr_ik], &x[k]);             rtmp-=A_p->val[iptr_ik]*x[k];
251          }          }
252          Paso_BlockOps_MV_1(&x[i], &diag[i], &x[i]);          x[i]=diag[i]*rtmp;
253        }        }
254                
255     } else if (n_block==2) {     } else if (n_block==2) {
256        for (i = n-2; i > -1; --i) {        for (i = n-2; i > -1; --i) {
257          mm=ptr_main[i];          mm=ptr_main[i];
258          Paso_BlockOps_MV_2(&x[2*i], &A_p->val[4*mm], &x[2*i]);          Paso_BlockOps_MViP_2(&A_p->val[4*mm], &x[2*i]);
259          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) {
260             k=A_p->pattern->index[iptr_ik];             k=A_p->pattern->index[iptr_ik];
261             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]);
262          }          }
263          Paso_BlockOps_MV_2(&x[2*i], &diag[i*4], &x[2*i]);          Paso_BlockOps_MViP_2(&diag[i*4], &x[2*i]);
264                    
265        }        }
266     } else if (n_block==3) {     } else if (n_block==3) {
267        for (i = n-2; i > -1; --i) {        for (i = n-2; i > -1; --i) {
268          mm=ptr_main[i];          mm=ptr_main[i];
269          Paso_BlockOps_MV_3(&x[3*i], &A_p->val[9*mm], &x[3*i]);          Paso_BlockOps_MViP_3(&A_p->val[9*mm], &x[3*i]);
       
270          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) {
271             k=A_p->pattern->index[iptr_ik];                 k=A_p->pattern->index[iptr_ik];    
272             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]);
273          }          }
274          Paso_BlockOps_MV_3(&x[3*i], &diag[i*9], &x[3*i]);          Paso_BlockOps_MViP_3(&diag[i*9], &x[3*i]);
275        }        }
276                
277     } /* add block size >3 */           } else {
278          double *y=TMPMEMALLOC(n_block, double);
279          for (i = n-2; i > -1; --i) {
280         mm=ptr_main[i];
281         Paso_BlockOps_MV_N(n_block, &y[0], &A_p->val[block_len*mm], &x[n_block*i]);
282         for (iptr_ik=mm+1; iptr_ik < A_p->pattern->ptr[i+1]; ++iptr_ik) {
283            k=A_p->pattern->index[iptr_ik];    
284            Paso_BlockOps_SMV_N(n_block, &y[0], &A_p->val[block_len*iptr_ik], &x[n_block*k]);
285         }
286         Paso_BlockOps_Cpy_N(n_block ,&x[n_block*i], &y[0]);
287         Paso_BlockOps_solve_N(n_block, &x[n_block*i], &diag[i*block_len], &pivot[i*n_block], &failed);
288          }  
289          TMPMEMFREE(y);
290       }
291      
292       if (failed > 0) {
293          Esys_setError(ZERO_DIVISION_ERROR, "Paso_Preconditioner_LocalSmoother_Sweep_sequential: non-regular main diagonal block.");
294       }
295        
296     return;     return;
297  }  }
# Line 273  void Paso_Preconditioner_LocalSmoother_S Line 302  void Paso_Preconditioner_LocalSmoother_S
302     const dim_t n_block=A_p->row_block_size;     const dim_t n_block=A_p->row_block_size;
303     const double *diag = smoother->diag;     const double *diag = smoother->diag;
304     index_t* pivot = smoother->pivot;     index_t* pivot = smoother->pivot;
305     const dim_t block_size=A_p->block_size;     const dim_t block_len=A_p->block_size;
306       double *y;
307        
308     register dim_t i,k;     register dim_t i,k;
309     register index_t color,iptr_ik, mm;     register index_t color,iptr_ik, mm;
310       register double rtmp;
311       int failed = 0;
312        
313     const index_t* coloring = Paso_Pattern_borrowColoringPointer(A_p->pattern);     const index_t* coloring = Paso_Pattern_borrowColoringPointer(A_p->pattern);
314     const dim_t num_colors = Paso_Pattern_getNumColors(A_p->pattern);     const dim_t num_colors = Paso_Pattern_getNumColors(A_p->pattern);
315     const index_t* ptr_main = Paso_SparseMatrix_borrowMainDiagonalPointer(A_p);     const index_t* ptr_main = Paso_SparseMatrix_borrowMainDiagonalPointer(A_p);
     
    /* forward substitution */  
316    
317         #pragma omp parallel  private(mm, i,iptr_ik,k,rtmp, color, y)
318     /* color = 0 */     {
319     if (n_block==1) {        if (n_block>3) {
320        #pragma omp parallel for schedule(static) private(i)       y=TMPMEMALLOC(n_block, double);
321        for (i = 0; i < n; ++i) {        } else {
322       if (coloring[i]== 0 ) Paso_BlockOps_MV_1(&x[i], &diag[i], &x[i]);       y=NULL;
323        }        }
324     } else if (n_block==2) {        /* forward substitution */
325           #pragma omp parallel for schedule(static) private(i)  
326          /* color = 0 */
327          if (n_block==1) {
328         #pragma omp  for schedule(static)
329       for (i = 0; i < n; ++i) {       for (i = 0; i < n; ++i) {
330          if (coloring[i]== 0 ) Paso_BlockOps_MV_2(&x[2*i], &diag[i*4], &x[2*i]);          if (coloring[i]== 0 ) x[i]*=diag[i];
331       }       }
332      } else if (n_block==3) {        } else if (n_block==2) {
333       #pragma omp parallel for schedule(static) private(i)           #pragma omp for schedule(static)
334       for (i = 0; i < n; ++i) {       for (i = 0; i < n; ++i) {
335          if (coloring[i]==0) Paso_BlockOps_MV_3(&x[3*i], &diag[i*9], &x[3*i]);          if (coloring[i]== 0 ) Paso_BlockOps_MViP_2(&diag[i*4], &x[2*i]);
336         }
337          } else if (n_block==3) {
338         #pragma omp for schedule(static)
339         for (i = 0; i < n; ++i) {
340            if (coloring[i]==0) Paso_BlockOps_MViP_3(&diag[i*9], &x[3*i]);
341         }
342          } else {
343         #pragma omp for schedule(static)
344         for (i = 0; i < n; ++i) {
345            if (coloring[i]==0) Paso_BlockOps_solve_N(n_block, &x[n_block*i], &diag[block_len*i], &pivot[n_block*i], &failed);
346       }       }
    } else {  
       #pragma omp parallel for schedule(static) private(i)  
       for (i = 0; i < n; ++i) {  
      if (coloring[i]==0) Paso_BlockOps_Solve_N(n_block, &x[n_block*i], &diag[block_size*i], &pivot[n_block*i], &x[n_block*i]);  
347        }        }
    }  
348        
349     for (color=1;color<num_colors;++color) {  
350          for (color=1;color<num_colors;++color) {
351    
352        if (n_block==1) {       if (n_block==1) {
353       #pragma omp parallel for schedule(static) private(i,iptr_ik,k)          #pragma omp for schedule(static)
354       for (i = 0; i < n; ++i) {          for (i = 0; i < n; ++i) {
355          if (coloring[i]==color) {             if (coloring[i]==color) {
356             /* x_i=x_i-a_ik*x_k */                                /* x_i=x_i-a_ik*x_k */  
357             for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<A_p->pattern->ptr[i+1]; ++iptr_ik) {            rtmp=x[i];
358            k=A_p->pattern->index[iptr_ik];            for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<A_p->pattern->ptr[i+1]; ++iptr_ik) {
359            if (coloring[k]<color) Paso_BlockOps_SMV_1(&x[i], &A_p->val[iptr_ik], &x[k]);               k=A_p->pattern->index[iptr_ik];
360                 if (coloring[k]<color) rtmp-=A_p->val[iptr_ik]*x[k];
361              }
362              x[i]=diag[i]*rtmp;
363             }             }
            Paso_BlockOps_MV_1(&x[i], &diag[i], &x[i]);  
364          }          }
365       }       } else if (n_block==2) {
366        } else if (n_block==2) {          #pragma omp for schedule(static)
367       #pragma omp parallel for schedule(static) private(i,iptr_ik,k)          for (i = 0; i < n; ++i) {
368       for (i = 0; i < n; ++i) {             if (coloring[i]==color) {
369          if (coloring[i]==color) {            for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<A_p->pattern->ptr[i+1]; ++iptr_ik) {
370             for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<A_p->pattern->ptr[i+1]; ++iptr_ik) {               k=A_p->pattern->index[iptr_ik];                          
371            k=A_p->pattern->index[iptr_ik];                                         if (coloring[k]<color) Paso_BlockOps_SMV_2(&x[2*i], &A_p->val[4*iptr_ik], &x[2*k]);
372            if (coloring[k]<color) Paso_BlockOps_SMV_2(&x[2*i], &A_p->val[4*iptr_ik], &x[2*k]);            }
373              Paso_BlockOps_MViP_2(&diag[4*i], &x[2*i]);
374             }             }
            Paso_BlockOps_MV_2(&x[2*i], &diag[4*i], &x[2*i]);  
375          }          }
376       }       } else if (n_block==3) {
377        } else if (n_block==3) {          #pragma omp for schedule(static)
378       #pragma omp parallel for schedule(static) private(i,iptr_ik,k)          for (i = 0; i < n; ++i) {
379       for (i = 0; i < n; ++i) {             if (coloring[i]==color) {
380          if (coloring[i]==color) {            for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<A_p->pattern->ptr[i+1]; ++iptr_ik) {
381             for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<A_p->pattern->ptr[i+1]; ++iptr_ik) {               k=A_p->pattern->index[iptr_ik];                          
382            k=A_p->pattern->index[iptr_ik];                                         if (coloring[k]<color) Paso_BlockOps_SMV_3(&x[3*i], &A_p->val[9*iptr_ik], &x[3*k]);
383            if (coloring[k]<color) Paso_BlockOps_SMV_3(&x[3*i], &A_p->val[9*iptr_ik], &x[3*k]);            }
384              Paso_BlockOps_MViP_3(&diag[9*i], &x[3*i]);
385             }             }
            Paso_BlockOps_MV_3(&x[3*i], &diag[9*i], &x[3*i]);  
386          }          }
387       }       } else {
388        } else {          #pragma omp for schedule(static)
389       #pragma omp parallel for schedule(static) private(i,iptr_ik,k)          for (i = 0; i < n; ++i) {
390       for (i = 0; i < n; ++i) {             if (coloring[i] == color) {
391          if (coloring[i] == color) {            for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<A_p->pattern->ptr[i+1]; ++iptr_ik) {
392             for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<A_p->pattern->ptr[i+1]; ++iptr_ik) {               k=A_p->pattern->index[iptr_ik];                          
393            k=A_p->pattern->index[iptr_ik];                                         if (coloring[k]<color) Paso_BlockOps_SMV_N(n_block, &x[n_block*i], &A_p->val[block_len*iptr_ik], &x[n_block*k]);
394            if (coloring[k]<color) Paso_BlockOps_SMV_N(n_block, &x[n_block*i], &A_p->val[block_size*iptr_ik], &x[n_block*k]);            }
395              Paso_BlockOps_solve_N(n_block, &x[n_block*i], &diag[block_len*i], &pivot[n_block*i], &failed);
396             }             }
            Paso_BlockOps_Solve_N(n_block, &x[n_block*i], &diag[block_size*i], &pivot[n_block*i], &x[n_block*i]);  
397          }          }
398       }       }
399        }        } /* end of coloring loop */
400     } /* end of coloring loop */        
401            /* backward substitution */
402     /* backward substitution */        for (color=(num_colors)-2 ;color>-1;--color) { /* Note: color=(num_colors)-1 is not required */
403     for (color=(num_colors)-2 ;color>-1;--color) { /* Note: color=(num_colors)-1 is not required */       if (n_block==1) {
404        if (n_block==1) {          #pragma omp for schedule(static)
405       #pragma omp parallel for schedule(static) private(mm, i,iptr_ik,k)          for (i = 0; i < n; ++i) {
406       for (i = 0; i < n; ++i) {             if (coloring[i]==color) {
407          if (coloring[i]==color) {            mm=ptr_main[i];
408             mm=ptr_main[i];            rtmp=A_p->val[mm]*x[i];
409             Paso_BlockOps_MV_1(&x[i], &A_p->val[mm], &x[i]);            for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<A_p->pattern->ptr[i+1]; ++iptr_ik) {
410             for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<A_p->pattern->ptr[i+1]; ++iptr_ik) {               k=A_p->pattern->index[iptr_ik];                          
411            k=A_p->pattern->index[iptr_ik];                                         if (coloring[k]>color) rtmp-=A_p->val[iptr_ik]*x[k];
412            if (coloring[k]>color) Paso_BlockOps_SMV_1(&x[i], &A_p->val[iptr_ik], &x[k]);            }
413              x[i]= rtmp*diag[i];
414             }             }
            Paso_BlockOps_MV_1(&x[i], &diag[i], &x[i]);  
415          }          }
416       }       } else if (n_block==2) {
417        } else if (n_block==2) {          #pragma omp for schedule(static)
418       #pragma omp parallel for schedule(static) private(mm, i,iptr_ik,k)          for (i = 0; i < n; ++i) {
419       for (i = 0; i < n; ++i) {             if (coloring[i]==color) {
420          if (coloring[i]==color) {            mm=ptr_main[i];
421             mm=ptr_main[i];            Paso_BlockOps_MViP_2(&A_p->val[4*mm], &x[2*i]);
422             Paso_BlockOps_MV_2(&x[2*i], &A_p->val[4*mm], &x[2*i]);            for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<A_p->pattern->ptr[i+1]; ++iptr_ik) {
423             for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<A_p->pattern->ptr[i+1]; ++iptr_ik) {               k=A_p->pattern->index[iptr_ik];                          
424            k=A_p->pattern->index[iptr_ik];                                         if (coloring[k]>color) Paso_BlockOps_SMV_2(&x[2*i], &A_p->val[4*iptr_ik], &x[2*k]);
425            if (coloring[k]>color) Paso_BlockOps_SMV_2(&x[2*i], &A_p->val[4*iptr_ik], &x[2*k]);            }
426              Paso_BlockOps_MViP_2(&diag[4*i], &x[2*i]);
427             }             }
            Paso_BlockOps_MV_2(&x[2*i], &diag[4*i], &x[2*i]);  
428          }          }
429       }       } else if (n_block==3) {
430        } else if (n_block==3) {          #pragma omp for schedule(static)
431       #pragma omp parallel for schedule(static) private(mm, i,iptr_ik,k)          for (i = 0; i < n; ++i) {
432       for (i = 0; i < n; ++i) {             if (coloring[i]==color) {
433          if (coloring[i]==color) {            mm=ptr_main[i];
434             mm=ptr_main[i];            Paso_BlockOps_MViP_3(&A_p->val[9*mm], &x[3*i]);
435             Paso_BlockOps_MV_3(&x[3*i], &A_p->val[9*mm], &x[3*i]);            for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<A_p->pattern->ptr[i+1]; ++iptr_ik) {
436             for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<A_p->pattern->ptr[i+1]; ++iptr_ik) {               k=A_p->pattern->index[iptr_ik];                          
437            k=A_p->pattern->index[iptr_ik];                                         if (coloring[k]>color) Paso_BlockOps_SMV_3(&x[3*i], &A_p->val[9*iptr_ik], &x[3*k]);
438            if (coloring[k]>color) Paso_BlockOps_SMV_3(&x[3*i], &A_p->val[9*iptr_ik], &x[3*k]);            }
439              Paso_BlockOps_MViP_3(&diag[9*i], &x[3*i]);
440             }             }
            Paso_BlockOps_MV_3(&x[3*i], &diag[9*i], &x[3*i]);  
441          }          }
442       }       } else {
443        } else {          #pragma omp for schedule(static)
444       #pragma omp parallel for schedule(static) private(mm, i,iptr_ik,k)          for (i = 0; i < n; ++i) {
445       for (i = 0; i < n; ++i) {             if (coloring[i]==color) {
446          if (coloring[i]==color) {            mm=ptr_main[i];
447             mm=ptr_main[i];            Paso_BlockOps_MV_N(n_block, &y[0], &A_p->val[block_len*mm], &x[n_block*i]);
448             Paso_BlockOps_MV_N( n_block, &x[n_block*i], &A_p->val[block_size*mm], &x[n_block*i] );            for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<A_p->pattern->ptr[i+1]; ++iptr_ik) {
449             for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<A_p->pattern->ptr[i+1]; ++iptr_ik) {               k=A_p->pattern->index[iptr_ik];                          
450            k=A_p->pattern->index[iptr_ik];                                         if (coloring[k]>color) Paso_BlockOps_SMV_N(n_block, &y[0], &A_p->val[block_len*iptr_ik], &x[n_block*k]);
451            if (coloring[k]>color) Paso_BlockOps_SMV_N(n_block, &x[n_block*i], &A_p->val[block_size*iptr_ik], &x[n_block*k]);            }
452              Paso_BlockOps_Cpy_N(n_block ,&x[n_block*i], &y[0]);
453              Paso_BlockOps_solve_N(n_block, &x[n_block*i], &diag[i*block_len], &pivot[i*n_block], &failed);
454             }             }
            Paso_BlockOps_Solve_N(n_block, &x[n_block*i], &diag[block_size*i], &pivot[n_block*i], &x[n_block*i]);  
455          }          }
456       }       }
457        }        }
458          TMPMEMFREE(y);
459       }
460       if (failed > 0) {
461          Esys_setError(ZERO_DIVISION_ERROR, "Paso_Preconditioner_LocalSmoother_Sweep_colored: non-regular main diagonal block.");
462     }     }
463     return;     return;
464  }  }

Legend:
Removed from v.3379  
changed lines
  Added in v.3402

  ViewVC Help
Powered by ViewVC 1.1.26