/[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 3120 by gross, Mon Aug 30 10:48:00 2010 UTC revision 3125 by gross, Tue Aug 31 05:37:30 2010 UTC
# Line 34  Line 34 
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_Preconditioner_GS_free(Paso_Preconditioner_GS * in) {
38       if (in!=NULL) {       if (in!=NULL) {
39      Paso_Solver_LocalGS_free(in->localGS);      Paso_Preconditioner_LocalGS_free(in->localGS);
40          MEMFREE(in);          MEMFREE(in);
41       }       }
42  }  }
43  void Paso_Solver_LocalGS_free(Paso_Solver_LocalGS * in) {  void Paso_Preconditioner_LocalGS_free(Paso_Preconditioner_LocalGS * in) {
44     if (in!=NULL) {     if (in!=NULL) {
45        MEMFREE(in->diag);        MEMFREE(in->diag);
46        MEMFREE(in->pivot);        MEMFREE(in->pivot);
# Line 53  void Paso_Solver_LocalGS_free(Paso_Solve Line 53  void Paso_Solver_LocalGS_free(Paso_Solve
53  /*   constructs the symmetric Gauss-Seidel preconditioner      /*   constructs the symmetric Gauss-Seidel preconditioner    
54    
55  */  */
56  Paso_Solver_GS* Paso_Solver_getGS(Paso_SystemMatrix * A_p, dim_t sweeps, bool_t is_local, bool_t verbose)  Paso_Preconditioner_GS* Paso_Preconditioner_GS_alloc(Paso_SystemMatrix * A_p, dim_t sweeps, bool_t is_local, bool_t verbose)
57  {  {
58        
59    /* allocations: */      /* allocations: */  
60    Paso_Solver_GS* out=MEMALLOC(1,Paso_Solver_GS);    Paso_Preconditioner_GS* out=MEMALLOC(1,Paso_Preconditioner_GS);
61    if (! Paso_checkPtr(out)) {    if (! Paso_checkPtr(out)) {
62       out->localGS=Paso_Solver_getLocalGS(A_p->mainBlock,sweeps,verbose);       out->localGS=Paso_Preconditioner_LocalGS_alloc(A_p->mainBlock,sweeps,verbose);
63       out->is_local=is_local;       out->is_local=is_local;
64    }    }
65    if (Paso_MPIInfo_noError(A_p->mpi_info)) {    if (Paso_MPIInfo_noError(A_p->mpi_info)) {
66       return out;       return out;
67    } else {    } else {
68       Paso_Solver_GS_free(out);       Paso_Preconditioner_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) {  Paso_Preconditioner_LocalGS* Paso_Preconditioner_LocalGS_alloc(Paso_SparseMatrix * A_p, dim_t sweeps, bool_t verbose) {
73        
74     dim_t n=A_p->numRows;     dim_t n=A_p->numRows;
75     dim_t n_block=A_p->row_block_size;     dim_t n_block=A_p->row_block_size;
# Line 77  Paso_Solver_LocalGS* Paso_Solver_getLoca Line 77  Paso_Solver_LocalGS* Paso_Solver_getLoca
77        
78     double time0=Paso_timer();     double time0=Paso_timer();
79     /* allocations: */       /* allocations: */  
80     Paso_Solver_LocalGS* out=MEMALLOC(1,Paso_Solver_LocalGS);     Paso_Preconditioner_LocalGS* out=MEMALLOC(1,Paso_Preconditioner_LocalGS);
81     if (! Paso_checkPtr(out)) {     if (! Paso_checkPtr(out)) {
82                
83        out->diag=MEMALLOC( ((size_t) n) * ((size_t) block_size),double);        out->diag=MEMALLOC( ((size_t) n) * ((size_t) block_size),double);
# Line 96  Paso_Solver_LocalGS* Paso_Solver_getLoca Line 96  Paso_Solver_LocalGS* Paso_Solver_getLoca
96        if (verbose) printf("timing: Gauss-Seidel preparation: elemination : %e\n",time0);        if (verbose) printf("timing: Gauss-Seidel preparation: elemination : %e\n",time0);
97        return out;        return out;
98     } else {     } else {
99        Paso_Solver_LocalGS_free(out);        Paso_Preconditioner_LocalGS_free(out);
100        return NULL;        return NULL;
101     }     }
102  }  }
# Line 114  if GS is local the defect b - A x_{k-1} Line 114  if GS is local the defect b - A x_{k-1}
114    
115  */  */
116    
117  void Paso_Solver_solveGS(Paso_SystemMatrix* A_p, Paso_Solver_GS * gs, double * x, const double * b)  void Paso_Preconditioner_GS_solve(Paso_SystemMatrix* A_p, Paso_Preconditioner_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) * (A_p->mainBlock->row_block_size);     const dim_t n= (A_p->mainBlock->numRows) * (A_p->mainBlock->row_block_size);
# Line 123  void Paso_Solver_solveGS(Paso_SystemMatr Line 123  void Paso_Solver_solveGS(Paso_SystemMatr
123     dim_t sweeps=gs->localGS->sweeps;     dim_t sweeps=gs->localGS->sweeps;
124        
125     if (gs->is_local) {     if (gs->is_local) {
126        Paso_Solver_solveLocalGS(A_p->mainBlock,gs->localGS,x,b);        Paso_Preconditioner_LocalGS_solve(A_p->mainBlock,gs->localGS,x,b);
127     } else {     } else {
128        #pragma omp parallel for private(i) schedule(static)        #pragma omp parallel for private(i) schedule(static)
129        for (i=0;i<n;++i) x[i]=b[i];        for (i=0;i<n;++i) x[i]=b[i];
130                
131        Paso_Solver_localGSSweep(A_p->mainBlock,gs->localGS,x);        Paso_Preconditioner_LocalGS_Sweep(A_p->mainBlock,gs->localGS,x);
132                
133        while (sweeps > 1 ) {        while (sweeps > 1 ) {
134       #pragma omp parallel for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
# Line 136  void Paso_Solver_solveGS(Paso_SystemMatr Line 136  void Paso_Solver_solveGS(Paso_SystemMatr
136    
137           Paso_SystemMatrix_MatrixVector((-1.), A_p, x, 1., b_new); /* b_new = b - A*x */           Paso_SystemMatrix_MatrixVector((-1.), A_p, x, 1., b_new); /* b_new = b - A*x */
138            
139       Paso_Solver_localGSSweep(A_p->mainBlock,gs->localGS,b_new);       Paso_Preconditioner_LocalGS_Sweep(A_p->mainBlock,gs->localGS,b_new);
140            
141       #pragma omp parallel for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
142       for (i=0;i<n;++i) x[i]+=b_new[i];       for (i=0;i<n;++i) x[i]+=b_new[i];
# Line 145  void Paso_Solver_solveGS(Paso_SystemMatr Line 145  void Paso_Solver_solveGS(Paso_SystemMatr
145                
146     }     }
147  }  }
148  void Paso_Solver_solveLocalGS(Paso_SparseMatrix* A_p, Paso_Solver_LocalGS * gs, double * x, const double * b)  void Paso_Preconditioner_LocalGS_solve(Paso_SparseMatrix* A_p, Paso_Preconditioner_LocalGS * gs, double * x, const double * b)
149  {  {
150     register dim_t i;     register dim_t i;
151     const dim_t n= (A_p->numRows) * (A_p->row_block_size);     const dim_t n= (A_p->numRows) * (A_p->row_block_size);
# Line 155  void Paso_Solver_solveLocalGS(Paso_Spars Line 155  void Paso_Solver_solveLocalGS(Paso_Spars
155     #pragma omp parallel for private(i) schedule(static)     #pragma omp parallel for private(i) schedule(static)
156     for (i=0;i<n;++i) x[i]=b[i];     for (i=0;i<n;++i) x[i]=b[i];
157        
158     Paso_Solver_localGSSweep(A_p,gs,x);     Paso_Preconditioner_LocalGS_Sweep(A_p,gs,x);
159        
160     while (sweeps > 1 ) {     while (sweeps > 1 ) {
161       #pragma omp parallel for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
# Line 163  void Paso_Solver_solveLocalGS(Paso_Spars Line 163  void Paso_Solver_solveLocalGS(Paso_Spars
163            
164       Paso_SparseMatrix_MatrixVector_CSC_OFFSET0((-1.), A_p, x, 1., b_new); /* b_new = b - A*x */       Paso_SparseMatrix_MatrixVector_CSC_OFFSET0((-1.), A_p, x, 1., b_new); /* b_new = b - A*x */
165            
166       Paso_Solver_localGSSweep(A_p,gs,b_new);       Paso_Preconditioner_LocalGS_Sweep(A_p,gs,b_new);
167            
168       #pragma omp parallel for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
169       for (i=0;i<n;++i) x[i]+=b_new[i];       for (i=0;i<n;++i) x[i]+=b_new[i];
# Line 172  void Paso_Solver_solveLocalGS(Paso_Spars Line 172  void Paso_Solver_solveLocalGS(Paso_Spars
172     }     }
173  }  }
174    
175  void Paso_Solver_localGSSweep(Paso_SparseMatrix* A, Paso_Solver_LocalGS * gs, double * x)  void Paso_Preconditioner_LocalGS_Sweep(Paso_SparseMatrix* A, Paso_Preconditioner_LocalGS * gs, double * x)
176  {  {
177     #ifdef _OPENMP     #ifdef _OPENMP
178     const dim_t nt=omp_get_max_threads();     const dim_t nt=omp_get_max_threads();
# Line 180  void Paso_Solver_localGSSweep(Paso_Spars Line 180  void Paso_Solver_localGSSweep(Paso_Spars
180     const dim_t nt = 1;     const dim_t nt = 1;
181     #endif     #endif
182     if (nt < 2) {     if (nt < 2) {
183        Paso_Solver_localGSSweep_sequential(A,gs,x);        Paso_Preconditioner_LocalGS_Sweep_sequential(A,gs,x);
184     } else {     } else {
185        Paso_Solver_localGSSweep_colored(A,gs,x);        Paso_Preconditioner_LocalGS_Sweep_colored(A,gs,x);
186     }     }
187  }  }
188    
189  /* inplace Gauss-Seidel sweep in seqential mode: */  /* inplace Gauss-Seidel sweep in seqential mode: */
190    
191  void Paso_Solver_localGSSweep_sequential(Paso_SparseMatrix* A_p, Paso_Solver_LocalGS * gs, double * x)  void Paso_Preconditioner_LocalGS_Sweep_sequential(Paso_SparseMatrix* A_p, Paso_Preconditioner_LocalGS * gs, double * x)
192  {  {
193     const dim_t n=A_p->numRows;     const dim_t n=A_p->numRows;
194     const dim_t n_block=A_p->row_block_size;     const dim_t n_block=A_p->row_block_size;
# Line 276  void Paso_Solver_localGSSweep_sequential Line 276  void Paso_Solver_localGSSweep_sequential
276     return;     return;
277  }  }
278                
279  void Paso_Solver_localGSSweep_colored(Paso_SparseMatrix* A_p, Paso_Solver_LocalGS * gs, double * x)  void Paso_Preconditioner_LocalGS_Sweep_colored(Paso_SparseMatrix* A_p, Paso_Preconditioner_LocalGS * gs, double * x)
280  {  {
281     const dim_t n=A_p->numRows;     const dim_t n=A_p->numRows;
282     const dim_t n_block=A_p->row_block_size;     const dim_t n_block=A_p->row_block_size;
# Line 424  void Paso_Solver_localGSSweep_colored(Pa Line 424  void Paso_Solver_localGSSweep_colored(Pa
424        }        }
425     }     }
426     return;     return;
 }  
   
   
   
427    }

Legend:
Removed from v.3120  
changed lines
  Added in v.3125

  ViewVC Help
Powered by ViewVC 1.1.26