/[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 3282 by jfenwick, Mon Oct 11 01:48:14 2010 UTC revision 3283 by gross, Mon Oct 18 22:39:28 2010 UTC
# Line 37  void Paso_Preconditioner_free(Paso_Preco Line 37  void Paso_Preconditioner_free(Paso_Preco
37        Paso_Preconditioner_Smoother_free(in->jacobi);        Paso_Preconditioner_Smoother_free(in->jacobi);
38        Paso_Preconditioner_Smoother_free(in->gs);        Paso_Preconditioner_Smoother_free(in->gs);
39        Paso_Preconditioner_LocalAMG_free(in->localamg);        Paso_Preconditioner_LocalAMG_free(in->localamg);
40                Paso_Preconditioner_LocalSmoother_free(in->localamgsubstitute);
41        /*********************************/        /*********************************/
42        Paso_Solver_ILU_free(in->ilu);        Paso_Solver_ILU_free(in->ilu);
43        Paso_Solver_RILU_free(in->rilu);        Paso_Solver_RILU_free(in->rilu);
       Paso_Solver_AMLI_free(in->amli);  
       Paso_Solver_AMLI_System_free(in->amliSystem);  
44        /*********************************/        /*********************************/
45                
46        MEMFREE(in);        MEMFREE(in);
# Line 52  void Paso_Preconditioner_free(Paso_Preco Line 50  void Paso_Preconditioner_free(Paso_Preco
50  Paso_Preconditioner* Paso_Preconditioner_alloc(Paso_SystemMatrix* A,Paso_Options* options) {  Paso_Preconditioner* Paso_Preconditioner_alloc(Paso_SystemMatrix* A,Paso_Options* options) {
51    
52      Paso_Preconditioner* prec=NULL;      Paso_Preconditioner* prec=NULL;
     dim_t i;  
53    
54      prec=MEMALLOC(1,Paso_Preconditioner);      prec=MEMALLOC(1,Paso_Preconditioner);
55    
# Line 64  Paso_Preconditioner* Paso_Preconditioner Line 61  Paso_Preconditioner* Paso_Preconditioner
61      prec->jacobi=NULL;      prec->jacobi=NULL;
62      prec->gs=NULL;      prec->gs=NULL;
63      prec->localamg=NULL;      prec->localamg=NULL;
64        prec->localamgsubstitute=NULL;
65        
66      /*********************************/      /*********************************/
67          prec->rilu=NULL;          prec->rilu=NULL;
68          prec->ilu=NULL;          prec->ilu=NULL;
   
         prec->amli=NULL;  
         prec->amliSystem=NULL;  
69      /*********************************/      /*********************************/
70            
71      if (options->verbose && options->use_local_preconditioner) printf("Apply preconditioner locally only.\n");      if (options->verbose && options->use_local_preconditioner) printf("Paso: Apply preconditioner locally only.\n");
72    
73          switch (options->preconditioner) {          switch (options->preconditioner) {
74             default:             default:
75             case PASO_JACOBI:             case PASO_JACOBI:
76            if (options->verbose) printf("Jacobi(%d) preconditioner is used.\n",options->sweeps);            if (options->verbose) printf("Paso: Jacobi(%d) preconditioner is used.\n",options->sweeps);
77            prec->jacobi=Paso_Preconditioner_Smoother_alloc(A, TRUE, options->use_local_preconditioner, options->verbose);            prec->jacobi=Paso_Preconditioner_Smoother_alloc(A, TRUE, options->use_local_preconditioner, options->verbose);
78                prec->type=PASO_JACOBI;                prec->type=PASO_JACOBI;
79            prec->sweeps=options->sweeps;            prec->sweeps=options->sweeps;
80                break;                break;
81         case PASO_GS:         case PASO_GS:
82            if (options->verbose) printf("Gauss-Seidel(%d) preconditioner is used.\n",options->sweeps);            if (options->verbose) printf("Paso: Gauss-Seidel(%d) preconditioner is used.\n",options->sweeps);
83            prec->gs=Paso_Preconditioner_Smoother_alloc(A, FALSE, options->use_local_preconditioner, options->verbose);            prec->gs=Paso_Preconditioner_Smoother_alloc(A, FALSE, options->use_local_preconditioner, options->verbose);
84            prec->type=PASO_GS;            prec->type=PASO_GS;
85            prec->sweeps=options->sweeps;            prec->sweeps=options->sweeps;
86            break;            break;
87           case PASO_AMLI:
88         case PASO_AMG:         case PASO_AMG:
89            if (options->verbose) printf("AMG preconditioner is used.\n");            if (options->verbose) printf("Paso: AMG preconditioner is used.\n");
90            prec->localamg=Paso_Preconditioner_LocalAMG_alloc(A->mainBlock,options->level_max,options);            prec->localamg=Paso_Preconditioner_LocalAMG_alloc(A->mainBlock,1,options);
91            Esys_MPIInfo_noError(A->mpi_info);            prec->sweeps=options->sweeps;
92              /* if NULL is returned (and no error) no AMG has been constructed because the system is too small or not big enough
93                 we now use the Smoother as a preconditioner                                                                      */
94              if ( (Esys_noError()) && (prec->localamg == NULL) ) {
95             if (options->verbose) {
96                if (options->smoother == PASO_JACOBI) {
97                printf("Paso: Jacobi(%d) preconditioner is used.\n",prec->sweeps);
98                } else {
99                   printf("Paso: Gauss-Seidel(%d) preconditioner is used.\n",prec->sweeps);
100                }
101             }
102             prec->localamgsubstitute=Paso_Preconditioner_LocalSmoother_alloc(A->mainBlock, (options->smoother == PASO_JACOBI), options->verbose);
103              } else {
104             if (options->verbose) printf("Paso: AMG preconditioner is used.\n");
105              }
106            prec->type=PASO_AMG;            prec->type=PASO_AMG;
107              Esys_MPIInfo_noError(A->mpi_info);
108            break;            break;
109                        
110         /***************************************************************************************/           /***************************************************************************************/  
111             case PASO_ILU0:             case PASO_ILU0:
112                if (options->verbose) printf("ILU preconditioner is used.\n");            if (options->verbose) printf("Paso: ILU preconditioner is used.\n");
113                prec->ilu=Paso_Solver_getILU(A->mainBlock,options->verbose);                prec->ilu=Paso_Solver_getILU(A->mainBlock,options->verbose);
114                prec->type=PASO_ILU0;                prec->type=PASO_ILU0;
115            Esys_MPIInfo_noError(A->mpi_info);            Esys_MPIInfo_noError(A->mpi_info);
116                break;                break;
117             case PASO_RILU:             case PASO_RILU:
118                if (options->verbose) printf("RILU preconditioner is used.\n");            if (options->verbose) printf("Paso: RILU preconditioner is used.\n");
119                prec->rilu=Paso_Solver_getRILU(A->mainBlock,options->verbose);                prec->rilu=Paso_Solver_getRILU(A->mainBlock,options->verbose);
120            Esys_MPIInfo_noError(A->mpi_info);            Esys_MPIInfo_noError(A->mpi_info);
121                prec->type=PASO_RILU;                prec->type=PASO_RILU;
122                break;                break;
   
   
   
             case PASO_AMLI:  
             
           prec->amliSystem=MEMALLOC(1,Paso_Solver_AMLI_System);  
               if (! Esys_checkPtr(prec->amliSystem)) {  
                   
               prec->amliSystem->block_size=A->row_block_size;  
                 
           for (i=0;i<A->row_block_size;++i) {  
           prec->amliSystem->amliblock[i]=NULL;  
           prec->amliSystem->block[i]=NULL;  
           }  
   
           if (options->verbose) printf("AMLI preconditioner is used.\n");  
                 
               /*For performace reasons we check if block_size is one. If yes, then we do not need to separate blocks.*/  
               if (A->row_block_size==1) {  
                 prec->amli=Paso_Solver_getAMLI(A->mainBlock,options->level_max,options);    
               }  
               else {  
                 for (i=0;i<A->row_block_size;++i) {  
                 prec->amliSystem->block[i]=Paso_SparseMatrix_getBlock(A->mainBlock,i+1);  
                 prec->amliSystem->amliblock[i]=Paso_Solver_getAMLI(prec->amliSystem->block[i],options->level_max,options);  
                 }  
               }  
               prec->type=PASO_AMLI;  
           }  
               break;  
   
123          }          }
124      }      }
125      if (! Esys_MPIInfo_noError(A->mpi_info ) ){      if (! Esys_MPIInfo_noError(A->mpi_info ) ){
# Line 153  Paso_Preconditioner* Paso_Preconditioner Line 134  Paso_Preconditioner* Paso_Preconditioner
134  /* has to be called within a parallel reqion */  /* has to be called within a parallel reqion */
135  /* 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 */
136  void Paso_Preconditioner_solve(Paso_Preconditioner* prec, Paso_SystemMatrix* A,double* x,double* b){  void Paso_Preconditioner_solve(Paso_Preconditioner* prec, Paso_SystemMatrix* A,double* x,double* b){
137      dim_t i,j;  
     dim_t n=A->mainBlock->numRows;  
138      switch (prec->type) {      switch (prec->type) {
139          default:          default:
140          case PASO_JACOBI:          case PASO_JACOBI:
# Line 164  void Paso_Preconditioner_solve(Paso_Prec Line 144  void Paso_Preconditioner_solve(Paso_Prec
144         Paso_Preconditioner_Smoother_solve(A, prec->gs,x,b,prec->sweeps, FALSE);         Paso_Preconditioner_Smoother_solve(A, prec->gs,x,b,prec->sweeps, FALSE);
145         break;         break;
146      case PASO_AMG:      case PASO_AMG:
147         Paso_Preconditioner_LocalAMG_solve(prec->localamg,x,b);         if (prec->localamg == NULL) {
148              Paso_Preconditioner_LocalSmoother_solve(A->mainBlock, prec->localamgsubstitute,x,b,prec->sweeps, FALSE);
149           } else {
150              Paso_Preconditioner_LocalAMG_solve(A->mainBlock, prec->localamg,x,b);
151           }
152         break;         break;
153                
154      /*=========================================================*/        /*=========================================================*/  
# Line 175  void Paso_Preconditioner_solve(Paso_Prec Line 159  void Paso_Preconditioner_solve(Paso_Prec
159         Paso_Solver_solveRILU(prec->rilu,x,b);         Paso_Solver_solveRILU(prec->rilu,x,b);
160             break;             break;
161    
   
         case PASO_AMLI:  
               
             /*For performace reasons we check if block_size is one. If yes, then we do not need to do unnecessary copying.*/  
             if (A->row_block_size==1) {  
                 Paso_Solver_solveAMLI(prec->amli,x,b);  
             }  
             else {  
            double **xx;  
            double **bb;  
            xx=MEMALLOC(A->row_block_size,double*);  
            bb=MEMALLOC(A->row_block_size,double*);  
            if (Esys_checkPtr(xx) || Esys_checkPtr(bb)) return;  
                  for (i=0;i<A->row_block_size;i++) {  
                     xx[i]=MEMALLOC(n,double);  
                     bb[i]=MEMALLOC(n,double);  
                     if (Esys_checkPtr(xx[i]) && Esys_checkPtr(bb[i])) return;  
                 }  
                   
                 /*#pragma omp parallel for private(i,j) schedule(static)*/  
                 for (i=0;i<n;i++) {  
                     for (j=0;j<A->row_block_size;j++) {  
                      bb[j][i]=b[A->row_block_size*i+j];  
                     }  
                  }  
                   
                   
                 for (i=0;i<A->row_block_size;i++) {  
                 Paso_Solver_solveAMLI(prec->amliSystem->amliblock[i],xx[i],bb[i]);  
                 }  
                                 
                 /*#pragma omp parallel for private(i,j) schedule(static)*/  
                 for (i=0;i<n;i++) {  
                     for (j=0;j<A->row_block_size;j++) {  
                     x[A->row_block_size*i+j]=xx[j][i];  
                     }  
                  }  
                   
                 for (i=0;i<A->row_block_size;i++) {  
                 MEMFREE(xx[i]);  
                 MEMFREE(bb[i]);  
                 }  
                 MEMFREE(xx);  
         MEMFREE(bb);  
             }  
         break;  
     /*=========================================================*/  
162      }      }
163    
164  }  }
Legend:
Removed from v.3282  
changed lines
  Added in v.3283
  ViewVC Help
Powered by ViewVC 1.1.26