/[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 3158 by gross, Mon Sep 6 06:09:11 2010 UTC revision 3159 by gross, Mon Sep 6 06:59:31 2010 UTC
# Line 121  if Smoother is local the defect b - A x_ Line 121  if Smoother is local the defect b - A x_
121    
122  */  */
123    
124  void Paso_Preconditioner_Smoother_solve(Paso_SystemMatrix* A_p, Paso_Preconditioner_Smoother * smoother, double * x, const double * b, const dim_t sweeps)  void Paso_Preconditioner_Smoother_solve(Paso_SystemMatrix* A_p, Paso_Preconditioner_Smoother * smoother, double * x, const double * b,
125                        const dim_t sweeps, const bool_t x_is_initial)
126  {  {
127     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);
128        
129     double *b_new = smoother->localSmoother->buffer;     double *b_new = smoother->localSmoother->buffer;
130     dim_t nsweeps=sweeps;     dim_t nsweeps=sweeps;
131     if (smoother->is_local) {     if (smoother->is_local) {
132        Paso_Preconditioner_LocalSmoother_solve(A_p->mainBlock,smoother->localSmoother,x,b,sweeps);        Paso_Preconditioner_LocalSmoother_solve(A_p->mainBlock,smoother->localSmoother,x,b,sweeps,x_is_initial);
133     } else {     } else {
134        Paso_Copy(n, x, b);        if (! x_is_initial) {
135               Paso_Copy(n, x, b);
136        Paso_Preconditioner_LocalSmoother_Sweep(A_p->mainBlock,smoother->localSmoother,x);       Paso_Preconditioner_LocalSmoother_Sweep(A_p->mainBlock,smoother->localSmoother,x);
137               nsweeps--;
138        while (nsweeps > 1 ) {        }
139          while (nsweeps > 0 ) {
140       Paso_Copy(n, b_new, b);       Paso_Copy(n, b_new, b);
   
141           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 */
142             Paso_Preconditioner_LocalSmoother_Sweep(A_p->mainBlock,smoother->localSmoother,b_new);  
      Paso_Preconditioner_LocalSmoother_Sweep(A_p->mainBlock,smoother->localSmoother,b_new);  
       
143       Paso_AXPY(n, x, 1., b_new);       Paso_AXPY(n, x, 1., b_new);
144       nsweeps--;       nsweeps--;
145        }        }
146                
147     }     }
148  }  }
149  void Paso_Preconditioner_LocalSmoother_solve(Paso_SparseMatrix* A_p, Paso_Preconditioner_LocalSmoother * smoother, double * x, const double * b, const dim_t sweeps)  void Paso_Preconditioner_LocalSmoother_solve(Paso_SparseMatrix* A_p, Paso_Preconditioner_LocalSmoother * smoother, double * x, const double * b,
150                             const dim_t sweeps, const bool_t x_is_initial)
151  {  {
152     const dim_t n= (A_p->numRows) * (A_p->row_block_size);     const dim_t n= (A_p->numRows) * (A_p->row_block_size);
153     double *b_new = smoother->buffer;     double *b_new = smoother->buffer;
154     dim_t nsweeps=sweeps;     dim_t nsweeps=sweeps;
155        
156     Paso_Copy(n, x, b);     if (! x_is_initial) {
157     Paso_Preconditioner_LocalSmoother_Sweep(A_p,smoother,x);        Paso_Copy(n, x, b);
158          Paso_Preconditioner_LocalSmoother_Sweep(A_p,smoother,x);
159          nsweeps--;
160       }
161        
162     while (nsweeps > 1 ) {     while (nsweeps > 0 ) {
   
163       Paso_Copy(n, b_new, b);       Paso_Copy(n, b_new, b);
       
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       Paso_Preconditioner_LocalSmoother_Sweep(A_p,smoother,b_new);       Paso_Preconditioner_LocalSmoother_Sweep(A_p,smoother,b_new);
166         Paso_AXPY(n, x, 1., b_new);
      Paso_AXPY(n, x, 1., b_new)  
       
167       nsweeps--;       nsweeps--;
168     }     }
169  }  }

Legend:
Removed from v.3158  
changed lines
  Added in v.3159

  ViewVC Help
Powered by ViewVC 1.1.26