/[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 3096 by gross, Fri Aug 13 08:38:06 2010 UTC revision 3097 by gross, Fri Aug 20 04:59:12 2010 UTC
# Line 14  Line 14 
14    
15  /**************************************************************/  /**************************************************************/
16    
17  /* Paso: GS preconditioner with reordering                 */  /* Paso: Gauss-Seidel                */
18    
19  /**************************************************************/  /**************************************************************/
20    
 /* Copyrights by ACcESS Australia 2003,2004,2005,2006,2007,2008  */  
21  /* Author: artak@uq.edu.au                                   */  /* Author: artak@uq.edu.au                                   */
22    
23  /**************************************************************/  /**************************************************************/
# Line 26  Line 25 
25  #include "Paso.h"  #include "Paso.h"
26  #include "Solver.h"  #include "Solver.h"
27  #include "PasoUtil.h"  #include "PasoUtil.h"
28    #include "BlockOps.h"
29    
30  #include <stdio.h>  #include <stdio.h>
31    
# Line 44  void Paso_Solver_LocalGS_free(Paso_Solve Line 44  void Paso_Solver_LocalGS_free(Paso_Solve
44     if (in!=NULL) {     if (in!=NULL) {
45        MEMFREE(in->diag);        MEMFREE(in->diag);
46        MEMFREE(in->pivot);        MEMFREE(in->pivot);
47          MEMFREE(in->buffer);
48        MEMFREE(in);        MEMFREE(in);
49     }     }
50  }  }
# Line 81  Paso_Solver_LocalGS* Paso_Solver_getLoca Line 82  Paso_Solver_LocalGS* Paso_Solver_getLoca
82                
83        out->diag=MEMALLOC( ((size_t) n) * ((size_t) block_size),double);        out->diag=MEMALLOC( ((size_t) n) * ((size_t) block_size),double);
84        out->pivot=MEMALLOC( ((size_t) n) * ((size_t)  n_block), index_t);        out->pivot=MEMALLOC( ((size_t) n) * ((size_t)  n_block), index_t);
85          out->buffer=MEMALLOC( ((size_t) n) * ((size_t)  n_block), double);
86        out->sweeps=sweeps;        out->sweeps=sweeps;
87                
88        if ( ! ( Paso_checkPtr(out->diag) || Paso_checkPtr(out->pivot) ) ) {        if ( ! ( Paso_checkPtr(out->diag) || Paso_checkPtr(out->pivot) ) ) {
# Line 99  Paso_Solver_LocalGS* Paso_Solver_getLoca Line 101  Paso_Solver_LocalGS* Paso_Solver_getLoca
101     }     }
102  }  }
103    
104  /**************************************************************/  /*
   
 /* Apply Gauss-Seidel                                  
   
    in fact it solves Ax=b in two steps:  
     
    step1: among different MPI ranks, we use block Jacobi  
   
     x{k} = x{k-1} + D{-1}(b-A*x{k-1})  
   
    => D*x{k} = b - (E+F)x{k-1}  
   
 where matrix D is (let p be the number of nodes):  
   
 --------------------  
 |A1|  |  | ...  |  |  
 --------------------  
 |  |A2|  | ...  |  |  
 --------------------  
 |  |  |A3| ...  |  |  
 --------------------  
 |          ...     |  
 --------------------  
 |  |  |  | ...  |Ap|  
 --------------------  
   
 and Ai (i \in [1,p]) represents the mainBlock of matrix  
 A on rank i. Matrix (E+F) is represented as the coupleBlock  
 of matrix A on each rank (annotated as ACi).  
   
 Therefore, step1 can be turned into the following for rank i:  
   
 => Ai * x{k} = b - ACi * x{k-1}  
105    
106  where both x{k} and b are the segment of x and b on node i,  performs a few sweeps of the  from
 and x{k-1} is the old segment values of x on all other nodes.  
107    
108  step2: inside rank i, we use Gauss-Seidel  S (x_{k} -  x_{k-1}) = b - A x_{k-1}
109    
110  let b'= b - ACi * x{k-1} we have Ai * x{k} = b' for rank i  where x_{0}=0 and S provides some approximatioon of A.
 by using symetrix Gauss-Seidel,  
111    
112  this can be solved in a forward phase and a backward phase:  Under MPI S is build on using A_p->mainBlock only.
113    if GS is local the defect b - A x_{k-1} is calculated using A_p->mainBlock only.
114    
    forward phase:  x{m} = diag(Ai){-1} (b' - E*x{m} - F*x{m-1})  
    backward phase: x{m+1} = diag(Ai){-1} (b' - F*{m+1} - E*x{m})  
     
115  */  */
116    
117  void Paso_Solver_solveGS(Paso_SystemMatrix* A_p, Paso_Solver_GS * gs, double * x, const double * b)  void Paso_Solver_solveGS(Paso_SystemMatrix* A_p, Paso_Solver_GS * gs, double * x, const double * b)
118  {  {
119     register dim_t i;     register dim_t i;
120     const dim_t n=A_p->mainBlock->numRows;     const dim_t n= (A_p->mainBlock->numRows) * (A_p->mainBlock->row_block_size);
121     const dim_t n_block=A_p->mainBlock->row_block_size;    
122     double *remote_x=NULL;     double *b_new = gs->localGS->buffer;
    const dim_t sweeps_ref=gs->localGS->sweeps;  
123     dim_t sweeps=gs->localGS->sweeps;     dim_t sweeps=gs->localGS->sweeps;
    const bool_t remote = (!gs->is_local) && (A_p->mpi_info->size > 1);  
    double *new_b = NULL;  
124        
125     if (remote) {     if (gs->is_local) {
126        new_b=TMPMEMALLOC(n*n_block,double);        Paso_Solver_solveLocalGS(A_p->mainBlock,gs->localGS,x,b);
127        gs->localGS->sweeps=(dim_t) ceil( sqrt(DBLE(gs->localGS->sweeps)) );     } else {
128     }        #pragma omp parallel for private(i) schedule(static)
129          for (i=0;i<n;++i) x[i]=b[i];
130          
131          Paso_Solver_localGSSweep(A_p->mainBlock,gs->localGS,x);
132          while (sweeps > 0 ) {
133         #pragma omp parallel for private(i) schedule(static)
134         for (i=0;i<n;++i) b_new[i]=b[i];
135    
136             Paso_SystemMatrix_MatrixVector((-1.), A_p, x, 1., b_new); /* b_new = b - A*x */
137        
138         Paso_Solver_localGSSweep(A_p->mainBlock,gs->localGS,b_new);
139        
140         #pragma omp parallel for private(i) schedule(static)
141         for (i=0;i<n;++i) x[i]+=b_new[i];
142        
143         sweeps--;
144          }
145          
146       }
147    }
148    void Paso_Solver_solveLocalGS(Paso_SparseMatrix* A_p, Paso_Solver_LocalGS * gs, double * x, const double * b)
149    {
150       register dim_t i;
151       const dim_t n= (A_p->numRows) * (A_p->row_block_size);
152       double *b_new = gs->buffer;
153       dim_t sweeps=gs->sweeps;
154      
155       #pragma omp parallel for private(i) schedule(static)
156       for (i=0;i<n;++i) x[i]=b[i];
157       Paso_Solver_localGSSweep(A_p,gs,x);
158        
159     Paso_Solver_solveLocalGS(A_p->mainBlock,gs->localGS,x,b);     while (sweeps > 0 ) {
    sweeps-=gs->localGS->sweeps;  
     
    while (sweeps > 0 && remote ) {  
          Paso_SystemMatrix_startCollect(A_p,x);  
      /* calculate new_b = b - ACi * x{k-1}, where x{k-1} are remote  
         value of x, which requires MPI communications */  
160       #pragma omp parallel for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
161       for (i=0;i<n*n_block;++i) new_b[i]=b[i];       for (i=0;i<n;++i) b_new[i]=b[i];
162                
163       remote_x=Paso_SystemMatrix_finishCollect(A_p);       Paso_SparseMatrix_MatrixVector_CSC_OFFSET0((-1.), A_p, x, 1., b_new); /* b_new = b - A*x */
164       /*new_b = (-1) * AC * x + 1. * new_b */      
165       Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1),A_p->col_coupleBlock,remote_x,DBLE(1), new_b);       Paso_Solver_localGSSweep(A_p,gs,b_new);
166        
167         #pragma omp parallel for private(i) schedule(static)
168         for (i=0;i<n;++i) x[i]+=b_new[i];
169            
170       Paso_Solver_solveLocalGS(A_p->mainBlock,gs->localGS,x,new_b);       sweeps--;
      sweeps-=gs->localGS->sweeps;  
171     }     }
    TMPMEMFREE(new_b);  
    gs->localGS->sweeps=sweeps_ref;  
    return;  
