/[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 2881 by jfenwick, Thu Jan 28 02:03:15 2010 UTC revision 3105 by gross, Wed Aug 25 04:37:52 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    
32    
33  /**************************************************************/  /**************************************************************/
34    
35  /* free all memory used by GS                                */  /* free all memory used by GS                                */
36    
37  void Paso_Solver_GS_free(Paso_Solver_GS * in) {  void Paso_Solver_GS_free(Paso_Solver_GS * in) {
38       if (in!=NULL) {       if (in!=NULL) {
39          MEMFREE(in->colorOf);      Paso_Solver_LocalGS_free(in->localGS);
         Paso_SparseMatrix_free(in->factors);  
         MEMFREE(in->diag);  
         MEMFREE(in->main_iptr);  
         Paso_Pattern_free(in->pattern);    
40          MEMFREE(in);          MEMFREE(in);
41       }       }
42  }  }
43    void Paso_Solver_LocalGS_free(Paso_Solver_LocalGS * in) {
44       if (in!=NULL) {
45          MEMFREE(in->diag);
46          MEMFREE(in->pivot);
47          MEMFREE(in->buffer);
48          MEMFREE(in);
49       }
50    }
51  /**************************************************************/  /**************************************************************/
52    
53  /*   constructs the incomplete block factorization of  /*   constructs the symmetric Gauss-Seidel preconditioner    
54    
55  */  */
56  Paso_Solver_GS* Paso_Solver_getGS(Paso_SparseMatrix * A,bool_t verbose) {  Paso_Solver_GS* Paso_Solver_getGS(Paso_SystemMatrix * A_p, dim_t sweeps, bool_t is_local, bool_t verbose)
57    dim_t n=A->numRows;  {
58    dim_t n_block=A->row_block_size;    
   dim_t block_size=A->block_size;  
   register index_t i,iptr_main=0,iPtr;  
   double time0=0,time_color=0,time_fac=0;  
59    /* allocations: */      /* allocations: */  
60    Paso_Solver_GS* out=MEMALLOC(1,Paso_Solver_GS);    Paso_Solver_GS* out=MEMALLOC(1,Paso_Solver_GS);
61    if (Paso_checkPtr(out)) return NULL;    if (! Paso_checkPtr(out)) {
62    out->colorOf=MEMALLOC(n,index_t);       out->localGS=Paso_Solver_getLocalGS(A_p->mainBlock,sweeps,verbose);
63    out->diag=MEMALLOC( ((size_t) n) * ((size_t) block_size),double);       out->is_local=is_local;
   /*out->diag=MEMALLOC(A->len,double);*/  
   out->main_iptr=MEMALLOC(n,index_t);  
   out->pattern=Paso_Pattern_getReference(A->pattern);  
   out->factors=Paso_SparseMatrix_getReference(A);  
   out->n_block=n_block;  
   out->n=n;  
   
   if ( !(Paso_checkPtr(out->colorOf) || Paso_checkPtr(out->main_iptr) || Paso_checkPtr(out->factors)) ) {  
     time0=Paso_timer();  
     Paso_Pattern_color(A->pattern,&out->num_colors,out->colorOf);  
     time_color=Paso_timer()-time0;  
   
     if (Paso_noError()) {  
        time0=Paso_timer();  
   
           if (! (Paso_checkPtr(out->diag))) {  
              if (n_block==1) {  
                 #pragma omp parallel for private(i,iPtr,iptr_main) schedule(static)  
                 for (i = 0; i < A->pattern->numOutput; i++) {  
                    iptr_main=0;  
                    out->diag[i]=1.;  
                    /* find main diagonal */  
                    for (iPtr = A->pattern->ptr[i]; iPtr < A->pattern->ptr[i + 1]; iPtr++) {  
                        if (A->pattern->index[iPtr]==i) {  
                            iptr_main=iPtr;  
                            out->diag[i]=A->val[iPtr];  
                            break;  
                        }  
                    }  
                    out->main_iptr[i]=iptr_main;  
                 }  
              } else if (n_block==2) {  
                 #pragma omp parallel for private(i,iPtr,iptr_main) schedule(static)  
                 for (i = 0; i < A->pattern->numOutput; i++) {  
                    out->diag[i*4+0]= 1.;  
                    out->diag[i*4+1]= 0.;  
                    out->diag[i*4+2]= 0.;  
                    out->diag[i*4+3]= 1.;  
                    iptr_main=0;  
                    /* find main diagonal */  
                    for (iPtr = A->pattern->ptr[i]; iPtr < A->pattern->ptr[i + 1]; iPtr++) {  
                        if (A->pattern->index[iPtr]==i) {  
                              iptr_main=iPtr;  
                              out->diag[i*4]= A->val[iPtr*4];  
                              out->diag[i*4+1]=A->val[iPtr*4+1];  
                              out->diag[i*4+2]=A->val[iPtr*4+2];  
                              out->diag[i*4+3]= A->val[iPtr*4+3];  
                           break;  
                        }  
                    }  
                    out->main_iptr[i]=iptr_main;  
                 }    
              } else if (n_block==3) {  
                 #pragma omp parallel for private(i, iPtr,iptr_main) schedule(static)  
                 for (i = 0; i < A->pattern->numOutput; i++) {  
                    out->diag[i*9  ]=1.;  
                    out->diag[i*9+1]=0.;  
                    out->diag[i*9+2]=0.;  
                    out->diag[i*9+3]=0.;  
                    out->diag[i*9+4]=1.;  
                    out->diag[i*9+5]=0.;  
                    out->diag[i*9+6]=0.;  
                    out->diag[i*9+7]=0.;  
                    out->diag[i*9+8]=1.;  
                    iptr_main=0;  
                    /* find main diagonal */  
                    for (iPtr = A->pattern->ptr[i]; iPtr < A->pattern->ptr[i + 1]; iPtr++) {  
                        if (A->pattern->index[iPtr]==i) {  
                            iptr_main=iPtr;  
                            out->diag[i*9  ]=A->val[iPtr*9  ];  
                            out->diag[i*9+1]=A->val[iPtr*9+1];  
                            out->diag[i*9+2]=A->val[iPtr*9+2];  
                            out->diag[i*9+3]=A->val[iPtr*9+3];  
                            out->diag[i*9+4]=A->val[iPtr*9+4];  
                            out->diag[i*9+5]=A->val[iPtr*9+5];  
                            out->diag[i*9+6]=A->val[iPtr*9+6];  
                            out->diag[i*9+7]=A->val[iPtr*9+7];  
                            out->diag[i*9+8]=A->val[iPtr*9+8];  
                            break;  
                        }  
                    }  
                    out->main_iptr[i]=iptr_main;  
                 }  
              }  
            }  
   
      time_fac=Paso_timer()-time0;  
      }  
64    }    }
65    if (Paso_noError()) {    if (Paso_MPIInfo_noError(A_p->mpi_info)) {
       if (verbose) {  
          printf("GS: %d color used \n",out->num_colors);  
          printf("timing: GS: coloring/elemination : %e/%e\n",time_color,time_fac);  
      }  
66       return out;       return out;
67    } else  {    } else {
68       Paso_Solver_GS_free(out);       Paso_Solver_GS_free(out);
69       return NULL;       return NULL;
70    }    }
71  }  }
72    Paso_Solver_LocalGS* Paso_Solver_getLocalGS(Paso_SparseMatrix * A_p, dim_t sweeps, bool_t verbose) {
73      
74       dim_t n=A_p->numRows;
75       dim_t n_block=A_p->row_block_size;
76       dim_t block_size=A_p->block_size;
77      
78       double time0=Paso_timer();
79       /* allocations: */  
80       Paso_Solver_LocalGS* out=MEMALLOC(1,Paso_Solver_LocalGS);
81       if (! Paso_checkPtr(out)) {
82          
83          out->diag=MEMALLOC( ((size_t) n) * ((size_t) block_size),double);
84          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;
87          
88          if ( ! ( Paso_checkPtr(out->diag) || Paso_checkPtr(out->pivot) ) ) {
89         Paso_SparseMatrix_invMain(A_p, out->diag, out->pivot );
90          }
91          
92       }
93       time0=Paso_timer()-time0;
94      
95       if (Paso_noError()) {
96          if (verbose) printf("timing: Gauss-Seidel preparation: elemination : %e\n",time0);
97          return out;
98       } else {
99          Paso_Solver_LocalGS_free(out);
100          return NULL;
101       }
102    }
103    
104  /**************************************************************/  /*
105    
106    performs a few sweeps of the  from
107    
108  /* apply GS precondition b-> x                                S (x_{k} -  x_{k-1}) = b - A x_{k-1}
109    
110       in fact it solves LD^{-1}Ux=b in the form x= U^{-1} D L^{-1}b  where x_{0}=0 and S provides some approximatioon of A.
111    
112   should be called within a parallel region                                                Under MPI S is build on using A_p->mainBlock only.
113   barrier synconization should be performed to make sure that the input vector available  if GS is local the defect b - A x_{k-1} is calculated using A_p->mainBlock only.
114    
115  */  */
116    
117  void Paso_Solver_solveGS(Paso_Solver_GS * gs, double * x, double * b) {  void Paso_Solver_solveGS(Paso_SystemMatrix* A_p, Paso_Solver_GS * gs, double * x, const double * b)
118       register dim_t i,k;  {
119       register index_t color,iptr_ik,iptr_main;     register dim_t i;
120       register double A11,A12,A21,A22,A13,A23,A33,A32,A31,S1,S2,S3,R1,R2,R3,D,S21,S22,S12,S11,S13,S23,S33,S32,S31;     const dim_t n= (A_p->mainBlock->numRows) * (A_p->mainBlock->row_block_size);
121       dim_t n_block=gs->n_block;    
122       dim_t n=gs->n;     double *b_new = gs->localGS->buffer;
123       index_t* pivot=NULL;     dim_t sweeps=gs->localGS->sweeps;
124          
125       /* copy x into b*/     if (gs->is_local) {
126       #pragma omp parallel for private(i) schedule(static)        Paso_Solver_solveLocalGS(A_p->mainBlock,gs->localGS,x,b);
127       for (i=0;i<n*n_block;++i) x[i]=b[i];     } else {
128       /* forward substitution */        #pragma omp parallel for private(i) schedule(static)
129       for (color=0;color<gs->num_colors;++color) {        for (i=0;i<n;++i) x[i]=b[i];
130             if (n_block==1) {        
131                #pragma omp parallel for schedule(static) private(i,iptr_ik,k,S1,R1,iptr_main)        Paso_Solver_localGSSweep(A_p->mainBlock,gs->localGS,x);
132                for (i = 0; i < n; ++i) {        
133                     if (gs->colorOf[i]==color) {        while (sweeps > 1 ) {
134                       /* x_i=x_i-a_ik*x_k */                           #pragma omp parallel for private(i) schedule(static)
135                       S1=x[i];       for (i=0;i<n;++i) b_new[i]=b[i];
136                       for (iptr_ik=gs->pattern->ptr[i];iptr_ik<gs->pattern->ptr[i+1]; ++iptr_ik) {  
137                            k=gs->pattern->index[iptr_ik];                                     Paso_SystemMatrix_MatrixVector((-1.), A_p, x, 1., b_new); /* b_new = b - A*x */
138                            if (gs->colorOf[k]<color) {      
139                               R1=x[k];                                     Paso_Solver_localGSSweep(A_p->mainBlock,gs->localGS,b_new);
140                               S1-=gs->factors->val[iptr_ik]*R1;      
141                             }       #pragma omp parallel for private(i) schedule(static)
142                       }       for (i=0;i<n;++i) x[i]+=b_new[i];
143                       iptr_main=gs->main_iptr[i];       sweeps--;
144                       x[i]=(1/gs->factors->val[iptr_main])*S1;        }
145                     }        
146                }     }
147             } else if (n_block==2) {  }
148                #pragma omp parallel for schedule(static) private(i,iptr_ik,k,iptr_main,S1,S2,R1,R2)  void Paso_Solver_solveLocalGS(Paso_SparseMatrix* A_p, Paso_Solver_LocalGS * gs, double * x, const double * b)
149                for (i = 0; i < n; ++i) {  {
150                     if (gs->colorOf[i]==color) {     register dim_t i;
151                       /* x_i=x_i-a_ik*x_k */     const dim_t n= (A_p->numRows) * (A_p->row_block_size);
152                       S1=x[2*i];     double *b_new = gs->buffer;
153                       S2=x[2*i+1];     dim_t sweeps=gs->sweeps;
154                       for (iptr_ik=gs->pattern->ptr[i];iptr_ik<gs->pattern->ptr[i+1]; ++iptr_ik) {    
155                            k=gs->pattern->index[iptr_ik];                               #pragma omp parallel for private(i) schedule(static)
156                            if (gs->colorOf[k]<color) {     for (i=0;i<n;++i) x[i]=b[i];
157                               R1=x[2*k];    
158                               R2=x[2*k+1];     Paso_Solver_localGSSweep(A_p,gs,x);
159                               S1-=gs->factors->val[4*iptr_ik  ]*R1+gs->factors->val[4*iptr_ik+2]*R2;    
160                               S2-=gs->factors->val[4*iptr_ik+1]*R1+gs->factors->val[4*iptr_ik+3]*R2;     while (sweeps > 1 ) {
161                            }       #pragma omp parallel for private(i) schedule(static)
162                       }       for (i=0;i<n;++i) b_new[i]=b[i];
163                       iptr_main=gs->main_iptr[i];      
164                       A11=gs->factors->val[iptr_main*4];       Paso_SparseMatrix_MatrixVector_CSC_OFFSET0((-1.), A_p, x, 1., b_new); /* b_new = b - A*x */
165                       A21=gs->factors->val[iptr_main*4+1];      
166                       A12=gs->factors->val[iptr_main*4+2];       Paso_Solver_localGSSweep(A_p,gs,b_new);
167                       A22=gs->factors->val[iptr_main*4+3];      
168                       D = A11*A22-A12*A21;       #pragma omp parallel for private(i) schedule(static)
169                       if (ABS(D)>0.) {       for (i=0;i<n;++i) x[i]+=b_new[i];
170                            D=1./D;      
171                            S11= A22*D;       sweeps--;
172                            S21=-A21*D;     }
173                            S12=-A12*D;  }
174                            S22= A11*D;  
175                            x[2*i  ]=S11*S1+S12*S2;  void Paso_Solver_localGSSweep(Paso_SparseMatrix* A, Paso_Solver_LocalGS * gs, double * x)
176                            x[2*i+1]=S21*S1+S22*S2;  {
177                       } else {     #ifdef _OPENMP
178                              Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getGS: non-regular main diagonal block.");     const dim_t nt=omp_get_max_threads();
179                         }     #else
180                     }     const dim_t nt = 1;
181       #endif
182                }     if (nt < 2) {
183             } else if (n_block==3) {        Paso_Solver_localGSSweep_sequential(A,gs,x);
184                #pragma omp parallel for schedule(static) private(i,iptr_ik,iptr_main,k,S1,S2,S3,R1,R2,R3)     } else {
185                for (i = 0; i < n; ++i) {        Paso_Solver_localGSSweep_colored(A,gs,x);
186                     if (gs->colorOf[i]==color) {     }
187                       /* x_i=x_i-a_ik*x_k */  }
188                       S1=x[3*i];  
189                       S2=x[3*i+1];  /* inplace Gauss-Seidel sweep in seqential mode: */
190                       S3=x[3*i+2];  
191                       for (iptr_ik=gs->pattern->ptr[i];iptr_ik<gs->pattern->ptr[i+1]; ++iptr_ik) {  void Paso_Solver_localGSSweep_sequential(Paso_SparseMatrix* A_p, Paso_Solver_LocalGS * gs, double * x)
192                            k=gs->pattern->index[iptr_ik];                            {
193                            if (gs->colorOf[k]<color) {     const dim_t n=A_p->numRows;
194                               R1=x[3*k];     const dim_t n_block=A_p->row_block_size;
195                               R2=x[3*k+1];     const double *diag = gs->diag;
196                               R3=x[3*k+2];     /* const index_t* pivot = gs->pivot;
197                               S1-=gs->factors->val[9*iptr_ik  ]*R1+gs->factors->val[9*iptr_ik+3]*R2+gs->factors->val[9*iptr_ik+6]*R3;     const dim_t block_size=A_p->block_size;  use for block size >3*/
198                               S2-=gs->factors->val[9*iptr_ik+1]*R1+gs->factors->val[9*iptr_ik+4]*R2+gs->factors->val[9*iptr_ik+7]*R3;    
199                               S3-=gs->factors->val[9*iptr_ik+2]*R1+gs->factors->val[9*iptr_ik+5]*R2+gs->factors->val[9*iptr_ik+8]*R3;     register dim_t i,k;
200                            }     register index_t iptr_ik, mm;
201                       }    
202                       iptr_main=gs->main_iptr[i];     const index_t* ptr_main = Paso_SparseMatrix_borrowMainDiagonalPointer(A_p);
203                       A11=gs->factors->val[iptr_main*9  ];     /* forward substitution */
204                       A21=gs->factors->val[iptr_main*9+1];    
205                       A31=gs->factors->val[iptr_main*9+2];     if (n_block==1) {
206                       A12=gs->factors->val[iptr_main*9+3];        Paso_BlockOps_MV_1(&x[0], &diag[0], &x[0]);
207                       A22=gs->factors->val[iptr_main*9+4];        for (i = 1; i < n; ++i) {
208                       A32=gs->factors->val[iptr_main*9+5];       mm=ptr_main[i];
209                       A13=gs->factors->val[iptr_main*9+6];       /* x_i=x_i-a_ik*x_k  (with k<i) */
210                       A23=gs->factors->val[iptr_main*9+7];       for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<mm; ++iptr_ik) {
211                       A33=gs->factors->val[iptr_main*9+8];          k=A_p->pattern->index[iptr_ik];  
212                       D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);          Paso_BlockOps_SMV_1(&x[i], &A_p->val[iptr_ik], &x[k]);
213                       if (ABS(D)>0.) {       }
214                            D=1./D;       Paso_BlockOps_MV_1(&x[i], &diag[i], &x[i]);
215                            S11=(A22*A33-A23*A32)*D;        }
216                            S21=(A31*A23-A21*A33)*D;     } else if (n_block==2) {
217                            S31=(A21*A32-A31*A22)*D;        Paso_BlockOps_MV_2(&x[0], &diag[0], &x[0]);
218                            S12=(A13*A32-A12*A33)*D;        for (i = 1; i < n; ++i) {
219                            S22=(A11*A33-A31*A13)*D;       mm=ptr_main[i];
220                            S32=(A12*A31-A11*A32)*D;       for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<mm; ++iptr_ik) {
221                            S13=(A12*A23-A13*A22)*D;          k=A_p->pattern->index[iptr_ik];                          
222                            S23=(A13*A21-A11*A23)*D;          Paso_BlockOps_SMV_2(&x[2*i], &A_p->val[4*iptr_ik], &x[2*k]);
223                            S33=(A11*A22-A12*A21)*D;       }
224                               /* a_ik=a_ii^{-1}*a_ik */       Paso_BlockOps_MV_2(&x[2*i], &diag[4*i], &x[2*i]);
225                            x[3*i  ]=S11*S1+S12*S2+S13*S3;        }
226                            x[3*i+1]=S21*S1+S22*S2+S23*S3;     } else if (n_block==3) {
227                            x[3*i+2]=S31*S1+S32*S2+S33*S3;        Paso_BlockOps_MV_3(&x[0], &diag[0], &x[0]);
228                         } else {        for (i = 1; i < n; ++i) {
229                              Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getGS: non-regular main diagonal block.");       mm=ptr_main[i];
230                         }       /* x_i=x_i-a_ik*x_k */
231                  }       for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<mm; ++iptr_ik) {
232                }          k=A_p->pattern->index[iptr_ik];
233             }          Paso_BlockOps_SMV_3(&x[3*i], &A_p->val[9*iptr_ik], &x[3*k]);
234       }       }
235             Paso_BlockOps_MV_3(&x[3*i], &diag[9*i], &x[3*i]);
236       /* Multipling with diag(A) */        }
237       Paso_Solver_applyBlockDiagonalMatrix(gs->n_block,gs->n,gs->diag,pivot,x,x);     } /* add block size >3 */
238    
239       /* backward substitution */     /* backward substitution */
240       for (color=(gs->num_colors)-1;color>-1;--color) {    
241             if (n_block==1) {     if (n_block==1) {
242                #pragma omp parallel for schedule(static) private(i,iptr_ik,k,S1,R1)        for (i = n-2; i > -1; --i) {        
243                for (i = 0; i < n; ++i) {          mm=ptr_main[i];
244                     if (gs->colorOf[i]==color) {          Paso_BlockOps_MV_1(&x[i], &A_p->val[mm], &x[i]);
245                       /* x_i=x_i-a_ik*x_k */          for (iptr_ik=mm+1; iptr_ik < A_p->pattern->ptr[i+1]; ++iptr_ik) {
246                       S1=x[i];             k=A_p->pattern->index[iptr_ik];  
247                       for (iptr_ik=gs->pattern->ptr[i];iptr_ik<gs->pattern->ptr[i+1]; ++iptr_ik) {             Paso_BlockOps_SMV_1(&x[i], &A_p->val[iptr_ik], &x[k]);
248                            k=gs->pattern->index[iptr_ik];                                    }
249                            if (gs->colorOf[k]>color) {          Paso_BlockOps_MV_1(&x[i], &diag[i], &x[i]);
250                               R1=x[k];        }
251                               S1-=gs->factors->val[iptr_ik]*R1;        
252                            }     } else if (n_block==2) {
253                       }        for (i = n-2; i > -1; --i) {
254                       /*x[i]=S1;*/          mm=ptr_main[i];
255                       iptr_main=gs->main_iptr[i];          Paso_BlockOps_MV_2(&x[2*i], &A_p->val[4*mm], &x[2*i]);
256                       x[i]=(1/gs->factors->val[iptr_main])*S1;          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             } else if (n_block==2) {          }
260                #pragma omp parallel for schedule(static) private(i,iptr_ik,k,S1,S2,R1,R2)          Paso_BlockOps_MV_2(&x[2*i], &diag[i*4], &x[2*i]);
261                for (i = 0; i < n; ++i) {        }
262                     if (gs->colorOf[i]==color) {     } else if (n_block==3) {
263                       /* x_i=x_i-a_ik*x_k */        for (i = n-2; i > -1; --i) {
264                       S1=x[2*i];          mm=ptr_main[i];
265                       S2=x[2*i+1];          Paso_BlockOps_MV_3(&x[3*i], &A_p->val[9*mm], &x[3*i]);
266                       for (iptr_ik=gs->pattern->ptr[i];iptr_ik<gs->pattern->ptr[i+1]; ++iptr_ik) {      
267                            k=gs->pattern->index[iptr_ik];                                    for (iptr_ik=mm+1; iptr_ik < A_p->pattern->ptr[i+1]; ++iptr_ik) {
268                            if (gs->colorOf[k]>color) {             k=A_p->pattern->index[iptr_ik];    
269                               R1=x[2*k];             Paso_BlockOps_SMV_3(&x[3*i], &A_p->val[9*iptr_ik], &x[3*k]);
270                               R2=x[2*k+1];          }
271                               S1-=gs->factors->val[4*iptr_ik  ]*R1+gs->factors->val[4*iptr_ik+2]*R2;          Paso_BlockOps_MV_3(&x[3*i], &diag[i*9], &x[3*i]);
272                               S2-=gs->factors->val[4*iptr_ik+1]*R1+gs->factors->val[4*iptr_ik+3]*R2;        }
273                            }        
274                       }     } /* add block size >3 */      
275                       /*x[2*i]=S1;    
276                       x[2*i+1]=S2;*/     return;
277                       iptr_main=gs->main_iptr[i];  }
278                       A11=gs->factors->val[iptr_main*4];        
279                       A21=gs->factors->val[iptr_main*4+1];  void Paso_Solver_localGSSweep_colored(Paso_SparseMatrix* A_p, Paso_Solver_LocalGS * gs, double * x)
280                       A12=gs->factors->val[iptr_main*4+2];  {
281                       A22=gs->factors->val[iptr_main*4+3];     const dim_t n=A_p->numRows;
282                       D = A11*A22-A12*A21;     const dim_t n_block=A_p->row_block_size;
283                       if (ABS(D)>0.) {     const double *diag = gs->diag;
284                            D=1./D;     index_t* pivot = gs->pivot;
285                            S11= A22*D;     const dim_t block_size=A_p->block_size;
286                            S21=-A21*D;    
287                            S12=-A12*D;     register dim_t i,k;
288                            S22= A11*D;     register index_t color,iptr_ik, mm;
289                            x[2*i  ]=S11*S1+S12*S2;    
290                            x[2*i+1]=S21*S1+S22*S2;     const index_t* coloring = Paso_Pattern_borrowColoringPointer(A_p->pattern);
291                       } else {     const dim_t num_colors = Paso_Pattern_getNumColors(A_p->pattern);
292                              Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getGS: non-regular main diagonal block.");     const index_t* ptr_main = Paso_SparseMatrix_borrowMainDiagonalPointer(A_p);
293                         }    
294       /* forward substitution */
295    
296      
297       /* color = 0 */
298       if (n_block==1) {
299          #pragma omp parallel for schedule(static) private(i)
300          for (i = 0; i < n; ++i) {
301         if (coloring[i]== 0 ) Paso_BlockOps_MV_1(&x[i], &diag[i], &x[i]);
302          }
303       } else if (n_block==2) {
304             #pragma omp parallel for schedule(static) private(i)
305         for (i = 0; i < n; ++i) {
306            if (coloring[i]== 0 ) Paso_BlockOps_MV_2(&x[2*i], &diag[i*4], &x[2*i]);
307         }
308        } else if (n_block==3) {
309         #pragma omp parallel for schedule(static) private(i)
310         for (i = 0; i < n; ++i) {
311            if (coloring[i]==0) Paso_BlockOps_MV_3(&x[3*i], &diag[i*9], &x[3*i]);
312         }
313       } 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) {
321    
322                      }        if (n_block==1) {
323                }       #pragma omp parallel for schedule(static) private(i,iptr_ik,k)
324             } else if (n_block==3) {       for (i = 0; i < n; ++i) {
325                #pragma omp parallel for schedule(static) private(i,iptr_ik,k,S1,S2,S3,R1,R2,R3)          if (coloring[i]==color) {
326                for (i = 0; i < n; ++i) {             /* x_i=x_i-a_ik*x_k */                    
327                     if (gs->colorOf[i]==color) {             for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<A_p->pattern->ptr[i+1]; ++iptr_ik) {
328                       /* x_i=x_i-a_ik*x_k */            k=A_p->pattern->index[iptr_ik];
329                       S1=x[3*i  ];            if (coloring[k]<color) Paso_BlockOps_SMV_1(&x[i], &A_p->val[iptr_ik], &x[k]);
330                       S2=x[3*i+1];             }
331                       S3=x[3*i+2];             Paso_BlockOps_MV_1(&x[i], &diag[i], &x[i]);
332                       for (iptr_ik=gs->pattern->ptr[i];iptr_ik<gs->pattern->ptr[i+1]; ++iptr_ik) {          }
333                            k=gs->pattern->index[iptr_ik];                                 }
334                            if (gs->colorOf[k]>color) {        } else if (n_block==2) {
335                               R1=x[3*k];       #pragma omp parallel for schedule(static) private(i,iptr_ik,k)
336                               R2=x[3*k+1];       for (i = 0; i < n; ++i) {
337                               R3=x[3*k+2];          if (coloring[i]==color) {
338                               S1-=gs->factors->val[9*iptr_ik  ]*R1+gs->factors->val[9*iptr_ik+3]*R2+gs->factors->val[9*iptr_ik+6]*R3;             for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<A_p->pattern->ptr[i+1]; ++iptr_ik) {
339                               S2-=gs->factors->val[9*iptr_ik+1]*R1+gs->factors->val[9*iptr_ik+4]*R2+gs->factors->val[9*iptr_ik+7]*R3;            k=A_p->pattern->index[iptr_ik];                          
340                               S3-=gs->factors->val[9*iptr_ik+2]*R1+gs->factors->val[9*iptr_ik+5]*R2+gs->factors->val[9*iptr_ik+8]*R3;            if (coloring[k]<color) Paso_BlockOps_SMV_2(&x[2*i], &A_p->val[4*iptr_ik], &x[2*k]);
341                            }             }
342                       }             Paso_BlockOps_MV_2(&x[2*i], &diag[4*i], &x[2*i]);
343  /*                     x[3*i]=S1;          }
344                       x[3*i+1]=S2;       }
345                       x[3*i+2]=S3;        } else if (n_block==3) {
346  */                   iptr_main=gs->main_iptr[i];       #pragma omp parallel for schedule(static) private(i,iptr_ik,k)
347                       A11=gs->factors->val[iptr_main*9  ];       for (i = 0; i < n; ++i) {
348                       A21=gs->factors->val[iptr_main*9+1];          if (coloring[i]==color) {
349                       A31=gs->factors->val[iptr_main*9+2];             for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<A_p->pattern->ptr[i+1]; ++iptr_ik) {
350                       A12=gs->factors->val[iptr_main*9+3];            k=A_p->pattern->index[iptr_ik];                          
351                       A22=gs->factors->val[iptr_main*9+4];            if (coloring[k]<color) Paso_BlockOps_SMV_3(&x[3*i], &A_p->val[9*iptr_ik], &x[3*k]);
352                       A32=gs->factors->val[iptr_main*9+5];             }
353                       A13=gs->factors->val[iptr_main*9+6];             Paso_BlockOps_MV_3(&x[3*i], &diag[9*i], &x[3*i]);
354                       A23=gs->factors->val[iptr_main*9+7];          }
355                       A33=gs->factors->val[iptr_main*9+8];       }
356                       D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);        } else {
357                       if (ABS(D)>0.) {       #pragma omp parallel for schedule(static) private(i,iptr_ik,k)
358                            D=1./D;       for (i = 0; i < n; ++i) {
359                            S11=(A22*A33-A23*A32)*D;          if (coloring[i] == color) {
360                            S21=(A31*A23-A21*A33)*D;             for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<A_p->pattern->ptr[i+1]; ++iptr_ik) {
361                            S31=(A21*A32-A31*A22)*D;            k=A_p->pattern->index[iptr_ik];                          
362                            S12=(A13*A32-A12*A33)*D;            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                            S22=(A11*A33-A31*A13)*D;             }
364                            S32=(A12*A31-A11*A32)*D;             Paso_BlockOps_Solve_N(n_block, &x[n_block*i], &diag[block_size*i], &pivot[n_block*i], &x[n_block*i]);
365                            S13=(A12*A23-A13*A22)*D;          }
366                            S23=(A13*A21-A11*A23)*D;       }
367                            S33=(A11*A22-A12*A21)*D;        }
368                            x[3*i  ]=S11*S1+S12*S2+S13*S3;     } /* end of coloring loop */
369                            x[3*i+1]=S21*S1+S22*S2+S23*S3;    
370                            x[3*i+2]=S31*S1+S32*S2+S33*S3;     /* backward substitution */
371                         } else {     for (color=(num_colors)-2 ;color>-1;--color) { /* Note: color=(num_colors)-1 is not required */
372                              Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getGS: non-regular main diagonal block.");        if (n_block==1) {
373                         }       #pragma omp parallel for schedule(static) private(mm, i,iptr_ik,k)
374                     }       for (i = 0; i < n; ++i) {
375                }          if (coloring[i]==color) {
376           }             mm=ptr_main[i];
377       }             Paso_BlockOps_MV_1(&x[i], &A_p->val[mm], &x[i]);
378       return;             for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<A_p->pattern->ptr[i+1]; ++iptr_ik) {
379              k=A_p->pattern->index[iptr_ik];                          
380              if (coloring[k]>color) Paso_BlockOps_SMV_1(&x[i], &A_p->val[iptr_ik], &x[k]);
381               }
382               Paso_BlockOps_MV_1(&x[i], &diag[i], &x[i]);
383            }
384         }
385          } else if (n_block==2) {
386         #pragma omp parallel for schedule(static) private(mm, i,iptr_ik,k)
387         for (i = 0; i < n; ++i) {
388            if (coloring[i]==color) {
389               mm=ptr_main[i];
390               Paso_BlockOps_MV_2(&x[2*i], &A_p->val[4*mm], &x[2*i]);
391               for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<A_p->pattern->ptr[i+1]; ++iptr_ik) {
392              k=A_p->pattern->index[iptr_ik];                          
393              if (coloring[k]>color) Paso_BlockOps_SMV_2(&x[2*i], &A_p->val[4*iptr_ik], &x[2*k]);
394               }
395               Paso_BlockOps_MV_2(&x[2*i], &diag[4*i], &x[2*i]);
396            }
397         }
398          } else if (n_block==3) {
399         #pragma omp parallel for schedule(static) private(mm, i,iptr_ik,k)
400         for (i = 0; i < n; ++i) {
401            if (coloring[i]==color) {
402               mm=ptr_main[i];
403               Paso_BlockOps_MV_3(&x[3*i], &A_p->val[9*mm], &x[3*i]);
404               for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<A_p->pattern->ptr[i+1]; ++iptr_ik) {
405              k=A_p->pattern->index[iptr_ik];                          
406              if (coloring[k]>color) Paso_BlockOps_SMV_3(&x[3*i], &A_p->val[9*iptr_ik], &x[3*k]);
407               }
408               Paso_BlockOps_MV_3(&x[3*i], &diag[9*i], &x[3*i]);
409            }
410         }
411          } else {
412         #pragma omp parallel for schedule(static) private(mm, i,iptr_ik,k)
413         for (i = 0; i < n; ++i) {
414            if (coloring[i]==color) {
415               mm=ptr_main[i];
416               Paso_BlockOps_MV_N( n_block, &x[n_block*i], &A_p->val[block_size*mm], &x[n_block*i] );
417               for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<A_p->pattern->ptr[i+1]; ++iptr_ik) {
418              k=A_p->pattern->index[iptr_ik];                          
419              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]);
420               }
421               Paso_BlockOps_Solve_N(n_block, &x[n_block*i], &diag[block_size*i], &pivot[n_block*i], &x[n_block*i]);
422            }
423         }
424          }
425       }
426       return;
427  }  }
428    
429    
430    

Legend:
Removed from v.2881  
changed lines
  Added in v.3105

  ViewVC Help
Powered by ViewVC 1.1.26