/[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

trunk/paso/src/Solver_GS.c revision 2998 by artak, Wed Mar 24 04:00:20 2010 UTC trunk/paso/src/Smoother.c revision 3402 by gross, Tue Dec 7 07:36: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  /**************************************************************/  /**************************************************************/
24    
25  #include "Paso.h"  #include "Paso.h"
26  #include "Solver.h"  #include "Preconditioner.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 Smoother                                */
36    
37  void Paso_Solver_GS_free(Paso_Solver_GS * in) {  void Paso_Preconditioner_Smoother_free(Paso_Preconditioner_Smoother * in) {
38       if (in!=NULL) {       if (in!=NULL) {
39          MEMFREE(in->colorOf);      Paso_Preconditioner_LocalSmoother_free(in->localSmoother);
         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_Preconditioner_LocalSmoother_free(Paso_Preconditioner_LocalSmoother * 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_Preconditioner_Smoother* Paso_Preconditioner_Smoother_alloc(Paso_SystemMatrix * A_p, const bool_t jacobi, const bool_t is_local, const 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_Preconditioner_Smoother* out=MEMALLOC(1,Paso_Preconditioner_Smoother);
61    if (Paso_checkPtr(out)) return NULL;    if (! Esys_checkPtr(out)) {
62    out->colorOf=MEMALLOC(n,index_t);       out->localSmoother=Paso_Preconditioner_LocalSmoother_alloc(A_p->mainBlock,jacobi,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 (Esys_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_Preconditioner_Smoother_free(out);
69       return NULL;       return NULL;
70    }    }
71  }  }
72    Paso_Preconditioner_LocalSmoother* Paso_Preconditioner_LocalSmoother_alloc(Paso_SparseMatrix * A_p, const bool_t jacobi, bool_t verbose)
73    {
74      
75       dim_t n=A_p->numRows;
76       dim_t n_block=A_p->row_block_size;
77       dim_t block_size=A_p->block_size;
78      
79       double time0=Esys_timer();
80       /* allocations: */  
81       Paso_Preconditioner_LocalSmoother* out=MEMALLOC(1,Paso_Preconditioner_LocalSmoother);
82       if (! Esys_checkPtr(out)) {
83          
84          out->diag=MEMALLOC( ((size_t) n) * ((size_t) block_size),double);
85          out->pivot=MEMALLOC( ((size_t) n) * ((size_t)  n_block), index_t);
86          out->buffer=MEMALLOC( ((size_t) n) * ((size_t)  n_block), double);
87          out->Jacobi=jacobi;
88          
89          if ( ! ( Esys_checkPtr(out->diag) || Esys_checkPtr(out->pivot) ) ) {
90         Paso_SparseMatrix_invMain(A_p, out->diag, out->pivot );
91          }
92          
93       }
94       time0=Esys_timer()-time0;
95      
96       if (Esys_noError()) {
97          return out;
98       } else {
99          Paso_Preconditioner_LocalSmoother_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 Smoother 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_Preconditioner_Smoother_solve(Paso_SystemMatrix* A_p, Paso_Preconditioner_Smoother * smoother, double * x, const double * b,
118       register dim_t i,k;                      const dim_t sweeps, const bool_t x_is_initial)
119       register index_t color,iptr_ik,iptr_main;  {
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 = smoother->localSmoother->buffer;
123       index_t* pivot=NULL;     dim_t nsweeps=sweeps;
124           if (smoother->is_local) {
125       /* copy x into b*/        Paso_Preconditioner_LocalSmoother_solve(A_p->mainBlock,smoother->localSmoother,x,b,sweeps,x_is_initial);
126       #pragma omp parallel for private(i) schedule(static)     } else {
127       for (i=0;i<n*n_block;++i) x[i]=b[i];        if (! x_is_initial) {
128       /* forward substitution */       Paso_Copy(n, x, b);
129       for (color=0;color<gs->num_colors;++color) {       Paso_Preconditioner_LocalSmoother_Sweep(A_p->mainBlock,smoother->localSmoother,x);
130             if (n_block==1) {       nsweeps--;
131                #pragma omp parallel for schedule(static) private(i,iptr_ik,k,S1,R1,iptr_main)        }
132                for (i = 0; i < n; ++i) {        while (nsweeps > 0 ) {
133                     if (gs->colorOf[i]==color) {       Paso_Copy(n, b_new, b);
134                       /* x_i=x_i-a_ik*x_k */                               Paso_SystemMatrix_MatrixVector((-1.), A_p, x, 1., b_new); /* b_new = b - A*x */
135                       S1=x[i];       Paso_Preconditioner_LocalSmoother_Sweep(A_p->mainBlock,smoother->localSmoother,b_new);  
136                       for (iptr_ik=gs->pattern->ptr[i];iptr_ik<gs->pattern->ptr[i+1]; ++iptr_ik) {       Paso_AXPY(n, x, 1., b_new);
137                            k=gs->pattern->index[iptr_ik];                                 nsweeps--;
138                            if (gs->colorOf[k]<color) {        }
139                               R1=x[k];                                      
140                               S1-=gs->factors->val[iptr_ik]*R1;     }
141                             }  }
142                       }  void Paso_Preconditioner_LocalSmoother_solve(Paso_SparseMatrix* A_p, Paso_Preconditioner_LocalSmoother * smoother, double * x, const double * b,
143                       iptr_main=gs->main_iptr[i];                           const dim_t sweeps, const bool_t x_is_initial)
144                       x[i]=(1/gs->factors->val[iptr_main])*S1;  {
145                     }     const dim_t n= (A_p->numRows) * (A_p->row_block_size);
146                }     double *b_new = smoother->buffer;
147             } else if (n_block==2) {     dim_t nsweeps=sweeps;
148                #pragma omp parallel for schedule(static) private(i,iptr_ik,k,iptr_main,S1,S2,R1,R2,A11,A21,A12,A22,D,S11,S21,S12,S22)    
149                for (i = 0; i < n; ++i) {     if (! x_is_initial) {
150                     if (gs->colorOf[i]==color) {        Paso_Copy(n, x, b);
151                       /* x_i=x_i-a_ik*x_k */        Paso_Preconditioner_LocalSmoother_Sweep(A_p,smoother,x);
152                       S1=x[2*i];        nsweeps--;
153                       S2=x[2*i+1];     }
154                       for (iptr_ik=gs->pattern->ptr[i];iptr_ik<gs->pattern->ptr[i+1]; ++iptr_ik) {    
155                            k=gs->pattern->index[iptr_ik];                               while (nsweeps > 0 ) {
156                            if (gs->colorOf[k]<color) {       Paso_Copy(n, b_new, b);
157                               R1=x[2*k];  
158                               R2=x[2*k+1];       Paso_SparseMatrix_MatrixVector_CSR_OFFSET0((-1.), A_p, x, 1., b_new); /* b_new = b - A*x */
159                               S1-=gs->factors->val[4*iptr_ik  ]*R1+gs->factors->val[4*iptr_ik+2]*R2;       Paso_Preconditioner_LocalSmoother_Sweep(A_p,smoother,b_new);
160                               S2-=gs->factors->val[4*iptr_ik+1]*R1+gs->factors->val[4*iptr_ik+3]*R2;       Paso_AXPY(n, x, 1., b_new);
161                            }       nsweeps--;
162                       }     }
163                       iptr_main=gs->main_iptr[i];  }
164                       A11=gs->factors->val[iptr_main*4];  
165                       A21=gs->factors->val[iptr_main*4+1];  void Paso_Preconditioner_LocalSmoother_Sweep(Paso_SparseMatrix* A, Paso_Preconditioner_LocalSmoother * smoother, double * x)
166                       A12=gs->factors->val[iptr_main*4+2];  {
167                       A22=gs->factors->val[iptr_main*4+3];     const dim_t nt=omp_get_max_threads();
168                       D = A11*A22-A12*A21;     if (smoother->Jacobi) {
169                       if (ABS(D)>0.) {        Paso_BlockOps_solveAll(A->row_block_size,A->numRows,smoother->diag,smoother->pivot,x);
170                            D=1./D;     } else {
171                            S11= A22*D;        if (nt < 2) {
172                            S21=-A21*D;       Paso_Preconditioner_LocalSmoother_Sweep_sequential(A,smoother,x);
173                            S12=-A12*D;        } else {
174                            S22= A11*D;       Paso_Preconditioner_LocalSmoother_Sweep_colored(A,smoother,x);
175                            x[2*i  ]=S11*S1+S12*S2;        }
176                            x[2*i+1]=S21*S1+S22*S2;     }
                      } else {  
                             Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getGS: non-regular main diagonal block.");  
                        }  
                    }  
   
               }  
            } else if (n_block==3) {  
               #pragma omp parallel for schedule(static) private(i,iptr_ik,iptr_main,k,S1,S2,S3,R1,R2,R3,A11,A21,A31,A12,A22,A32,A13,A23,A33,D,S11,S21,S31,S12,S22,S32,S13,S23,S33)  
               for (i = 0; i < n; ++i) {  
                    if (gs->colorOf[i]==color) {  
                      /* x_i=x_i-a_ik*x_k */  
                      S1=x[3*i];  
                      S2=x[3*i+1];  
                      S3=x[3*i+2];  
                      for (iptr_ik=gs->pattern->ptr[i];iptr_ik<gs->pattern->ptr[i+1]; ++iptr_ik) {  
                           k=gs->pattern->index[iptr_ik];                            
                           if (gs->colorOf[k]<color) {  
                              R1=x[3*k];  
                              R2=x[3*k+1];  
                              R3=x[3*k+2];  
                              S1-=gs->factors->val[9*iptr_ik  ]*R1+gs->factors->val[9*iptr_ik+3]*R2+gs->factors->val[9*iptr_ik+6]*R3;  
                              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;  
                              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;  
                           }  
                      }  
                      iptr_main=gs->main_iptr[i];  
                      A11=gs->factors->val[iptr_main*9  ];  
                      A21=gs->factors->val[iptr_main*9+1];  
                      A31=gs->factors->val[iptr_main*9+2];  
                      A12=gs->factors->val[iptr_main*9+3];  
                      A22=gs->factors->val[iptr_main*9+4];  
                      A32=gs->factors->val[iptr_main*9+5];  
                      A13=gs->factors->val[iptr_main*9+6];  
                      A23=gs->factors->val[iptr_main*9+7];  
                      A33=gs->factors->val[iptr_main*9+8];  
                      D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);  
                      if (ABS(D)>0.) {  
                           D=1./D;  
                           S11=(A22*A33-A23*A32)*D;  
                           S21=(A31*A23-A21*A33)*D;  
                           S31=(A21*A32-A31*A22)*D;  
                           S12=(A13*A32-A12*A33)*D;  
                           S22=(A11*A33-A31*A13)*D;  
                           S32=(A12*A31-A11*A32)*D;  
                           S13=(A12*A23-A13*A22)*D;  
                           S23=(A13*A21-A11*A23)*D;  
                           S33=(A11*A22-A12*A21)*D;  
                              /* a_ik=a_ii^{-1}*a_ik */  
                           x[3*i  ]=S11*S1+S12*S2+S13*S3;  
                           x[3*i+1]=S21*S1+S22*S2+S23*S3;  
                           x[3*i+2]=S31*S1+S32*S2+S33*S3;  
                        } else {  
                             Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getGS: non-regular main diagonal block.");  
                        }  
                 }  
               }  
            }  
      }  
       
      /* Multipling with diag(A) */  
      Paso_Solver_applyBlockDiagonalMatrix(gs->n_block,gs->n,gs->diag,pivot,x,x);  
   
      /* backward substitution */  
      for (color=(gs->num_colors)-1;color>-1;--color) {  
            if (n_block==1) {  
               /*#pragma omp parallel for schedule(static) private(i,iptr_ik,k,S1,R1,iptr_main)*/  
               for (i = 0; i < n; ++i) {  
                    if (gs->colorOf[i]==color) {  
                      /* x_i=x_i-a_ik*x_k */  
                      S1=x[i];  
                      for (iptr_ik=gs->pattern->ptr[i];iptr_ik<gs->pattern->ptr[i+1]; ++iptr_ik) {  
                           k=gs->pattern->index[iptr_ik];                            
                           if (gs->colorOf[k]>color) {  
                              R1=x[k];  
                              S1-=gs->factors->val[iptr_ik]*R1;  
                           }  
                      }  
                      /*x[i]=S1;*/  
                      iptr_main=gs->main_iptr[i];  
                      x[i]=(1/gs->factors->val[iptr_main])*S1;  
                    }  
               }  
            } else if (n_block==2) {  
               #pragma omp parallel for schedule(static) private(i,iptr_ik,k,S1,S2,R1,R2,iptr_main,D,A11,A21,A12,A22,S11,S21,S12,S22)  
               for (i = 0; i < n; ++i) {  
                    if (gs->colorOf[i]==color) {  
                      /* x_i=x_i-a_ik*x_k */  
                      S1=x[2*i];  
                      S2=x[2*i+1];  
                      for (iptr_ik=gs->pattern->ptr[i];iptr_ik<gs->pattern->ptr[i+1]; ++iptr_ik) {  
                           k=gs->pattern->index[iptr_ik];                            
                           if (gs->colorOf[k]>color) {  
                              R1=x[2*k];  
                              R2=x[2*k+1];  
                              S1-=gs->factors->val[4*iptr_ik  ]*R1+gs->factors->val[4*iptr_ik+2]*R2;  
                              S2-=gs->factors->val[4*iptr_ik+1]*R1+gs->factors->val[4*iptr_ik+3]*R2;  
                           }  
                      }  
                      /*x[2*i]=S1;  
                      x[2*i+1]=S2;*/  
                      iptr_main=gs->main_iptr[i];  
                      A11=gs->factors->val[iptr_main*4];  
                      A21=gs->factors->val[iptr_main*4+1];  
                      A12=gs->factors->val[iptr_main*4+2];  
                      A22=gs->factors->val[iptr_main*4+3];  
                      D = A11*A22-A12*A21;  
                      if (ABS(D)>0.) {  
                           D=1./D;  
                           S11= A22*D;  
                           S21=-A21*D;  
                           S12=-A12*D;  
                           S22= A11*D;  
                           x[2*i  ]=S11*S1+S12*S2;  
                           x[2*i+1]=S21*S1+S22*S2;  
                      } else {  
                             Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getGS: non-regular main diagonal block.");  
                        }  
   
                     }  
               }  
            } else if (n_block==3) {  
               #pragma omp parallel for schedule(static) private(i,iptr_ik,k,S1,S2,S3,R1,R2,R3,iptr_main,D,A11,A21,A31,A12,A22,A32,A13,A23,A33,S11,S21,S31,S12,S22,S32,S13,S23,S33)  
               for (i = 0; i < n; ++i) {  
                    if (gs->colorOf[i]==color) {  
                      /* x_i=x_i-a_ik*x_k */  
                      S1=x[3*i  ];  
                      S2=x[3*i+1];  
                      S3=x[3*i+2];  
                      for (iptr_ik=gs->pattern->ptr[i];iptr_ik<gs->pattern->ptr[i+1]; ++iptr_ik) {  
                           k=gs->pattern->index[iptr_ik];                            
                           if (gs->colorOf[k]>color) {  
                              R1=x[3*k];  
                              R2=x[3*k+1];  
                              R3=x[3*k+2];  
                              S1-=gs->factors->val[9*iptr_ik  ]*R1+gs->factors->val[9*iptr_ik+3]*R2+gs->factors->val[9*iptr_ik+6]*R3;  
                              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;  
                              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;  
                           }  
                      }  
 /*                     x[3*i]=S1;  
                      x[3*i+1]=S2;  
                      x[3*i+2]=S3;  
 */                   iptr_main=gs->main_iptr[i];  
                      A11=gs->factors->val[iptr_main*9  ];  
                      A21=gs->factors->val[iptr_main*9+1];  
                      A31=gs->factors->val[iptr_main*9+2];  
                      A12=gs->factors->val[iptr_main*9+3];  
                      A22=gs->factors->val[iptr_main*9+4];  
                      A32=gs->factors->val[iptr_main*9+5];  
                      A13=gs->factors->val[iptr_main*9+6];  
                      A23=gs->factors->val[iptr_main*9+7];  
                      A33=gs->factors->val[iptr_main*9+8];  
                      D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);  
                      if (ABS(D)>0.) {  
                           D=1./D;  
                           S11=(A22*A33-A23*A32)*D;  
                           S21=(A31*A23-A21*A33)*D;  
                           S31=(A21*A32-A31*A22)*D;  
                           S12=(A13*A32-A12*A33)*D;  
                           S22=(A11*A33-A31*A13)*D;  
                           S32=(A12*A31-A11*A32)*D;  
                           S13=(A12*A23-A13*A22)*D;  
                           S23=(A13*A21-A11*A23)*D;  
                           S33=(A11*A22-A12*A21)*D;  
                           x[3*i  ]=S11*S1+S12*S2+S13*S3;  
                           x[3*i+1]=S21*S1+S22*S2+S23*S3;  
                           x[3*i+2]=S31*S1+S32*S2+S33*S3;  
                        } else {  
                             Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getGS: non-regular main diagonal block.");  
                        }  
                    }  
               }  
          }  
      }  
      return;  
177  }  }
178    
179    /* inplace Gauss-Seidel sweep in seqential mode: */
180    
181    void Paso_Preconditioner_LocalSmoother_Sweep_sequential(Paso_SparseMatrix* A_p, Paso_Preconditioner_LocalSmoother * smoother, double * x)
182    {
183       const dim_t n=A_p->numRows;
184       const dim_t n_block=A_p->row_block_size;
185       const double *diag = smoother->diag;
186       const index_t* pivot = smoother->pivot;
187       const dim_t block_len=A_p->block_size;
188      
189       register dim_t i,k;
190       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);
195       /* forward substitution */
196      
197       if (n_block==1) {
198          x[0]*=diag[0];
199          for (i = 1; i < n; ++i) {
200         mm=ptr_main[i];
201         /* 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) {
204            k=A_p->pattern->index[iptr_ik];
205            rtmp-=A_p->val[iptr_ik]*x[k];
206         }
207         x[i]=rtmp*diag[i];
208          }
209       } else if (n_block==2) {
210          Paso_BlockOps_MViP_2(&diag[0], &x[0]);
211          for (i = 1; i < n; ++i) {
212         mm=ptr_main[i];
213         for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<mm; ++iptr_ik) {
214            k=A_p->pattern->index[iptr_ik];                          
215            Paso_BlockOps_SMV_2(&x[2*i], &A_p->val[4*iptr_ik], &x[2*k]);
216         }
217         Paso_BlockOps_MViP_2(&diag[4*i], &x[2*i]);
218          }
219       } else if (n_block==3) {
220          Paso_BlockOps_MViP_3(&diag[0], &x[0]);
221          for (i = 1; i < n; ++i) {
222         mm=ptr_main[i];
223         /* x_i=x_i-a_ik*x_k */
224         for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<mm; ++iptr_ik) {
225            k=A_p->pattern->index[iptr_ik];
226            Paso_BlockOps_SMV_3(&x[3*i], &A_p->val[9*iptr_ik], &x[3*k]);
227         }
228         Paso_BlockOps_MViP_3(&diag[9*i], &x[3*i]);
229          }
230       } 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 */
243      
244       if (n_block==1) {
245          for (i = n-2; i > -1; --i) {        
246            mm=ptr_main[i];
247            rtmp=x[i]*A_p->val[mm];
248            for (iptr_ik=mm+1; iptr_ik < A_p->pattern->ptr[i+1]; ++iptr_ik) {
249               k=A_p->pattern->index[iptr_ik];
250               rtmp-=A_p->val[iptr_ik]*x[k];
251            }
252            x[i]=diag[i]*rtmp;
253          }
254          
255       } else if (n_block==2) {
256          for (i = n-2; i > -1; --i) {
257            mm=ptr_main[i];
258            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) {
260               k=A_p->pattern->index[iptr_ik];
261               Paso_BlockOps_SMV_2(&x[2*i], &A_p->val[4*iptr_ik], &x[2*k]);
262            }
263            Paso_BlockOps_MViP_2(&diag[i*4], &x[2*i]);
264            
265          }
266       } else if (n_block==3) {
267          for (i = n-2; i > -1; --i) {
268            mm=ptr_main[i];
269            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) {
271               k=A_p->pattern->index[iptr_ik];    
272               Paso_BlockOps_SMV_3(&x[3*i], &A_p->val[9*iptr_ik], &x[3*k]);
273            }
274            Paso_BlockOps_MViP_3(&diag[i*9], &x[3*i]);
275          }
276          
277       } 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;
297    }
298          
299    void Paso_Preconditioner_LocalSmoother_Sweep_colored(Paso_SparseMatrix* A_p, Paso_Preconditioner_LocalSmoother * smoother, double * x)
300    {
301       const dim_t n=A_p->numRows;
302       const dim_t n_block=A_p->row_block_size;
303       const double *diag = smoother->diag;
304       index_t* pivot = smoother->pivot;
305       const dim_t block_len=A_p->block_size;
306       double *y;
307      
308       register dim_t i,k;
309       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);
314       const dim_t num_colors = Paso_Pattern_getNumColors(A_p->pattern);
315       const index_t* ptr_main = Paso_SparseMatrix_borrowMainDiagonalPointer(A_p);
316    
317       #pragma omp parallel  private(mm, i,iptr_ik,k,rtmp, color, y)
318       {
319          if (n_block>3) {
320         y=TMPMEMALLOC(n_block, double);
321          } else {
322         y=NULL;
323          }
324          /* forward substitution */
325    
326          /* color = 0 */
327          if (n_block==1) {
328         #pragma omp  for schedule(static)
329         for (i = 0; i < n; ++i) {
330            if (coloring[i]== 0 ) x[i]*=diag[i];
331         }
332          } else if (n_block==2) {
333             #pragma omp for schedule(static)
334         for (i = 0; i < n; ++i) {
335            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         }
347          }
348      
349    
350          for (color=1;color<num_colors;++color) {
351    
352         if (n_block==1) {
353            #pragma omp for schedule(static)
354            for (i = 0; i < n; ++i) {
355               if (coloring[i]==color) {
356              /* x_i=x_i-a_ik*x_k */  
357              rtmp=x[i];
358              for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<A_p->pattern->ptr[i+1]; ++iptr_ik) {
359                 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               }
364            }
365         } else if (n_block==2) {
366            #pragma omp for schedule(static)
367            for (i = 0; i < n; ++i) {
368               if (coloring[i]==color) {
369              for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<A_p->pattern->ptr[i+1]; ++iptr_ik) {
370                 k=A_p->pattern->index[iptr_ik];                          
371                 if (coloring[k]<color) Paso_BlockOps_SMV_2(&x[2*i], &A_p->val[4*iptr_ik], &x[2*k]);
372              }
373              Paso_BlockOps_MViP_2(&diag[4*i], &x[2*i]);
374               }
375            }
376         } else if (n_block==3) {
377            #pragma omp for schedule(static)
378            for (i = 0; i < n; ++i) {
379               if (coloring[i]==color) {
380              for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<A_p->pattern->ptr[i+1]; ++iptr_ik) {
381                 k=A_p->pattern->index[iptr_ik];                          
382                 if (coloring[k]<color) Paso_BlockOps_SMV_3(&x[3*i], &A_p->val[9*iptr_ik], &x[3*k]);
383              }
384              Paso_BlockOps_MViP_3(&diag[9*i], &x[3*i]);
385               }
386            }
387         } else {
388            #pragma omp for schedule(static)
389            for (i = 0; i < n; ++i) {
390               if (coloring[i] == color) {
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_N(n_block, &x[n_block*i], &A_p->val[block_len*iptr_ik], &x[n_block*k]);
394              }
395              Paso_BlockOps_solve_N(n_block, &x[n_block*i], &diag[block_len*i], &pivot[n_block*i], &failed);
396               }
397            }
398         }
399          } /* end of coloring loop */
400          
401          /* backward substitution */
402          for (color=(num_colors)-2 ;color>-1;--color) { /* Note: color=(num_colors)-1 is not required */
403         if (n_block==1) {
404            #pragma omp for schedule(static)
405            for (i = 0; i < n; ++i) {
406               if (coloring[i]==color) {
407              mm=ptr_main[i];
408              rtmp=A_p->val[mm]*x[i];
409              for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<A_p->pattern->ptr[i+1]; ++iptr_ik) {
410                 k=A_p->pattern->index[iptr_ik];                          
411                 if (coloring[k]>color) rtmp-=A_p->val[iptr_ik]*x[k];
412              }
413              x[i]= rtmp*diag[i];
414               }
415            }
416         } else if (n_block==2) {
417            #pragma omp for schedule(static)
418            for (i = 0; i < n; ++i) {
419               if (coloring[i]==color) {
420              mm=ptr_main[i];
421              Paso_BlockOps_MViP_2(&A_p->val[4*mm], &x[2*i]);
422              for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<A_p->pattern->ptr[i+1]; ++iptr_ik) {
423                 k=A_p->pattern->index[iptr_ik];                          
424                 if (coloring[k]>color) Paso_BlockOps_SMV_2(&x[2*i], &A_p->val[4*iptr_ik], &x[2*k]);
425              }
426              Paso_BlockOps_MViP_2(&diag[4*i], &x[2*i]);
427               }
428            }
429         } else if (n_block==3) {
430            #pragma omp for schedule(static)
431            for (i = 0; i < n; ++i) {
432               if (coloring[i]==color) {
433              mm=ptr_main[i];
434              Paso_BlockOps_MViP_3(&A_p->val[9*mm], &x[3*i]);
435              for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<A_p->pattern->ptr[i+1]; ++iptr_ik) {
436                 k=A_p->pattern->index[iptr_ik];                          
437                 if (coloring[k]>color) Paso_BlockOps_SMV_3(&x[3*i], &A_p->val[9*iptr_ik], &x[3*k]);
438              }
439              Paso_BlockOps_MViP_3(&diag[9*i], &x[3*i]);
440               }
441            }
442         } else {
443            #pragma omp for schedule(static)
444            for (i = 0; i < n; ++i) {
445               if (coloring[i]==color) {
446              mm=ptr_main[i];
447              Paso_BlockOps_MV_N(n_block, &y[0], &A_p->val[block_len*mm], &x[n_block*i]);
448              for (iptr_ik=A_p->pattern->ptr[i];iptr_ik<A_p->pattern->ptr[i+1]; ++iptr_ik) {
449                 k=A_p->pattern->index[iptr_ik];                          
450                 if (coloring[k]>color) Paso_BlockOps_SMV_N(n_block, &y[0], &A_p->val[block_len*iptr_ik], &x[n_block*k]);
451              }
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               }
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;
464    }

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

  ViewVC Help
Powered by ViewVC 1.1.26