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

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

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

revision 3641 by gross, Mon Aug 30 10:48:00 2010 UTC revision 3642 by caltinay, Thu Oct 27 03:41:51 2011 UTC
# Line 14  Line 14 
14    
15  /**************************************************************/  /**************************************************************/
16    
17  /* Paso: GS preconditioner with reordering                 */  /* Paso: GS preconditioner with reordering                    */
18    
19  /**************************************************************/  /**************************************************************/
20    
21  /* Copyrights by ACcESS Australia 2003-2010 */  /* Copyrights by ACcESS Australia 2003-2010 */
22  /* Author: l.gao@uq.edu.au                                   */  /* Author: l.gao@uq.edu.au                                    */
23    
24  /**************************************************************/  /**************************************************************/
25    
# Line 45  void Paso_Solver_GSMPI_free(Paso_Solver_ Line 45  void Paso_Solver_GSMPI_free(Paso_Solver_
45       }       }
46  }  }
47    
48  =============================================================================================  ==========================================================================
49  /**************************************************************/  /**************************************************************/
50    
51  /*   gs->diag saves the matrix of D{-1}  /*   gs->diag saves the matrix of D{-1}
52       this is different from Paso_Solver_getGS(), in which, gs->diag       This is different from Paso_Solver_getGS(), in which, gs->diag
53       is the matrix D.       is the matrix D.
54  */  */
55  Paso_Solver_GS* Paso_Solver_getGSMPI(Paso_SparseMatrix * A,bool_t verbose) {  Paso_Solver_GS* Paso_Solver_getGSMPI(Paso_SparseMatrix * A,bool_t verbose) {
# Line 200  void Paso_Solver_GS_local(Paso_SystemMat Line 200  void Paso_Solver_GS_local(Paso_SystemMat
200    
201  /**************************************************************/  /**************************************************************/
202    
203  /* apply MPI versioned GS                                  /* Applies MPI versioned GS
204    
205       in fact it solves Ax=b in two steps:       In fact it solves Ax=b in two steps:
206       step1: among different nodes (MPI ranks), we use block Jacobi       step 1: among different nodes (MPI ranks), we use block Jacobi
207         x{k} = x{k-1} + D{-1}(b-A*x{k-1})         x{k} = x{k-1} + D{-1}(b-A*x{k-1})
208          => D*x{k} = b - (E+F)x{k-1}          => D*x{k} = b - (E+F)x{k-1}
209            where matrix D is (let p be the number of nodes):            where matrix D is (let p be the number of nodes):
# Line 221  void Paso_Solver_GS_local(Paso_SystemMat Line 221  void Paso_Solver_GS_local(Paso_SystemMat
221            and Ai (i \in [1,p]) represents the mainBlock of matrix            and Ai (i \in [1,p]) represents the mainBlock of matrix
222            A on node i. Matrix (E+F) is represented as the coupleBlock            A on node i. Matrix (E+F) is represented as the coupleBlock
223                of matrix A on each node (annotated as ACi).                of matrix A on each node (annotated as ACi).
224             Therefore, step1 can be turned into the following for node i:             Therefore, step 1 can be turned into the following for node i:
225         => Ai * x{k} = b - ACi * x{k-1}         => Ai * x{k} = b - ACi * x{k-1}
226             where both x{k} and b are the segment of x and b on node i,             where both x{k} and b are the segment of x and b on node i,
227         and x{k-1} is the old segment values of x on all other nodes.         and x{k-1} is the old segment values of x on all other nodes.
228    
229       step2: inside node i, we use Gauss-Seidel       step 2: inside node i, we use Gauss-Seidel
230           let b'= b - ACi * x{k-1} we have Ai * x{k} = b' for node i           let b'= b - ACi * x{k-1} we have Ai * x{k} = b' for node i
231           by using symetrix Gauss-Seidel, this can be solved in a forward           by using symmetric Gauss-Seidel, this can be solved in a forward
232           phase and a backward phase:           phase and a backward phase:
233         forward phase:  x{m} = diag(Ai){-1} (b' - E*x{m} - F*x{m-1})         forward phase:  x{m} = diag(Ai){-1} (b' - E*x{m} - F*x{m-1})
234             backward phase: x{m+1} = diag(Ai){-1} (b' - F*{m+1} - E*x{m})             backward phase: x{m+1} = diag(Ai){-1} (b' - F*{m+1} - E*x{m})
# Line 256  void Paso_Solver_solveGSMPI(Paso_SystemM Line 256  void Paso_Solver_solveGSMPI(Paso_SystemM
256    
257            while (sweeps > 1) {            while (sweeps > 1) {
258                 /* calculate new_b = b - ACi * x{k-1}, where x{k-1} are remote                 /* calculate new_b = b - ACi * x{k-1}, where x{k-1} are remote
259                    value of x, which requires MPI communications */                    values of x, which requires MPI communication */
260                 #pragma omp parallel for private(i) schedule(static)                 #pragma omp parallel for private(i) schedule(static)
261                 for (i=0;i<n*n_block;++i) new_b[i]=b[i];                 for (i=0;i<n*n_block;++i) new_b[i]=b[i];
262    
# Line 292  void Paso_Solver_GS_local(Paso_SystemMat Line 292  void Paso_Solver_GS_local(Paso_SystemMat
292       len=n/nt;       len=n/nt;
293       rest=n-len*nt;       rest=n-len*nt;
294  #endif  #endif
295       /* TO BE DONE: add handler to deal with the case "n" is too small       /* TO BE DONE: add handler to deal with the case "n is too small"
296                      to be worth run in threads. */                      to be worth run in threads. */
297    
298  #ifdef _OPENMP  #ifdef _OPENMP
299       /* calculate new_b = b - ACi * x{k-1}, where x{k-1} are x value       /* calculate new_b = b - ACi * x{k-1}, where x{k-1} are x values
300          computed by other threads in previous sweep */          computed by other threads in previous sweep */
301       if (nt > 1) {       if (nt > 1) {
302       if (n_block == 1){       if (n_block == 1){
# Line 371  void Paso_Solver_GS_local(Paso_SystemMat Line 371  void Paso_Solver_GS_local(Paso_SystemMat
371       }       }
372  #endif  #endif
373    
374       /* step1: forward iteration       /* step 1: forward iteration
375                 x{k} = D{-1}(b - E*x{k} - F*x{k-1}) */                 x{k} = D{-1}(b - E*x{k} - F*x{k-1}) */
376       /* One Guass-Seidel itertion       /* One Guass-Seidel iteration
377          In case of forward iteration x{k} = D{-1}(b - E*x{k} - F*x{k-1})          In case of forward iteration x{k} = D{-1}(b - E*x{k} - F*x{k-1})
378             => into a loop (without coloring):             => into a loop (without coloring):
379              for i in [0,n-1] do                  for i in [0,n-1] do    
# Line 491  void Paso_Solver_GS_local(Paso_SystemMat Line 491  void Paso_Solver_GS_local(Paso_SystemMat
491  #endif  #endif
492       }       }
493    
494       /* step2: backward iteration       /* step 2: backward iteration
495                 x{k} = D{-1}(b - F*x{k} - E*x{k-1}) */                 x{k} = D{-1}(b - F*x{k} - E*x{k-1}) */
496       if (n_block == 1){       if (n_block == 1){
497  #ifdef _OPENMP  #ifdef _OPENMP
# Line 600  void Paso_Solver_GS_local(Paso_SystemMat Line 600  void Paso_Solver_GS_local(Paso_SystemMat
600           }           }
601  #endif  #endif
602       }       }
   
      return;  
603  }  }
604    

Legend:
Removed from v.3641  
changed lines
  Added in v.3642

  ViewVC Help
Powered by ViewVC 1.1.26