172  }  }
173    
174  void Paso_Solver_solveLocalGS(Paso_SparseMatrix* A, Paso_Solver_LocalGS * gs, double * x, const double * b)  void Paso_Solver_localGSSweep(Paso_SparseMatrix* A, Paso_Solver_LocalGS * gs, double * x)
175  {  {
    dim_t i;  
176     #ifdef _OPENMP     #ifdef _OPENMP
177     const dim_t nt=omp_get_max_threads();     const dim_t nt=omp_get_max_threads();
178     #else     #else
179     const dim_t nt = 1;     const dim_t nt = 1;
180     #endif     #endif
181     gs->sweeps=MAX(gs->sweeps,1);     if (nt > 1) {
182            Paso_Solver_localGSSweep_sequential(A,gs,x);
183     for (i =0 ; i<gs->sweeps; i++) {     } else {
184        if (nt > 1) {        Paso_Solver_localGSSweep_colored(A,gs,x);
      Paso_Solver_solveLocalGS_sequential(A,gs,x,b);  
       } else {  
      Paso_Solver_solveLocalGS_colored(A,gs,x,b);  
         /* Paso_Solver_solveLocalGS_tiled(A,gs,x,b); LIN: ADD YOUR STUFF */  
       }  
185     }     }
186  }  }
187    
188  /*  /* inplace Gauss-Seidel sweep in seqential mode: */
   
    applies symmetric Gauss Seidel with coloring = (U+D)^{-1}*D* (L+D)^{-1}  
