/[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 2555 by gross, Thu Jul 23 09:19:15 2009 UTC revision 2556 by artak, Mon Jul 27 02:05:20 2009 UTC
# Line 43  void Paso_Solver_AMG_free(Paso_Solver_AM Line 43  void Paso_Solver_AMG_free(Paso_Solver_AM
43          MEMFREE(in->A_FF_pivot);          MEMFREE(in->A_FF_pivot);
44          Paso_SparseMatrix_free(in->A_FC);          Paso_SparseMatrix_free(in->A_FC);
45          Paso_SparseMatrix_free(in->A_CF);          Paso_SparseMatrix_free(in->A_CF);
46            Paso_SparseMatrix_free(in->A);
47          MEMFREE(in->rows_in_F);          MEMFREE(in->rows_in_F);
48          MEMFREE(in->rows_in_C);          MEMFREE(in->rows_in_C);
49          MEMFREE(in->mask_F);          MEMFREE(in->mask_F);
# Line 84  to Line 85  to
85  */  */
86  Paso_Solver_AMG* Paso_Solver_getAMG(Paso_SparseMatrix *A_p,dim_t level,Paso_Options* options) {  Paso_Solver_AMG* Paso_Solver_getAMG(Paso_SparseMatrix *A_p,dim_t level,Paso_Options* options) {
87    Paso_Solver_AMG* out=NULL;    Paso_Solver_AMG* out=NULL;
88      Paso_Pattern* temp1=NULL;
89      Paso_Pattern* temp2=NULL;
90    bool_t verbose=options->verbose;    bool_t verbose=options->verbose;
91    dim_t n=A_p->numRows;    dim_t n=A_p->numRows;
92    dim_t n_block=A_p->row_block_size;    dim_t n_block=A_p->row_block_size;
# Line 94  Paso_Solver_AMG* Paso_Solver_getAMG(Paso Line 97  Paso_Solver_AMG* Paso_Solver_getAMG(Paso
97    Paso_SparseMatrix * schur=NULL;    Paso_SparseMatrix * schur=NULL;
98    Paso_SparseMatrix * schur_withFillIn=NULL;    Paso_SparseMatrix * schur_withFillIn=NULL;
99    double S=0;    double S=0;
100      char fname[6];
101        
102    /*Make sure we have block sizes 1*/    /*Make sure we have block sizes 1*/
103    if (A_p->col_block_size>1) {    if (A_p->col_block_size>1) {
# Line 128  Paso_Solver_AMG* Paso_Solver_getAMG(Paso Line 132  Paso_Solver_AMG* Paso_Solver_getAMG(Paso
132       out->solver=NULL;       out->solver=NULL;
133       /*out->GS=Paso_Solver_getGS(A_p,verbose);*/       /*out->GS=Paso_Solver_getGS(A_p,verbose);*/
134       out->level=level;       out->level=level;
135        
136           sprintf(fname,"A%d.mat",level);
137      
138           Paso_SparseMatrix_saveMM(A_p,fname);
139        
140       if (level==0 || n<=options->min_coarse_matrix_size) {       if (level==0 || n<=options->min_coarse_matrix_size) {
141           out->coarsest_level=TRUE;           out->coarsest_level=TRUE;
# Line 157  Paso_Solver_AMG* Paso_Solver_getAMG(Paso Line 165  Paso_Solver_AMG* Paso_Solver_getAMG(Paso
165          }          }
166          else {          else {
167             /*Default coarseneing*/             /*Default coarseneing*/
168             /* Paso_Pattern_RS(A_p,mis_marker,options->coarsening_threshold); */              Paso_Pattern_RS(A_p,mis_marker,options->coarsening_threshold);
169    
170               Paso_Pattern_Aggregiation(A_p,mis_marker,options->coarsening_threshold);               /*Paso_Pattern_Aggregiation(A_p,mis_marker,options->coarsening_threshold);*/
171          }          }
172            
173          if (Paso_noError()) {          if (Paso_noError()) {
# Line 191  Paso_Solver_AMG* Paso_Solver_getAMG(Paso Line 199  Paso_Solver_AMG* Paso_Solver_getAMG(Paso
199                #pragma omp parallel for private(i,iPtr,j,S) schedule(static)                #pragma omp parallel for private(i,iPtr,j,S) schedule(static)
200                for (i = 0; i < out->n_F; ++i) {                for (i = 0; i < out->n_F; ++i) {
201                  S=0;                  S=0;
202  printf("%d : %d -> ",i, out->rows_in_F[i]);  /*printf("%d : %d -> ",i, out->rows_in_F[i]);*/
203                  for (iPtr=A_p->pattern->ptr[out->rows_in_F[i]];iPtr<A_p->pattern->ptr[out->rows_in_F[i] + 1]; ++iPtr) {                  for (iPtr=A_p->pattern->ptr[out->rows_in_F[i]];iPtr<A_p->pattern->ptr[out->rows_in_F[i] + 1]; ++iPtr) {
204                   j=A_p->pattern->index[iPtr];                   j=A_p->pattern->index[iPtr];
205  if (j==out->rows_in_F[i]) printf("%e",A_p->val[iPtr]);  /*if (j==out->rows_in_F[i]) printf("diagonal %e",A_p->val[iPtr]);*/
206                   if (mis_marker[j])                   if (mis_marker[j] && j==out->rows_in_F[i])
207                       S+=A_p->val[iPtr];                       S=A_p->val[iPtr];
208                  }                  }
209  printf("-> %e\n",S);  /*printf("-> %e\n",S);*/
210                  out->inv_A_FF[i]=1./S;                  out->inv_A_FF[i]=1./S;
211                }                }
212             }             }
# Line 240  printf("-> %e\n",S); Line 248  printf("-> %e\n",S);
248          }          }
249          if ( Paso_noError()) {          if ( Paso_noError()) {
250                 /*find the pattern of the schur complement with fill in*/                 /*find the pattern of the schur complement with fill in*/
251                      temp1=Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,out->A_CF->pattern,out->A_FC->pattern);
252                schur_withFillIn=Paso_SparseMatrix_alloc(A_p->type,Paso_Pattern_binop(PATTERN_FORMAT_DEFAULT, schur->pattern, Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,out->A_CF->pattern,out->A_FC->pattern)),1,1, TRUE);                temp2=Paso_Pattern_binop(PATTERN_FORMAT_DEFAULT, schur->pattern, temp1);
253                  schur_withFillIn=Paso_SparseMatrix_alloc(A_p->type,temp2,1,1, TRUE);
254                  Paso_Pattern_free(temp1);
255                  Paso_Pattern_free(temp2);
256          }          }
257          if ( Paso_noError()) {          if ( Paso_noError()) {
258                /* copy values over*/                /* copy values over*/

Legend:
Removed from v.2555  
changed lines
  Added in v.2556

  ViewVC Help
Powered by ViewVC 1.1.26