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

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

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

revision 1938 by artak, Mon Oct 27 01:31:28 2008 UTC revision 1939 by artak, Tue Oct 28 03:47:08 2008 UTC
# Line 35  Line 35 
35  void Paso_Solver_AMG_free(Paso_Solver_AMG * in) {  void Paso_Solver_AMG_free(Paso_Solver_AMG * in) {
36       if (in!=NULL) {       if (in!=NULL) {
37          Paso_Solver_AMG_free(in->AMG_of_Schur);          Paso_Solver_AMG_free(in->AMG_of_Schur);
38            Paso_Solver_GS_free(in->GS);
39          MEMFREE(in->inv_A_FF);          MEMFREE(in->inv_A_FF);
40          MEMFREE(in->A_FF_pivot);          MEMFREE(in->A_FF_pivot);
41          Paso_SparseMatrix_free(in->A_FC);          Paso_SparseMatrix_free(in->A_FC);
# Line 100  Paso_Solver_AMG* Paso_Solver_getAMG(Paso Line 101  Paso_Solver_AMG* Paso_Solver_getAMG(Paso
101    out->x_C=NULL;    out->x_C=NULL;
102    out->b_C=NULL;    out->b_C=NULL;
103    out->A=Paso_SparseMatrix_getReference(A_p);    out->A=Paso_SparseMatrix_getReference(A_p);
104      out->GS=Paso_Solver_getGS(A_p,verbose);
105      out->GS->sweeps=2;
106    out->level=level;    out->level=level;
107        
108    if ( !(Paso_checkPtr(mis_marker) || Paso_checkPtr(out) || Paso_checkPtr(counter) ) ) {    if ( !(Paso_checkPtr(mis_marker) || Paso_checkPtr(out) || Paso_checkPtr(counter) ) ) {
# Line 107  Paso_Solver_AMG* Paso_Solver_getAMG(Paso Line 110  Paso_Solver_AMG* Paso_Solver_getAMG(Paso
110       time0=Paso_timer();       time0=Paso_timer();
111       #pragma omp parallel for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
112       for (i=0;i<n;++i) mis_marker[i]=-1;       for (i=0;i<n;++i) mis_marker[i]=-1;
113       Paso_Pattern_RS(A_p,mis_marker,0.25);       /*Paso_Pattern_RS(A_p,mis_marker,0.25);*/
114       /*Paso_Pattern_coup(A_p,mis_marker,0.05);*/       Paso_Pattern_coup(A_p,mis_marker,1/n);
115       time2=Paso_timer()-time0;       time2=Paso_timer()-time0;
116       if (Paso_noError()) {       if (Paso_noError()) {
117          #pragma omp parallel for private(i) schedule(static)          #pragma omp parallel for private(i) schedule(static)
# Line 185  Paso_Solver_AMG* Paso_Solver_getAMG(Paso Line 188  Paso_Solver_AMG* Paso_Solver_getAMG(Paso
188             if( Paso_noError()) {             if( Paso_noError()) {
189                /* if there are no nodes in the coarse level there is no more work to do */                /* if there are no nodes in the coarse level there is no more work to do */
190                out->n_C=n-out->n_F;                out->n_C=n-out->n_F;
191                 if (level>0) {                /*if (level>0) {*/
192                   if (out->n_C>10) {
193                     out->rows_in_C=MEMALLOC(out->n_C,index_t);                     out->rows_in_C=MEMALLOC(out->n_C,index_t);
194                     out->mask_C=MEMALLOC(n,index_t);                     out->mask_C=MEMALLOC(n,index_t);
195                     if (! (Paso_checkPtr(out->mask_C) || Paso_checkPtr(out->rows_in_C) ) ) {                     if (! (Paso_checkPtr(out->mask_C) || Paso_checkPtr(out->rows_in_C) ) ) {
# Line 241  Paso_Solver_AMG* Paso_Solver_getAMG(Paso Line 245  Paso_Solver_AMG* Paso_Solver_getAMG(Paso
245                                  /* update A_CC block to get Schur complement and then apply AMG to it */                                  /* update A_CC block to get Schur complement and then apply AMG to it */
246                                  Paso_Solver_updateIncompleteSchurComplement(schur_withFillIn,out->A_CF,out->inv_A_FF,out->A_FF_pivot,out->A_FC);                                  Paso_Solver_updateIncompleteSchurComplement(schur_withFillIn,out->A_CF,out->inv_A_FF,out->A_FF_pivot,out->A_FC);
247                                  time1=Paso_timer()-time1;                                  time1=Paso_timer()-time1;
248                                  out->AMG_of_Schur=Paso_Solver_getAMG(schur_withFillIn,verbose,level-1);                                  out->AMG_of_Schur=Paso_Solver_getAMG(schur_withFillIn,verbose,level+1);
249                                                                    
250                                  /*Paso_SparseMatrix_free(schur);*/                                  /*Paso_SparseMatrix_free(schur);*/
251                                  /* Paso_SparseMatrix_free(schur_withFillIn);*/                                  /* Paso_SparseMatrix_free(schur_withFillIn);*/
# Line 327  void Paso_Solver_solveAMG(Paso_Solver_AM Line 331  void Paso_Solver_solveAMG(Paso_Solver_AM
331       dim_t i,k,oldsweeps;       dim_t i,k,oldsweeps;
332       dim_t n_block=amg->n_block;       dim_t n_block=amg->n_block;
333       double *r=MEMALLOC(amg->n,double);       double *r=MEMALLOC(amg->n,double);
334       Paso_Solver_GS* GS=NULL;       /*Paso_Solver_GS* GS=NULL;*/
335       double *bold=MEMALLOC(amg->n*amg->n_block,double);       double *bold=MEMALLOC(amg->n*amg->n_block,double);
336       double *bnew=MEMALLOC(amg->n*amg->n_block,double);       double *bnew=MEMALLOC(amg->n*amg->n_block,double);
       
      if (amg->level==0) {  
         Paso_UMFPACK1(amg->A,x,b,1);  
      } else {  
337    
338         /*if (amg->level==0) {*/
339         if (amg->n_C<=10) {
340            Paso_UMFPACK1(amg->A,x,b,0);
341         } else {
342          /* presmoothing */          /* presmoothing */
343           GS=Paso_Solver_getGS(amg->A,-1);           Paso_Solver_solveGS(amg->GS,x,b);
344           Paso_Solver_solveGS(GS,x,b);           oldsweeps=amg->GS->sweeps;
345           oldsweeps=GS->sweeps;           if (amg->GS->sweeps>1) {
          if (GS->sweeps>1) {  
346                        
347             #pragma omp parallel for private(i) schedule(static)             #pragma omp parallel for private(i) schedule(static)
348             for (i=0;i<GS->n*GS->n_block;++i) bold[i]=b[i];             for (i=0;i<amg->GS->n*amg->GS->n_block;++i) bold[i]=b[i];
349                        
350             while(GS->sweeps>1) {             while(amg->GS->sweeps>1) {
351                 #pragma omp parallel for private(i) schedule(static)                 #pragma omp parallel for private(i) schedule(static)
352                 for (i=0;i<GS->n*GS->n_block;++i) bnew[i]=bold[i]+b[i];                 for (i=0;i<amg->GS->n*amg->GS->n_block;++i) bnew[i]=bold[i]+b[i];
353                  /* Compute the residual b=b-Ax*/                  /* Compute the residual b=b-Ax*/
354                 Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), amg->A, x, DBLE(1), bnew);                 Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), amg->A, x, DBLE(1), bnew);
355                 /* Go round again*/                 /* Go round again*/
356                 Paso_Solver_solveGS(GS,x,bnew);                 Paso_Solver_solveGS(amg->GS,x,bnew);
357                 #pragma omp parallel for private(i) schedule(static)                 #pragma omp parallel for private(i) schedule(static)
358                 for (i=0;i<GS->n*GS->n_block;++i) bold[i]=bnew[i];                 for (i=0;i<amg->GS->n*amg->GS->n_block;++i) bold[i]=bnew[i];
359                 GS->sweeps=GS->sweeps-1;                 amg->GS->sweeps=amg->GS->sweeps-1;
360             }             }
361             }             }
362             GS->sweeps=oldsweeps;             amg->GS->sweeps=oldsweeps;
363          /* end of presmoothing */          /* end of presmoothing */
364            
365           #pragma omp parallel for private(i) schedule(static)           #pragma omp parallel for private(i) schedule(static)
366           for (i=0;i<amg->n;++i) r[i]=b[i];           for (i=0;i<amg->n;++i) r[i]=b[i];
367                    
# Line 416  void Paso_Solver_solveAMG(Paso_Solver_AM Line 419  void Paso_Solver_solveAMG(Paso_Solver_AM
419          }          }
420    
421       /*postsmoothing*/       /*postsmoothing*/
422       Paso_Solver_solveGS1(GS,x,b);       Paso_Solver_solveGS1(amg->GS,x,b);
423       if (GS->sweeps>1) {       if (amg->GS->sweeps>1) {
424             #pragma omp parallel for private(i) schedule(static)             #pragma omp parallel for private(i) schedule(static)
425             for (i=0;i<GS->n*GS->n_block;++i) bold[i]=b[i];             for (i=0;i<amg->GS->n*amg->GS->n_block;++i) bold[i]=b[i];
426                        
427             while(GS->sweeps>1) {             while(amg->GS->sweeps>1) {
428                 #pragma omp parallel for private(i) schedule(static)                 #pragma omp parallel for private(i) schedule(static)
429                 for (i=0;i<GS->n*GS->n_block;++i) bnew[i]=bold[i]+b[i];                 for (i=0;i<amg->GS->n*amg->GS->n_block;++i) bnew[i]=bold[i]+b[i];
430                  /* Compute the residual b=b-Ax*/                  /* Compute the residual b=b-Ax*/
431                 Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), amg->A, x, DBLE(1), bnew);                 Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), amg->A, x, DBLE(1), bnew);
432                 /* Go round again*/                 /* Go round again*/
433                 Paso_Solver_solveGS(GS,x,bnew);                 Paso_Solver_solveGS1(amg->GS,x,bnew);
434                 #pragma omp parallel for private(i) schedule(static)                 #pragma omp parallel for private(i) schedule(static)
435                 for (i=0;i<GS->n*GS->n_block;++i) bold[i]=bnew[i];                 for (i=0;i<amg->GS->n*amg->GS->n_block;++i) bold[i]=bnew[i];
436                 GS->sweeps=GS->sweeps-1;                 amg->GS->sweeps=amg->GS->sweeps-1;
437             }             }
438             }             }
439       /*end of postsmoothing*/       /*end of postsmoothing*/
440            
441       Paso_Solver_GS_free(GS);       /*Paso_Solver_GS_free(GS);*/
442       }       }
443       MEMFREE(r);       MEMFREE(r);
444       MEMFREE(bold);       MEMFREE(bold);

Legend:
Removed from v.1938  
changed lines
  Added in v.1939

  ViewVC Help
Powered by ViewVC 1.1.26