189    
190  */  void Paso_Solver_localGSSweep_sequential(Paso_SparseMatrix* A_p, Paso_Solver_LocalGS * gs, double * x)
191    {
192       const dim_t n=A_p->numRows;
193       const dim_t n_block=A_p->row_block_size;
194       const double *diag = gs->diag;
195       /* const index_t* pivot = gs->pivot;
196       const dim_t block_size=A_p->block_size;  use for block size >3*/
197      
198       register dim_t i,k;
199       register index_t iptr_ik, mm;
200      
201       const index_t* ptr_main = Paso_SparseMatrix_borrowMainDiagonalPointer(A_p);
202      
203       /* forward substitution */
204      
205       if (n_block==1) {
206          for (i = 0; i < n; ++i) {
207         mm=ptr_main[i];
208         /* x_i=x_i-a_ik*x_k  (with k<i) */                    
209         for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<mm; ++iptr_ik) {
210            k=A_p->pattern->index[iptr_ik];                          
211            Paso_BlockOps_SMV_1(&x[i], &A_p->val[iptr_ik], &x[k]);
212         }
213         Paso_BlockOps_MV_1(&x[i], &A_p->val[i], &x[i]);
214          }
215       } else if (n_block==2) {
216          for (i = 0; i < n; ++i) {
217         mm=ptr_main[i];
218        
219         for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<mm; ++iptr_ik) {
220            k=A_p->pattern->index[iptr_ik];                          
221            Paso_BlockOps_SMV_2(&x[2*i], &A_p->val[4*iptr_ik], &x[2*k]);
222         }
223         Paso_BlockOps_MV_2(&x[2*i], &A_p->val[4*i], &x[2*i]);
224          }
225       } else if (n_block==3) {
226          for (i = 0; i < n; ++i) {
227         mm=ptr_main[i];
228         /* x_i=x_i-a_ik*x_k */
229         for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<mm; ++iptr_ik) {
230            k=A_p->pattern->index[iptr_ik];
231            Paso_BlockOps_SMV_3(&x[3*i], &A_p->val[9*iptr_ik], &x[3*k]);
232         }
233         Paso_BlockOps_MV_3(&x[3*i], &A_p->val[9*i], &x[3*i]);
234          }
235          
236       } /* add block size >3 */
237      
238        
239      
240       /* backward substitution */
241      
242       if (n_block==1) {
243          for (i = n-2; i > -1; ++i) {        
244         mm=ptr_main[i];
245         Paso_BlockOps_MV_1(&x[i], &A_p->val[mm], &x[i]);
246         for (iptr_ik=mm+1; iptr_ik < A_p->pattern->ptr[i+1]; ++iptr_ik) {
247            k=A_p->pattern->index[iptr_ik];  
248            Paso_BlockOps_SMV_1(&x[i], &A_p->val[iptr_ik], &x[k]);
249         }
250         Paso_BlockOps_MV_1(&x[i], &diag[i], &x[i]);
251          }
252       } else if (n_block==2) {
253          for (i = n-2; i > -1; ++i) {
254         mm=ptr_main[i];
255         Paso_BlockOps_MV_3(&x[2*i], &A_p->val[4*mm], &x[2*i]);
256         for (iptr_ik=mm+1; iptr_ik < A_p->pattern->ptr[i+1]; ++iptr_ik) {
257            k=A_p->pattern->index[iptr_ik];
258            Paso_BlockOps_SMV_2(&x[2*i], &A_p->val[4*iptr_ik], &x[2*k]);
259         }
260         Paso_BlockOps_MV_2(&x[2*i], &diag[i*4], &x[2*i]);      
261          }
262       } else if (n_block==3) {
263          for (i = n-2; i > -1; ++i) {
264        
265         mm=ptr_main[i];
266         Paso_BlockOps_MV_3(&x[3*i], &A_p->val[9*mm], &x[3*i]);
267        
268         for (iptr_ik=mm+1; iptr_ik < A_p->pattern->ptr[i+1]; ++iptr_ik) {
269            k=A_p->pattern->index[iptr_ik];    
270            Paso_BlockOps_SMV_3(&x[3*i], &A_p->val[9*iptr_ik], &x[3*k]);
271         }
272         Paso_BlockOps_MV_3(&x[3*i], &diag[i*9], &x[3*i]);
273          }
274          
275       } /* add block size >3 */      
276      
277       return;
278    }
279                
280  void Paso_Solver_solveLocalGS_colored(Paso_SparseMatrix* A_p, Paso_Solver_LocalGS * gs, double * x, const double * b)  void Paso_Solver_localGSSweep_colored(Paso_SparseMatrix* A_p, Paso_Solver_LocalGS * gs, double * x)
281  {  {
282     const dim_t n=A_p->numRows;     const dim_t n=A_p->numRows;
283     const dim_t n_block=A_p->row_block_size;     const dim_t n_block=A_p->row_block_size;
284     const double *diag = gs->diag;     const double *diag = gs->diag;
285     /* const index_t* pivot = gs->pivot;     index_t* pivot = gs->pivot;
286        const dim_t block_size=A_p->block_size;  use for block size >3*/     const dim_t block_size=A_p->block_size;
287        
288     register dim_t i,k;     register dim_t i,k;
289     register index_t color,iptr_ik, mm;     register index_t color,iptr_ik, mm;
    register double A11,A12,A21,A22,A13,A23,A33,A32,A31,S1,S2,S3,R1,R2,R3;  
290        
291     const index_t* coloring = Paso_Pattern_borrowColoringPointer(A_p->pattern);     const index_t* coloring = Paso_Pattern_borrowColoringPointer(A_p->pattern);
292     const dim_t num_colors = Paso_Pattern_getNumColors(A_p->pattern);     const dim_t num_colors = Paso_Pattern_getNumColors(A_p->pattern);
# Line 234  void Paso_Solver_solveLocalGS_colored(Pa Line 296  void Paso_Solver_solveLocalGS_colored(Pa
296        
297     /* color = 0 */     /* color = 0 */
298     if (n_block==1) {     if (n_block==1) {
299        #pragma omp parallel for schedule(static) private(i,S1, A11)        #pragma omp parallel for schedule(static) private(i)
300        for (i = 0; i < n; ++i) {        for (i = 0; i < n; ++i) {
301         if (coloring[i]==0) {       if (coloring[i]== 0 ) Paso_BlockOps_MV_1(&x[i], &diag[i], &x[i]);
302             /* x_i=x_i-a_ik*x_k */                            }
            S1=b[i];  
            A11=diag[i];  
            x[i]=A11*S1;  
         }  
      }  
303     } else if (n_block==2) {     } else if (n_block==2) {
304       #pragma omp parallel for schedule(static) private(i,S1,S2,A11,A21,A12,A22)           #pragma omp parallel for schedule(static) private(i)
305       for (i = 0; i < n; ++i) {       for (i = 0; i < n; ++i) {
306          if (coloring[i]== 0 ) {          if (coloring[i]== 0 ) Paso_BlockOps_MV_2(&x[2*i], &diag[i*4], &x[2*i]);
            /* x_i=x_i-a_ik*x_k */  
            S1=b[2*i];  
            S2=b[2*i+1];  
   
            A11=diag[i*4];  
            A12=diag[i*4+2];  
            A21=diag[i*4+1];  
            A22=diag[i*4+3];  
             
            x[2*i  ]=A11 * S1 + A12 * S2;  
            x[2*i+1]=A21 * S1 + A22 * S2;  
             
         }  
307       }       }
308        } else if (n_block==3) {      } else if (n_block==3) {
309       #pragma omp parallel for schedule(static) private(i,S1,S2,S3,A11,A21,A31,A12,A22,A32,A13,A23,A33)       #pragma omp parallel for schedule(static) private(i)
310       for (i = 0; i < n; ++i) {       for (i = 0; i < n; ++i) {
311          if (coloring[i]==0) {          if (coloring[i]==0) Paso_BlockOps_MV_3(&x[3*i], &diag[i*9], &x[3*i]);
            /* x_i=x_i-a_ik*x_k */  
            S1=b[3*i];  
            S2=b[3*i+1];  
            S3=b[3*i+2];  
            A11=diag[i*9  ];  
            A21=diag[i*9+1];  
            A31=diag[i*9+2];  
            A12=diag[i*9+3];  
            A22=diag[i*9+4];  
            A32=diag[i*9+5];  
            A13=diag[i*9+6];  
            A23=diag[i*9+7];  
            A33=diag[i*9+8];  
            x[3*i  ]=A11 * S1 + A12 * S2 + A13 * S3;  
            x[3*i+1]=A21 * S1 + A22 * S2 + A23 * S3;  
            x[3*i+2]=A31 * S1 + A32 * S2 + A33 * S3;  
         }  
312       }       }
313     } /* add block size >3 */     } else {
314          #pragma omp parallel for schedule(static) private(i)
315          for (i = 0; i < n; ++i) {
316         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]);
317          }
318       }
319        
320     for (color=1;color<num_colors;++color) {     for (color=1;color<num_colors;++color) {
321        if (n_block==1) {        if (n_block==1) {
322       #pragma omp parallel for schedule(static) private(i,iptr_ik,k,S1, A11,R1)       #pragma omp parallel for schedule(static) private(i,iptr_ik,k)
323       for (i = 0; i < n; ++i) {       for (i = 0; i < n; ++i) {
324          if (coloring[i]==color) {          if (coloring[i]==color) {
325             /* x_i=x_i-a_ik*x_k */                                 /* x_i=x_i-a_ik*x_k */                    
            S1=b[i];  
326             for (iptr_ik=A_p->pattern->ptr[i];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) {
327            k=A_p->pattern->index[iptr_ik];                                      k=A_p->pattern->index[iptr_ik];                          
328            if (coloring[k]<color) {            if (coloring[k]<color) Paso_BlockOps_SMV_1(&x[i], &A_p->val[iptr_ik], &x[k]);
              R1=x[k];                                
              S1-=A_p->val[iptr_ik]*R1;  
           }  
329             }             }
330             A11=diag[i];             Paso_BlockOps_MV_1(&x[i], &diag[i], &x[i]);
            x[i]=A11*S1;  
331          }          }
332       }       }
333        } else if (n_block==2) {        } else if (n_block==2) {
334       #pragma omp parallel for schedule(static) private(i,iptr_ik,k,S1,S2,R1,R2,A11,A21,A12,A22)       #pragma omp parallel for schedule(static) private(i,iptr_ik,k)
335       for (i = 0; i < n; ++i) {       for (i = 0; i < n; ++i) {
336          if (coloring[i]==color) {          if (coloring[i]==color) {
            /* x_i=x_i-a_ik*x_k */  
            S1=b[2*i];  
            S2=b[2*i+1];  
337             for (iptr_ik=A_p->pattern->ptr[i];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) {
338            k=A_p->pattern->index[iptr_ik];                                      k=A_p->pattern->index[iptr_ik];                          
339            if (coloring[k]<color) {            if (coloring[k]<color) Paso_BlockOps_SMV_2(&x[2*i], &A_p->val[4*iptr_ik], &x[2*k]);
              R1=x[2*k];  
              R2=x[2*k+1];  
              S1-=A_p->val[4*iptr_ik  ]*R1+A_p->val[4*iptr_ik+2]*R2;  
              S2-=A_p->val[4*iptr_ik+1]*R1+A_p->val[4*iptr_ik+3]*R2;  
           }  
340             }             }
341             A11=diag[i*4];             Paso_BlockOps_MV_2(&x[2*i], &diag[4*i], &x[2*i]);
            A12=diag[i*4+2];  
            A21=diag[i*4+1];  
            A22=diag[i*4+3];  
   
                x[2*i  ]=A11 * S1 + A12 * S2;  
                x[2*i+1]=A21 * S1 + A22 * S2;  
   
342          }          }
343                    
344       }       }
345        } else if (n_block==3) {        } else if (n_block==3) {
346       #pragma omp parallel for schedule(static) private(i,iptr_ik,k,S1,S2,S3,R1,R2,R3,A11,A21,A31,A12,A22,A32,A13,A23,A33)       #pragma omp parallel for schedule(static) private(i,iptr_ik,k)
347       for (i = 0; i < n; ++i) {       for (i = 0; i < n; ++i) {
348          if (coloring[i]==color) {          if (coloring[i]==color) {
            /* x_i=x_i-a_ik*x_k */  
            S1=b[3*i];  
            S2=b[3*i+1];  
            S3=b[3*i+2];  
349             for (iptr_ik=A_p->pattern->ptr[i];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) {
350            k=A_p->pattern->index[iptr_ik];                                      k=A_p->pattern->index[iptr_ik];                          
351            if (coloring[k]<color) {            if (coloring[k]<color) Paso_BlockOps_SMV_3(&x[3*i], &A_p->val[9*iptr_ik], &x[3*k]);
352               R1=x[3*k];             }
353               R2=x[3*k+1];             Paso_BlockOps_MV_3(&x[3*i], &diag[9*i], &x[3*i]);
354               R3=x[3*k+2];          }
355               S1-=A_p->val[9*iptr_ik  ]*R1+A_p->val[9*iptr_ik+3]*R2+A_p->val[9*iptr_ik+6]*R3;       }
356               S2-=A_p->val[9*iptr_ik+1]*R1+A_p->val[9*iptr_ik+4]*R2+A_p->val[9*iptr_ik+7]*R3;        } else {
357               S3-=A_p->val[9*iptr_ik+2]*R1+A_p->val[9*iptr_ik+5]*R2+A_p->val[9*iptr_ik+8]*R3;       #pragma omp parallel for schedule(static) private(i,iptr_ik,k)
358            }       for (i = 0; i < n; ++i) {
359            if (coloring[i] == color) {
360               for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<A_p->pattern->ptr[i+1]; ++iptr_ik) {
361              k=A_p->pattern->index[iptr_ik];                          
362              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]);
363             }             }
364             A11=diag[i*9  ];             Paso_BlockOps_Solve_N(n_block, &x[n_block*i], &diag[block_size*i], &pivot[n_block*i], &x[n_block*i]);
            A21=diag[i*9+1];  
            A31=diag[i*9+2];  
            A12=diag[i*9+3];  
            A22=diag[i*9+4];  
            A32=diag[i*9+5];  
            A13=diag[i*9+6];  
            A23=diag[i*9+7];  
            A33=diag[i*9+8];  
            x[3*i  ]=A11 * S1 + A12 * S2 + A13 * S3;  
            x[3*i+1]=A21 * S1 + A22 * S2 + A23 * S3;  
            x[3*i+2]=A31 * S1 + A32 * S2 + A33 * S3;  
365          }          }
366       }       }
367        } /* add block size >3 */        }
368     } /* end of coloring loop */     } /* end of coloring loop */
369        
370    
371     /* backward substitution */     /* backward substitution */
372     for (color=(num_colors)-2 ;color>-1;--color) { /* Note: color=(num_colors)-1 is not required */     for (color=(num_colors)-2 ;color>-1;--color) { /* Note: color=(num_colors)-1 is not required */
373        if (n_block==1) {        if (n_block==1) {
374       #pragma omp parallel for schedule(static) private(mm, i,iptr_ik,k,S1,R1)       #pragma omp parallel for schedule(static) private(mm, i,iptr_ik,k)
375       for (i = 0; i < n; ++i) {       for (i = 0; i < n; ++i) {
376          if (coloring[i]==color) {          if (coloring[i]==color) {
             
377             mm=ptr_main[i];             mm=ptr_main[i];
378             R1=x[i];             Paso_BlockOps_MV_1(&x[i], &A_p->val[mm], &x[i]);
            A11=A_p->val[mm];  
            S1 = A11 * R1;  
             
379             for (iptr_ik=A_p->pattern->ptr[i];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) {
380            k=A_p->pattern->index[iptr_ik];                                      k=A_p->pattern->index[iptr_ik];                          
381            if (coloring[k]>color) {            if (coloring[k]>color) Paso_BlockOps_SMV_1(&x[i], &A_p->val[iptr_ik], &x[k]);
              R1=x[k];  
              S1-=A_p->val[iptr_ik]*R1;  
           }  
382             }             }
383                         Paso_BlockOps_MV_1(&x[i], &diag[i], &x[i]);
            A11=diag[i];  
            x[i]=A11*S1;  
             
384          }          }
385       }       }
386        } else if (n_block==2) {        } else if (n_block==2) {
387       #pragma omp parallel for schedule(static) private(mm, i,iptr_ik,k,S1,S2,R1,R2,A11,A21,A12,A22)       #pragma omp parallel for schedule(static) private(mm, i,iptr_ik,k)
388       for (i = 0; i < n; ++i) {       for (i = 0; i < n; ++i) {
389          if (coloring[i]==color) {          if (coloring[i]==color) {
             
390             mm=ptr_main[i];             mm=ptr_main[i];
391                         Paso_BlockOps_MV_2(&x[2*i], &A_p->val[4*mm], &x[2*i]);
            R1=x[2*i];  
            R2=x[2*i+1];  
             
            A11=A_p->val[mm*4  ];  
            A21=A_p->val[mm*4+1];  
            A12=A_p->val[mm*4+2];  
            A22=A_p->val[mm*4+3];  
             
            S1 = A11 * R1 + A12 * R2;  
            S2 = A21 * R1 + A22 * R2;  
             
392             for (iptr_ik=A_p->pattern->ptr[i];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) {
393            k=A_p->pattern->index[iptr_ik];                                      k=A_p->pattern->index[iptr_ik];                          
394            if (coloring[k]>color) {            if (coloring[k]>color) Paso_BlockOps_SMV_2(&x[2*i], &A_p->val[4*iptr_ik], &x[2*k]);
              R1=x[2*k];  
              R2=x[2*k+1];  
              S1-=A_p->val[4*iptr_ik  ]*R1+A_p->val[4*iptr_ik+2]*R2;  
              S2-=A_p->val[4*iptr_ik+1]*R1+A_p->val[4*iptr_ik+3]*R2;  
           }  
395             }             }
396                         Paso_BlockOps_MV_2(&x[2*i], &diag[4*i], &x[2*i]);
            A11=diag[i*4];  
            A12=diag[i*4+2];  
            A21=diag[i*4+1];  
            A22=diag[i*4+3];  
             
            x[2*i  ]=A11 * S1 + A12 * S2;  
            x[2*i+1]=A21 * S1 + A22 * S2;  
             
397          }          }
398       }       }
399        } else if (n_block==3) {        } else if (n_block==3) {
400       #pragma omp parallel for schedule(static) private(mm, i,iptr_ik,k,S1,S2,S3,R1,R2,R3,A11,A21,A31,A12,A22,A32,A13,A23,A33)       #pragma omp parallel for schedule(static) private(mm, i,iptr_ik,k)
401       for (i = 0; i < n; ++i) {       for (i = 0; i < n; ++i) {
402          if (coloring[i]==color) {          if (coloring[i]==color) {
   
403             mm=ptr_main[i];             mm=ptr_main[i];
404             R1=x[3*i];             Paso_BlockOps_MV_3(&x[3*i], &A_p->val[9*mm], &x[3*i]);
            R2=x[3*i+1];  
            R3=x[3*i+2];  
             
            A11=A_p->val[mm*9  ];  
            A21=A_p->val[mm*9+1];  
            A31=A_p->val[mm*9+2];  
            A12=A_p->val[mm*9+3];  
            A22=A_p->val[mm*9+4];  
            A32=A_p->val[mm*9+5];  
            A13=A_p->val[mm*9+6];  
            A23=A_p->val[mm*9+7];  
            A33=A_p->val[mm*9+8];  
             
            S1 =A11 * R1 + A12 * R2 + A13 * R3;  
            S2 =A21 * R1 + A22 * R2 + A23 * R3;  
            S3 =A31 * R1 + A32 * R2 + A33 * R3;  
   
405             for (iptr_ik=A_p->pattern->ptr[i];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) {
406            k=A_p->pattern->index[iptr_ik];                                      k=A_p->pattern->index[iptr_ik];                          
407            if (coloring[k]>color) {            if (coloring[k]>color) Paso_BlockOps_SMV_3(&x[3*i], &A_p->val[9*iptr_ik], &x[3*k]);
              R1=x[3*k];  
              R2=x[3*k+1];  
              R3=x[3*k+2];  
              S1-=A_p->val[9*iptr_ik  ]*R1+A_p->val[9*iptr_ik+3]*R2+A_p->val[9*iptr_ik+6]*R3;  
              S2-=A_p->val[9*iptr_ik+1]*R1+A_p->val[9*iptr_ik+4]*R2+A_p->val[9*iptr_ik+7]*R3;  
              S3-=A_p->val[9*iptr_ik+2]*R1+A_p->val[9*iptr_ik+5]*R2+A_p->val[9*iptr_ik+8]*R3;  
           }  
408             }             }
409                         Paso_BlockOps_MV_3(&x[3*i], &diag[9*i], &x[3*i]);
            A11=diag[i*9  ];  
            A21=diag[i*9+1];  
            A31=diag[i*9+2];  
            A12=diag[i*9+3];  
            A22=diag[i*9+4];  
            A32=diag[i*9+5];  
            A13=diag[i*9+6];  
            A23=diag[i*9+7];  
            A33=diag[i*9+8];  
             
            x[3*i  ]=A11 * S1 + A12 * S2 + A13 * S3;  
            x[3*i+1]=A21 * S1 + A22 * S2 + A23 * S3;  
            x[3*i+2]=A31 * S1 + A32 * S2 + A33 * S3;  
             
410          }          }
411       }       }
412        } /* add block size >3 */              } else {
413         #pragma omp parallel for schedule(static) private(mm, i,iptr_ik,k)
414         for (i = 0; i < n; ++i) {
415            if (coloring[i]==color) {
416               mm=ptr_main[i];
417               Paso_BlockOps_MV_N( n_block, &x[n_block*i], &A_p->val[block_size*mm], &x[n_block*i] );
418               for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<A_p->pattern->ptr[i+1]; ++iptr_ik) {
419              k=A_p->pattern->index[iptr_ik];                          
420              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]);
421               }
422               Paso_BlockOps_Solve_N(n_block, &x[n_block*i], &diag[block_size*i], &pivot[n_block*i], &x[n_block*i]);
423            }
424         }
425          }
426     }     }
427     return;     return;
428  }  }
429    
 void Paso_Solver_solveLocalGS_sequential(Paso_SparseMatrix* A_p, Paso_Solver_LocalGS * gs, double * x, const double * b)  
 {  
       const dim_t n=A_p->numRows;  
       const dim_t n_block=A_p->row_block_size;  
       const double *diag = gs->diag;  
       /* const index_t* pivot = gs->pivot;  
       const dim_t block_size=A_p->block_size;  use for block size >3*/  
         
       register dim_t i,k;  
       register index_t iptr_ik, mm;  
       register double A11,A12,A21,A22,A13,A23,A33,A32,A31,S1,S2,S3,R1,R2,R3;  
         
       const index_t* ptr_main = Paso_SparseMatrix_borrowMainDiagonalPointer(A_p);  
         
       /* forward substitution */  
         
      if (n_block==1) {  
         for (i = 0; i < n; ++i) {  
           /* x_i=x_i-a_ik*x_k */                      
           S1=b[i];  
           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];                            
              if (k<i) {  
             R1=x[k];                                
             S1-=A_p->val[iptr_ik]*R1;  
              } else {  
             break; /* index is ordered */  
              }  
           }  
           A11=diag[i];  
           x[i]=A11*S1;  
            }  
      } else if (n_block==2) {  
         for (i = 0; i < n; ++i) {  
           S1=b[2*i];  
           S2=b[2*i+1];  
           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];                            
              if (k<i) {  
             R1=x[2*k];  
             R2=x[2*k+1];  
             S1-=A_p->val[4*iptr_ik  ]*R1+A_p->val[4*iptr_ik+2]*R2;  
             S2-=A_p->val[4*iptr_ik+1]*R1+A_p->val[4*iptr_ik+3]*R2;  
              } else {  
             break; /* index is ordered */  
              }  
           }  
           A11=diag[i*4];  
           A12=diag[i*4+2];  
           A21=diag[i*4+1];  
           A22=diag[i*4+3];  
             
           x[2*i  ]=A11 * S1 + A12 * S2;  
           x[2*i+1]=A21 * S1 + A22 * S2;  
             
         }  
      } else if (n_block==3) {  
         for (i = 0; i < n; ++i) {  
           /* x_i=x_i-a_ik*x_k */  
           S1=b[3*i];  
           S2=b[3*i+1];  
           S3=b[3*i+2];  
           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];                            
              if ( k<i ) {  
             R1=x[3*k];  
             R2=x[3*k+1];  
             R3=x[3*k+2];  
             S1-=A_p->val[9*iptr_ik  ]*R1+A_p->val[9*iptr_ik+3]*R2+A_p->val[9*iptr_ik+6]*R3;  
             S2-=A_p->val[9*iptr_ik+1]*R1+A_p->val[9*iptr_ik+4]*R2+A_p->val[9*iptr_ik+7]*R3;  
             S3-=A_p->val[9*iptr_ik+2]*R1+A_p->val[9*iptr_ik+5]*R2+A_p->val[9*iptr_ik+8]*R3;  
              } else {  
             break; /* index is ordered */  
              }  
           }  
           A11=diag[i*9  ];  
           A21=diag[i*9+1];  
           A31=diag[i*9+2];  
           A12=diag[i*9+3];  
           A22=diag[i*9+4];  
           A32=diag[i*9+5];  
           A13=diag[i*9+6];  
           A23=diag[i*9+7];  
           A33=diag[i*9+8];  
           x[3*i  ]=A11 * S1 + A12 * S2 + A13 * S3;  
           x[3*i+1]=A21 * S1 + A22 * S2 + A23 * S3;  
           x[3*i+2]=A31 * S1 + A32 * S2 + A33 * S3;  
            }  
         
      } /* add block size >3 */  
