/[escript]/trunk/paso/src/Preconditioner.c
ViewVC logotype

Diff of /trunk/paso/src/Preconditioner.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 2992 by artak, Tue Mar 23 04:16:28 2010 UTC revision 2993 by artak, Tue Mar 23 04:29:30 2010 UTC
# Line 200  void Paso_Solver_solvePreconditioner(Pas Line 200  void Paso_Solver_solvePreconditioner(Pas
200                            
201             Paso_Solver_solveGS(prec->gs,x,b);             Paso_Solver_solveGS(prec->gs,x,b);
202             if (prec->gs->sweeps>1) {             if (prec->gs->sweeps>1) {
                 double *bold=MEMALLOC(prec->gs->n*prec->gs->n_block,double);  
203                  double *bnew=MEMALLOC(prec->gs->n*prec->gs->n_block,double);                  double *bnew=MEMALLOC(prec->gs->n*prec->gs->n_block,double);
204                  dim_t sweeps=prec->gs->sweeps;                  dim_t sweeps=prec->gs->sweeps;
205                        
206                  #pragma omp parallel for private(i) schedule(static)                  #pragma omp parallel for private(i) schedule(static)
207                  for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bold[i]=b[i];                  for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bnew[i]=b[i];
208                        
209                  while(sweeps>1) {                  while(sweeps>1) {
210                     #pragma omp parallel for private(i) schedule(static)                     #pragma omp parallel for private(i) schedule(static)
211                     for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bnew[i]=bold[i]+b[i];                     for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bnew[i]+=b[i];
212                      /* Compute the residual b=b-Ax*/                      /* Compute the residual b=b-Ax*/
213                     Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), bnew);                     Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), bnew);
214                     /* Go round again*/                     /* Go round again*/
215                     Paso_Solver_solveGS(prec->gs,x,bnew);                     Paso_Solver_solveGS(prec->gs,x,bnew);
                    #pragma omp parallel for private(i) schedule(static)  
                    for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bold[i]=bnew[i];  
216                     sweeps=sweeps-1;                     sweeps=sweeps-1;
217                  }                  }
                 MEMFREE(bold);  
218                  MEMFREE(bnew);                  MEMFREE(bnew);
219             }             }
220             break;             break;

Legend:
Removed from v.2992  
changed lines
  Added in v.2993

  ViewVC Help
Powered by ViewVC 1.1.26