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

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

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

revision 2661 by artak, Fri Sep 11 00:59:59 2009 UTC revision 2662 by artak, Tue Sep 15 03:05:23 2009 UTC
# Line 47  void Paso_Preconditioner_free(Paso_Solve Line 47  void Paso_Preconditioner_free(Paso_Solve
47  void Paso_Solver_setPreconditioner(Paso_SystemMatrix* A,Paso_Options* options) {  void Paso_Solver_setPreconditioner(Paso_SystemMatrix* A,Paso_Options* options) {
48    
49      Paso_Solver_Preconditioner* prec=NULL;      Paso_Solver_Preconditioner* prec=NULL;
50        dim_t i;
51      if (A->solver==NULL) {      if (A->solver==NULL) {
52          /* allocate structure to hold preconditioner */          /* allocate structure to hold preconditioner */
53          prec=MEMALLOC(1,Paso_Solver_Preconditioner);          prec=MEMALLOC(1,Paso_Solver_Preconditioner);
# Line 58  void Paso_Solver_setPreconditioner(Paso_ Line 59  void Paso_Solver_setPreconditioner(Paso_
59          prec->jacobi=NULL;          prec->jacobi=NULL;
60          prec->gs=NULL;          prec->gs=NULL;
61          prec->amg=NULL;          prec->amg=NULL;
62          prec->amgSystem->amgblock1=NULL;          
63          prec->amgSystem->amgblock2=NULL;          prec->amgSystem->block_size=A->row_block_size;
64          prec->amgSystem->amgblock3=NULL;          for (i=0;i<A->row_block_size;++i) {
65          prec->amgSystem->block1=NULL;            prec->amgSystem->amgblock[i]=NULL;
66          prec->amgSystem->block2=NULL;            prec->amgSystem->block[i]=NULL;
67          prec->amgSystem->block3=NULL;          }
68            
69    
70          A->solver=prec;          A->solver=prec;
71          switch (options->preconditioner) {          switch (options->preconditioner) {
72             default:             default:
# Line 91  void Paso_Solver_setPreconditioner(Paso_ Line 94  void Paso_Solver_setPreconditioner(Paso_
94              case PASO_AMG:              case PASO_AMG:
95                if (options->verbose) printf("AMG preconditioner is used.\n");                if (options->verbose) printf("AMG preconditioner is used.\n");
96                                
97                  /*For performace reasons we check if block_size is one. If yes, then we do not need to separate blocks.*/
98                if (A->row_block_size==1) {                if (A->row_block_size==1) {
99                  prec->amg=Paso_Solver_getAMG(A->mainBlock,options->level_max,options);                    prec->amg=Paso_Solver_getAMG(A->mainBlock,options->level_max,options);  
100                }                }
101                else if (A->row_block_size==2) {                else {
102                                    #pragma omp parallel for private(i) schedule(static)
103                  prec->amgSystem->block1=Paso_SparseMatrix_getBlock(A->mainBlock,1);                  for (i=0;i<A->row_block_size;++i) {
104                  prec->amgSystem->block2=Paso_SparseMatrix_getBlock(A->mainBlock,2);                  prec->amgSystem->block[i]=Paso_SparseMatrix_getBlock(A->mainBlock,i+1);
105                    prec->amgSystem->amgblock[i]=Paso_Solver_getAMG(prec->amgSystem->block[i],options->level_max,options);
106                  prec->amgSystem->amgblock1=Paso_Solver_getAMG(prec->amgSystem->block1,options->level_max,options);                  }
                 prec->amgSystem->amgblock2=Paso_Solver_getAMG(prec->amgSystem->block2,options->level_max,options);  
                   
107                }                }
               else if (A->row_block_size==3)  
               {  
108    
                 prec->amgSystem->block1=Paso_SparseMatrix_getBlock(A->mainBlock,1);  
                 prec->amgSystem->block2=Paso_SparseMatrix_getBlock(A->mainBlock,2);  
                 prec->amgSystem->block3=Paso_SparseMatrix_getBlock(A->mainBlock,3);  
                   
                 /*  
                 Paso_SparseMatrix_saveMM(prec->amgSystem->block1,"amgtemp1.mat");  
                 Paso_SparseMatrix_saveMM(prec->amgSystem->block2,"amgtemp2.mat");  
                 Paso_SparseMatrix_saveMM(prec->amgSystem->block3,"amgtemp3.mat");  
                 */  
                 prec->amgSystem->amgblock1=Paso_Solver_getAMG(prec->amgSystem->block1,options->level_max,options);  
                 prec->amgSystem->amgblock2=Paso_Solver_getAMG(prec->amgSystem->block2,options->level_max,options);  
                 prec->amgSystem->amgblock3=Paso_Solver_getAMG(prec->amgSystem->block3,options->level_max,options);  
               }  
109                prec->type=PASO_AMG;                prec->type=PASO_AMG;
110                break;                break;
111    
# Line 135  void Paso_Solver_setPreconditioner(Paso_ Line 122  void Paso_Solver_setPreconditioner(Paso_
122  /* barrier synchronization is performed before the evaluation to make sure that the input vector is available */  /* barrier synchronization is performed before the evaluation to make sure that the input vector is available */
123  void Paso_Solver_solvePreconditioner(Paso_SystemMatrix* A,double* x,double* b){  void Paso_Solver_solvePreconditioner(Paso_SystemMatrix* A,double* x,double* b){
124      Paso_Solver_Preconditioner* prec=(Paso_Solver_Preconditioner*) A->solver;      Paso_Solver_Preconditioner* prec=(Paso_Solver_Preconditioner*) A->solver;
125            dim_t i,j;
126        dim_t n=Paso_SystemMatrix_getGlobalNumRows(A);
127        double **xx;
128        double **bb;
129    
130      switch (prec->type) {      switch (prec->type) {
131          default:          default:
132          case PASO_JACOBI:          case PASO_JACOBI:
# Line 147  void Paso_Solver_solvePreconditioner(Pas Line 138  void Paso_Solver_solvePreconditioner(Pas
138          case PASO_RILU:          case PASO_RILU:
139             Paso_Solver_solveRILU(prec->rilu,x,b);             Paso_Solver_solveRILU(prec->rilu,x,b);
140             break;             break;
141           case PASO_GS:          case PASO_GS:
142              /* Gauss-Seidel preconditioner P=U^{-1}DL^{-1} is used here with sweeps paramenter.              /* Gauss-Seidel preconditioner P=U^{-1}DL^{-1} is used here with sweeps paramenter.
143                We want to solve x_new=x_old+P^{-1}(b-Ax_old). So for fisrt 3 we will have the following:                We want to solve x_new=x_old+P^{-1}(b-Ax_old). So for fisrt 3 we will have the following:
144                x_0=0                x_0=0
# Line 169  void Paso_Solver_solvePreconditioner(Pas Line 160  void Paso_Solver_solvePreconditioner(Pas
160    
161             Paso_Solver_solveGS(prec->gs,x,b);             Paso_Solver_solveGS(prec->gs,x,b);
162             if (prec->gs->sweeps>1) {             if (prec->gs->sweeps>1) {
163             double *bold=MEMALLOC(prec->gs->n*prec->gs->n_block,double);                  double *bold=MEMALLOC(prec->gs->n*prec->gs->n_block,double);
164             double *bnew=MEMALLOC(prec->gs->n*prec->gs->n_block,double);                  double *bnew=MEMALLOC(prec->gs->n*prec->gs->n_block,double);
            dim_t i;  
            #pragma omp parallel for private(i) schedule(static)  
            for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bold[i]=b[i];  
             
            while(prec->gs->sweeps>1) {  
                #pragma omp parallel for private(i) schedule(static)  
                for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bnew[i]=bold[i]+b[i];  
                 /* Compute the residual b=b-Ax*/  
                Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), bnew);  
                /* Go round again*/  
                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];  
                prec->gs->sweeps=prec->gs->sweeps-1;  
            }  
165                        
166             MEMFREE(bold);                  #pragma omp parallel for private(i) schedule(static)
167             MEMFREE(bnew);                  for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bold[i]=b[i];
168                        
169                    while(prec->gs->sweeps>1) {
170                       #pragma omp parallel for private(i) schedule(static)
171                       for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bnew[i]=bold[i]+b[i];
172                        /* Compute the residual b=b-Ax*/
173                       Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), bnew);
174                       /* Go round again*/
175                       Paso_Solver_solveGS(prec->gs,x,bnew);
176                       #pragma omp parallel for private(i) schedule(static)
177                       for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bold[i]=bnew[i];
178                       prec->gs->sweeps=prec->gs->sweeps-1;
179                    }
180                    MEMFREE(bold);
181                    MEMFREE(bnew);
182             }             }
183             break;             break;
184           case PASO_AMG:          case PASO_AMG:
185                
186                /*For performace reasons we check if block_size is one. If yes, then we do not need to do unnecessary copying.*/
187              if (A->row_block_size==1) {              if (A->row_block_size==1) {
188                  Paso_Solver_solveAMG(prec->amg,x,b);                  Paso_Solver_solveAMG(prec->amg,x,b);
189              }              }
190              else if (A->row_block_size==2) {              else {
191                  dim_t i;            
192                  double *x1=MEMALLOC(prec->amgSystem->amgblock1->n,double);                  xx=MEMALLOC(A->row_block_size,double*);
193                  double *x2=MEMALLOC(prec->amgSystem->amgblock2->n,double);                  bb=MEMALLOC(A->row_block_size,double*);
                 double *b1=MEMALLOC(prec->amgSystem->amgblock1->n,double);  
                 double *b2=MEMALLOC(prec->amgSystem->amgblock2->n,double);  
                   
                 #pragma omp parallel for private(i) schedule(static)  
                 for (i=0;i<prec->amgSystem->amgblock1->n;i++) {  
                     b1[i]=b[2*i];  
                     b2[i]=b[2*i+1];  
                 }  
                   
                   
                 Paso_Solver_solveAMG(prec->amgSystem->amgblock1,x1,b1);  
                 Paso_Solver_solveAMG(prec->amgSystem->amgblock2,x2,b2);  
194                                    
195                  #pragma omp parallel for private(i) schedule(static)                  for (i=0;i<A->row_block_size;i++) {
196                  for (i=0;i<prec->amgSystem->amgblock1->n;i++) {                      xx[i]=MEMALLOC(n,double);
197                      x[2*i]=x1[i];                      bb[i]=MEMALLOC(n,double);
                     x[2*i+1]=x2[i];  
198                  }                  }
                   
                 MEMFREE(x1);  
                 MEMFREE(x2);  
                 MEMFREE(b1);  
                 MEMFREE(b2);  
199    
200              }                  #pragma omp parallel for private(i,j) schedule(static)
201              else if (A->row_block_size==3) {                  for (i=0;i<n;i++) {
202                  dim_t i;                      for (j=0;j<A->row_block_size;j++) {
203                                         bb[j][i]=b[A->row_block_size*i+j];
204                  double *x1=MEMALLOC(prec->amgSystem->amgblock1->n,double);                      }
                 double *x2=MEMALLOC(prec->amgSystem->amgblock2->n,double);  
                 double *x3=MEMALLOC(prec->amgSystem->amgblock3->n,double);  
                 double *b1=MEMALLOC(prec->amgSystem->amgblock1->n,double);  
                 double *b2=MEMALLOC(prec->amgSystem->amgblock2->n,double);  
                 double *b3=MEMALLOC(prec->amgSystem->amgblock3->n,double);  
                   
                 #pragma omp parallel for private(i) schedule(static)  
                 for (i=0;i<prec->amgSystem->amgblock1->n;i++) {  
                     b1[i]=b[3*i];  
                     b2[i]=b[3*i+1];  
                     b3[i]=b[3*i+2];  
205                   }                   }
206                                                    
                 Paso_Solver_solveAMG(prec->amgSystem->amgblock1,x1,b1);  
                 Paso_Solver_solveAMG(prec->amgSystem->amgblock2,x2,b2);  
                 Paso_Solver_solveAMG(prec->amgSystem->amgblock3,x3,b3);  
                   
207                  #pragma omp parallel for private(i) schedule(static)                  #pragma omp parallel for private(i) schedule(static)
208                  for (i=0;i<prec->amgSystem->amgblock1->n;i++) {                  for (i=0;i<A->row_block_size;i++) {
209                      x[3*i]=x1[i];                  Paso_Solver_solveAMG(prec->amgSystem->amgblock[i],xx[i],bb[i]);
210                      x[3*i+1]=x2[i];                  }
211                      x[3*i+2]=x3[i];                    
212                    #pragma omp parallel for private(i,j) schedule(static)
213                    for (i=0;i<n;i++) {
214                        for (j=0;j<A->row_block_size;j++) {
215                        x[A->row_block_size*i+j]=xx[j][i];
216                        }
217                     }
218                                  
219                    for (i=0;i<A->row_block_size;i++) {
220                    MEMFREE(xx[i]);
221                    MEMFREE(bb[i]);
222                  }                  }
223                                    
224                                    MEMFREE(xx);
225                  MEMFREE(x1);                  MEMFREE(bb);
                 MEMFREE(x2);  
                 MEMFREE(x3);  
                 MEMFREE(b1);  
                 MEMFREE(b2);  
                 MEMFREE(b3);  
226              }              }
227                
228          break;          break;
229      }      }
230    

Legend:
Removed from v.2661  
changed lines
  Added in v.2662

  ViewVC Help
Powered by ViewVC 1.1.26