430    
         
         
       /* backward substitution */  
   
      if (n_block==1) {  
         for (i = n-2; i > -1; ++i) {  
             
           mm=ptr_main[i];  
           R1=x[i];  
           A11=A_p->val[mm];  
           S1 = A11 * R1;  
             
           for (iptr_ik=A_p->pattern->ptr[i+1]-1; iptr_ik > A_p->pattern->ptr[i]-1; ++iptr_ik) {  
              k=A_p->pattern->index[iptr_ik];                            
              if (k > i) {  
             R1=x[k];  
             S1-=A_p->val[iptr_ik]*R1;  
              } else {  
             break ;  
              }  
           }  
             
           A11=diag[i];  
           x[i]=A11*S1;  
             
         }  
      } else if (n_block==2) {  
         for (i = n-2; i > -1; ++i) {  
             
           mm=ptr_main[i];  
             
           R1=x[2*i];  
           R2=x[2*i+1];  
             
           A11=A_p->val[mm*4  ];  
           A21=A_p->val[mm*4+1];  
           A12=A_p->val[mm*4+2];  
           A22=A_p->val[mm*4+3];  
             
           S1 = A11 * R1 + A12 * R2;  
           S2 = A21 * R1 + A22 * R2;  
             
           for (iptr_ik=A_p->pattern->ptr[i+1]-1; iptr_ik > A_p->pattern->ptr[i]-1; ++iptr_ik) {  
              k=A_p->pattern->index[iptr_ik];                            
              if (k > i) {  
             R1=x[2*k];  
             R2=x[2*k+1];  
             S1-=A_p->val[4*iptr_ik  ]*R1+A_p->val[4*iptr_ik+2]*R2;  
             S2-=A_p->val[4*iptr_ik+1]*R1+A_p->val[4*iptr_ik+3]*R2;  
              } else {  
             break ;  
              }  
           }  
             
           A11=diag[i*4];  
           A12=diag[i*4+2];  
           A21=diag[i*4+1];  
           A22=diag[i*4+3];  
             
           x[2*i  ]=A11 * S1 + A12 * S2;  
           x[2*i+1]=A21 * S1 + A22 * S2;  
             
           
         }  
      } else if (n_block==3) {  
         for (i = n-2; i > -1; ++i) {  
             
           mm=ptr_main[i];  
           R1=x[3*i];  
           R2=x[3*i+1];  
           R3=x[3*i+2];  
             
           A11=A_p->val[mm*9  ];  
           A21=A_p->val[mm*9+1];  
           A31=A_p->val[mm*9+2];  
           A12=A_p->val[mm*9+3];  
           A22=A_p->val[mm*9+4];  
           A32=A_p->val[mm*9+5];  
           A13=A_p->val[mm*9+6];  
           A23=A_p->val[mm*9+7];  
           A33=A_p->val[mm*9+8];  
             
           S1 =A11 * R1 + A12 * R2 + A13 * R3;  
           S2 =A21 * R1 + A22 * R2 + A23 * R3;  
           S3 =A31 * R1 + A32 * R2 + A33 * R3;  
             
           for (iptr_ik=A_p->pattern->ptr[i+1]-1; iptr_ik > A_p->pattern->ptr[i]-1; ++iptr_ik) {  
              k=A_p->pattern->index[iptr_ik];                            
              if (k > i) {  
             R1=x[3*k];  
             R2=x[3*k+1];  
             R3=x[3*k+2];  
             S1-=A_p->val[9*iptr_ik  ]*R1+A_p->val[9*iptr_ik+3]*R2+A_p->val[9*iptr_ik+6]*R3;  
             S2-=A_p->val[9*iptr_ik+1]*R1+A_p->val[9*iptr_ik+4]*R2+A_p->val[9*iptr_ik+7]*R3;  
             S3-=A_p->val[9*iptr_ik+2]*R1+A_p->val[9*iptr_ik+5]*R2+A_p->val[9*iptr_ik+8]*R3;  
              } else {  
             break ;  
              }  
           }  
             
           A11=diag[i*9  ];  
           A21=diag[i*9+1];  
           A31=diag[i*9+2];  
           A12=diag[i*9+3];  
           A22=diag[i*9+4];  
           A32=diag[i*9+5];  
           A13=diag[i*9+6];  
           A23=diag[i*9+7];  
           A33=diag[i*9+8];  
             
           x[3*i  ]=A11 * S1 + A12 * S2 + A13 * S3;  
           x[3*i+1]=A21 * S1 + A22 * S2 + A23 * S3;  
           x[3*i+2]=A31 * S1 + A32 * S2 + A33 * S3;  
             
            }  
           
      } /* add block size >3 */        
   
       return;  
 }  
431    

Legend:
Removed from v.3096  
changed lines
  Added in v.3097

  ViewVC Help
Powered by ViewVC 1.1.26