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

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

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

revision 2803 by artak, Thu Dec 3 01:51:55 2009 UTC revision 2804 by artak, Fri Dec 4 02:06:51 2009 UTC
# Line 162  void Paso_SparseMatrix_updateWeights(Pas Line 162  void Paso_SparseMatrix_updateWeights(Pas
162    
163    double *alpha;    double *alpha;
164    double *beta;    double *beta;
   double *alpha_den;  
   double *beta_den;  
165    double a_ii=0;    double a_ii=0;
166      double sum_all_neg,sum_all_pos,sum_strong_neg,sum_strong_pos;
167        
168    dim_t n=A->numRows;    dim_t n=A->numRows;
169    dim_t F=W_FC->numRows;    dim_t F=W_FC->numRows;
170    dim_t i,j,k;    dim_t i,j,k;
171        
172    index_t iPtr,*index, *where_p;    index_t iPtr,dptr=0;
173      /*index_t *index, *where_p;*/
174    
175    alpha=TMPMEMALLOC(F,double);    alpha=TMPMEMALLOC(F,double);
176    beta=TMPMEMALLOC(F,double);    beta=TMPMEMALLOC(F,double);
   alpha_den=TMPMEMALLOC(F,double);  
   beta_den=TMPMEMALLOC(F,double);  
177    
178    k=0;    
179      k=0;
180    for (i = 0; i < n; ++i) {    for (i = 0; i < n; ++i) {
181        if(mis_marker[i]) {        if(mis_marker[i]) {
             iPtr=A->pattern->ptr[i];  
             index=&(A->pattern->index[iPtr]);  
             where_p=(index_t*)bsearch(&i,  
                                      index,  
                                      A->pattern->ptr[i + 1]-A->pattern->ptr[i],  
                                      sizeof(index_t),  
                                      Paso_comparIndex);  
             if (where_p!=NULL) {  
                    a_ii=A->val[iPtr+(index_t)(where_p-index)];  
             }  
             else {  
              Paso_setError(VALUE_ERROR,"Paso_Solver_getAMG: Main diagonal element of system matrix is missing.");  
             }  
182              alpha[k]=0;              alpha[k]=0;
183              beta[k]=0;              beta[k]=0;
184              alpha_den[k]=0;              sum_all_neg=0;
185              beta_den[k]=0;              sum_all_pos=0;
186                sum_strong_neg=0;
187                sum_strong_pos=0;
188              /*compute beta_i=sum_{j\inN_i}a^{+}_ij and alpha_i=sum_{j\inN_i}a^{-}_ij and denominators for both alpha_i and beta_i*/              /*compute beta_i=sum_{j\inN_i}a^{+}_ij and alpha_i=sum_{j\inN_i}a^{-}_ij and denominators for both alpha_i and beta_i*/
189              for (iPtr=A->pattern->ptr[i];iPtr<A->pattern->ptr[i + 1]; ++iPtr) {              for (iPtr=A->pattern->ptr[i];iPtr<A->pattern->ptr[i + 1]; ++iPtr) {
190                   j=A->pattern->index[iPtr];                   j=A->pattern->index[iPtr];
191                   if(j!=i) {                   if(j!=i) {
192                          if(A->val[iPtr]<0) {                          if(A->val[iPtr]<0) {
193                            alpha[k]+=A->val[iPtr];                            sum_all_neg+=A->val[iPtr];
194                            if(!mis_marker[j]) {                            if(!mis_marker[j]) {
195                              alpha_den[k]+=A->val[iPtr];                              sum_strong_neg+=A->val[iPtr];
196                            }                            }
197                          }                          }
198                          else if(A->val[iPtr]>0) {                          else if(A->val[iPtr]>0) {
199                            beta[k]+=A->val[iPtr];                            sum_all_pos+=A->val[iPtr];
200                            if(!mis_marker[j]) {                            if(!mis_marker[j]) {
201                              beta_den[k]+=A->val[iPtr];                              sum_strong_pos+=A->val[iPtr];
202                            }                            }
203                          }                          }
204                   }                   }
205                     else {
206                          a_ii=A->val[iPtr];
207                          dptr=iPtr;
208                     }
209               }               }
210              if(alpha_den[k]!=0) {  
211                   alpha[k]=alpha[k]/(alpha_den[k]*a_ii);              if(sum_strong_neg!=0) {
212                     alpha[k]=sum_all_neg/(sum_strong_neg);
213                }
214                if(sum_strong_pos!=0) {
215                     beta[k]=sum_all_pos/(sum_strong_pos);
216              }              }
217              if(beta_den[k]!=0) {  /*            else {
218                   beta[k]=beta[k]/(beta_den[k]*a_ii);                a_ii+=beta[k];
219                  A->val[dptr]=a_ii;
220                  beta[k]=0;
221              }              }
222              k++;  */           alpha[k]=alpha[k]/(a_ii);
223                beta[k]=beta[k]/(a_ii);
224              k++;
225         /*printf("Got in row=%d, alpha[%d]=%e, beta[%d]=%e, a_den=%e, b_den=%e \n",i,k-1,alpha[k-1],k-1,beta[k-1],alpha_den[k-1],beta_den[k-1]);*/         /*printf("Got in row=%d, alpha[%d]=%e, beta[%d]=%e, a_den=%e, b_den=%e \n",i,k-1,alpha[k-1],k-1,beta[k-1],alpha_den[k-1],beta_den[k-1]);*/
226        }        }
227     }     }
# Line 229  void Paso_SparseMatrix_updateWeights(Pas Line 229  void Paso_SparseMatrix_updateWeights(Pas
229        for (i = 0; i < W_FC->numRows; ++i) {        for (i = 0; i < W_FC->numRows; ++i) {
230              for (iPtr=W_FC->pattern->ptr[i];iPtr<W_FC->pattern->ptr[i + 1]; ++iPtr) {              for (iPtr=W_FC->pattern->ptr[i];iPtr<W_FC->pattern->ptr[i + 1]; ++iPtr) {
231                   j=W_FC->pattern->index[iPtr];                   j=W_FC->pattern->index[iPtr];
                  /*if(!mis_marker[j]) {*/  
232                     if(W_FC->val[iPtr]<0) {                     if(W_FC->val[iPtr]<0) {
233                        W_FC->val[iPtr]=-alpha[i]*W_FC->val[iPtr];                        W_FC->val[iPtr]=-alpha[i]*W_FC->val[iPtr];
234                      }                      }
235                     else if(W_FC->val[iPtr]>0) {                     else if(W_FC->val[iPtr]>0) {
236                        W_FC->val[iPtr]=-beta[i]*W_FC->val[iPtr];                        W_FC->val[iPtr]=-beta[i]*W_FC->val[iPtr];
237                     }                     }
                  /*}*/  
238              }              }
239           }           }
240                
241        TMPMEMFREE(alpha);        TMPMEMFREE(alpha);
242        TMPMEMFREE(beta);        TMPMEMFREE(beta);
       TMPMEMFREE(alpha_den);  
       TMPMEMFREE(beta_den);  
243  }  }
244    
245        

Legend:
Removed from v.2803  
changed lines
  Added in v.2804

  ViewVC Help
Powered by ViewVC 1.1.26