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

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

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

revision 2806 by artak, Fri Dec 4 05:21:31 2009 UTC revision 2807 by artak, Mon Dec 7 00:02:55 2009 UTC
# Line 159  Paso_Solver_AMLI* Paso_Solver_getAMLI(Pa Line 159  Paso_Solver_AMLI* Paso_Solver_getAMLI(Pa
159       out->n=n;       out->n=n;
160       out->n_F=n+1;       out->n_F=n+1;
161       out->n_block=n_block;       out->n_block=n_block;
162         out->post_sweeps=options->post_sweeps;
163         out->pre_sweeps=options->pre_sweeps;
164            
165       if (level==0 || n<=options->min_coarse_matrix_size) {       if (level==0 || n<=options->min_coarse_matrix_size) {
166           out->coarsest_level=TRUE;           out->coarsest_level=TRUE;
# Line 406  Paso_Solver_AMLI* Paso_Solver_getAMLI(Pa Line 408  Paso_Solver_AMLI* Paso_Solver_getAMLI(Pa
408  void Paso_Solver_solveAMLI(Paso_Solver_AMLI * amli, double * x, double * b) {  void Paso_Solver_solveAMLI(Paso_Solver_AMLI * amli, double * x, double * b) {
409       dim_t i;       dim_t i;
410       double time0=0;       double time0=0;
411       double *r=NULL, *x0=NULL;       double *r=NULL, *x0=NULL,*x_F_temp=NULL;
412       bool_t verbose=0;       bool_t verbose=0;
413        
414         dim_t post_sweeps=amli->post_sweeps;
415         dim_t pre_sweeps=amli->pre_sweeps;
416        
417       #ifdef UMFPACK       #ifdef UMFPACK
418            Paso_UMFPACK_Handler * ptr=NULL;            Paso_UMFPACK_Handler * ptr=NULL;
419       #endif       #endif
420            
421       r=MEMALLOC(amli->n,double);       r=MEMALLOC(amli->n,double);
422       x0=MEMALLOC(amli->n,double);       x0=MEMALLOC(amli->n,double);
423         x_F_temp=MEMALLOC(amli->n_F,double);
424            
425       if (amli->coarsest_level) {       if (amli->coarsest_level) {
426                
# Line 442  void Paso_Solver_solveAMLI(Paso_Solver_A Line 449  void Paso_Solver_solveAMLI(Paso_Solver_A
449          /* presmoothing */          /* presmoothing */
450           time0=Paso_timer();           time0=Paso_timer();
451           Paso_Solver_solveJacobi(amli->GS,x,b);           Paso_Solver_solveJacobi(amli->GS,x,b);
452            
453             /***************/
454             #pragma omp parallel for private(i) schedule(static)
455             for (i=0;i<amli->n;++i) r[i]=b[i];
456      
457             while(pre_sweeps>1) {
458                 #pragma omp parallel for private(i) schedule(static)
459                 for (i=0;i<amli->n;++i) r[i]+=b[i];
460                
461                  /* Compute the residual b=b-Ax*/
462                 Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amli->A,x,1.,r);
463                 /* Go round again*/
464                 Paso_Solver_solveJacobi(amli->GS,x,r);
465                 pre_sweeps-=1;
466             }
467             /***************/  
468            
469           time0=Paso_timer()-time0;           time0=Paso_timer()-time0;
470           if (verbose) fprintf(stderr,"timing: Presmooting: %e\n",time0);           if (verbose) fprintf(stderr,"timing: Presmooting: %e\n",time0);
471          /* end of presmoothing */          /* end of presmoothing */
# Line 475  void Paso_Solver_solveAMLI(Paso_Solver_A Line 499  void Paso_Solver_solveAMLI(Paso_Solver_A
499                    
500          time0=Paso_timer();          time0=Paso_timer();
501                    
502          /* b_F=b_F-A_FC*x_C */          /* b_F=-A_FC*x_C */
503          Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amli->A_FC,amli->x_C,1.,amli->b_F);          Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amli->A_FC,amli->x_C,0.,amli->b_F);
504          /* x_F=invA_FF*b_F  */          /* x_F_temp=invA_FF*b_F  */
505          Paso_Solver_applyBlockDiagonalMatrix(1,amli->n_F,amli->inv_A_FF,amli->A_FF_pivot,amli->x_F,amli->b_F);          Paso_Solver_applyBlockDiagonalMatrix(1,amli->n_F,amli->inv_A_FF,amli->A_FF_pivot,x_F_temp,amli->b_F);
506            
507            #pragma omp parallel for private(i) schedule(static)
508            for (i=0;i<amli->n_F;++i) {
509                     amli->x_F[i]+=x_F_temp[i];
510            }
511                    
512          /* x<-[x_F,x_C]     */          /* x<-[x_F,x_C]     */
513          #pragma omp parallel for private(i) schedule(static)          #pragma omp parallel for private(i) schedule(static)
# Line 505  void Paso_Solver_solveAMLI(Paso_Solver_A Line 534  void Paso_Solver_solveAMLI(Paso_Solver_A
534       #pragma omp parallel for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
535       for (i=0;i<amli->n;++i) x[i]+=x0[i];       for (i=0;i<amli->n;++i) x[i]+=x0[i];
536            
537        
538         /***************/
539           while(post_sweeps>1) {
540              
541              #pragma omp parallel for private(i) schedule(static)
542              for (i=0;i<amli->n;++i) r[i]=b[i];
543              
544              Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amli->A,x,1.,r);
545              Paso_Solver_solveJacobi(amli->GS,x0,r);
546              #pragma omp parallel for private(i) schedule(static)
547               for (i=0;i<amli->n;++i)  {
548                x[i]+=x0[i];
549               }
550              post_sweeps-=1;
551           }
552           /**************/
553        
554       time0=Paso_timer()-time0;       time0=Paso_timer()-time0;
555       if (verbose) fprintf(stderr,"timing: Postsmoothing: %e\n",time0);       if (verbose) fprintf(stderr,"timing: Postsmoothing: %e\n",time0);
556    
# Line 513  void Paso_Solver_solveAMLI(Paso_Solver_A Line 559  void Paso_Solver_solveAMLI(Paso_Solver_A
559       }       }
560       MEMFREE(r);       MEMFREE(r);
561       MEMFREE(x0);       MEMFREE(x0);
562         MEMFREE(x_F_temp);
563       return;       return;
564  }  }

Legend:
Removed from v.2806  
changed lines
  Added in v.2807

  ViewVC Help
Powered by ViewVC 1.1.26