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

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

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

revision 3009 by jfenwick, Thu Jan 28 02:03:15 2010 UTC revision 3010 by artak, Tue Apr 27 05:10:46 2010 UTC
# Line 79  void Paso_Pattern_Read(char *fileName,di Line 79  void Paso_Pattern_Read(char *fileName,di
79       for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=1);       for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=1);
80  }  }
81    
82    void Paso_Pattern_Write(char *fileName,dim_t n,index_t* mis_marker) {
83        
84        dim_t i;
85        
86        FILE *fileHandle_p = NULL;
87        
88        fileHandle_p = fopen( fileName, "w+" );
89        if( fileHandle_p == NULL )
90        {
91            Paso_setError(IO_ERROR, "Paso_Pattern_Write: Cannot open file for writing.");
92            return;
93        }
94        
95        for (i=0;i<n;++i) {    
96            fprintf(fileHandle_p, "%d\n", mis_marker[i]);
97        }
98        
99        fclose(fileHandle_p);
100        
101    }
102    
103    
104  void Paso_Pattern_YS(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {  void Paso_Pattern_YS(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {
105    
106    dim_t i,j;    dim_t i,j,bi;
107    /*double sum;*/    /*double sum;*/
108    index_t iptr;    index_t iptr;
109    /*index_t *index,*where_p;*/    /*index_t *index,*where_p;*/
110    bool_t passed=FALSE;    bool_t passed=FALSE;
111    dim_t n=A->numRows;    dim_t n=A->numRows;
112      dim_t block_size=A->row_block_size;
113    double *diags;    double *diags;
114      double fnorm;
115    diags=MEMALLOC(n,double);    diags=MEMALLOC(n,double);
116    
117    
118    if (A->pattern->type & PATTERN_FORMAT_SYM) {    if (A->pattern->type & PATTERN_FORMAT_SYM) {
119      Paso_setError(TYPE_ERROR,"Paso_Pattern_coup: symmetric matrix pattern is not supported yet");      Paso_setError(TYPE_ERROR,"Paso_Pattern_coup: symmetric matrix pattern is not supported yet");
120      return;      return;
# Line 102  void Paso_Pattern_YS(Paso_SparseMatrix* Line 125  void Paso_Pattern_YS(Paso_SparseMatrix*
125          if(mis_marker[i]==IS_AVAILABLE)          if(mis_marker[i]==IS_AVAILABLE)
126                      mis_marker[i]=IS_IN_C;                      mis_marker[i]=IS_IN_C;
127    
128      #pragma omp parallel for private(i,j) schedule(static)      #pragma omp parallel for private(i,j,fnorm,bi) schedule(static)
129      for (i=0;i<n;++i) {      for (i=0;i<n;++i) {
130           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) {
131              j=A->pattern->index[iptr];              j=A->pattern->index[iptr];
132              if(i==j) {              if(i==j) {
133                  diags[i]=A->val[iptr];                  if(block_size==1) {
134                      diags[i]=A->val[iptr];    
135                    } else {
136                       fnorm=0;
137                       for(bi=0;bi<block_size*block_size;++bi)
138                       {
139                        fnorm+=A->val[iptr*block_size*block_size+bi]*A->val[iptr*block_size*block_size+bi];
140                       }
141                       fnorm=sqrt(fnorm);
142                       diags[i]=fnorm;
143                    }
144                    
145              }              }
146          }          }
147      }      }
# Line 117  void Paso_Pattern_YS(Paso_SparseMatrix* Line 151  void Paso_Pattern_YS(Paso_SparseMatrix*
151        if (mis_marker[i]==IS_IN_C) {        if (mis_marker[i]==IS_IN_C) {
152          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) {
153               j=A->pattern->index[iptr];               j=A->pattern->index[iptr];
154               if (j!=i && ABS(A->val[iptr])>=threshold*ABS(diags[i])) {               if(block_size==1) {
155                  mis_marker[j]=IS_IN_F;                  if (j!=i && ABS(A->val[iptr])>=threshold*ABS(diags[i])) {
156                       mis_marker[j]=IS_IN_F;
157                    }
158                 } else {
159                       if (j!=i) {
160                            fnorm=0;
161                            for(bi=0;bi<block_size*block_size;++bi)
162                            {
163                             fnorm+=A->val[iptr*block_size*block_size+bi]*A->val[iptr*block_size*block_size+bi];
164                            }
165                            fnorm=sqrt(fnorm);
166                            if (fnorm>=threshold*diags[i]) {
167                               mis_marker[j]=IS_IN_F;
168                            }
169                       }
170               }               }
171          }          }
172        }        }
173      }      }
       
       
174            
175        /*This loop cannot be parallelized, as order matters here.*/        /*This loop cannot be parallelized, as order matters here.*/
176      for (i=0;i<n;i++) {      for (i=0;i<n;i++) {
# Line 133  void Paso_Pattern_YS(Paso_SparseMatrix* Line 179  void Paso_Pattern_YS(Paso_SparseMatrix*
179             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) {
180                j=A->pattern->index[iptr];                j=A->pattern->index[iptr];
181                if (mis_marker[j]==IS_IN_C) {                if (mis_marker[j]==IS_IN_C) {
182                  if ((A->val[iptr]/diags[i])>=-threshold) {                  if(block_size==1) {
183                      passed=TRUE;                      if ((A->val[iptr]/diags[i])>=-threshold) {
184                  }                          passed=TRUE;
185                  else {                      }
186                      passed=FALSE;                      else {
187                      break;                          passed=FALSE;
188                            break;
189                        }
190                    } else {
191                        fnorm=0;
192                        for(bi=0;bi<block_size*block_size;++bi)
193                        {
194                           fnorm+=A->val[iptr*block_size*block_size+bi]*A->val[iptr*block_size*block_size+bi];
195                        }
196                        fnorm=sqrt(fnorm);
197                        if ((fnorm/diags[i])<=threshold) {
198                            passed=TRUE;
199                        }
200                        else {
201                            passed=FALSE;
202                            break;
203                        }
204                  }                  }
205                }                      
206                  }
207                  
208             }             }
209             if (passed) mis_marker[i]=IS_IN_C;             if (passed) mis_marker[i]=IS_IN_C;
210          }          }
# Line 160  void Paso_Pattern_YS(Paso_SparseMatrix* Line 224  void Paso_Pattern_YS(Paso_SparseMatrix*
224   */   */
225  void Paso_Pattern_RS(Paso_SparseMatrix* A, index_t* mis_marker, double theta)  void Paso_Pattern_RS(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
226  {  {
227    dim_t i,n,j;    dim_t i,n,j,bi;
228    index_t iptr;    index_t iptr;
229    double threshold,max_offdiagonal;    double threshold,max_offdiagonal;
230      double fnorm;
231      dim_t block_size=A->row_block_size;
232        
233    Paso_Pattern *out=NULL;    Paso_Pattern *out=NULL;
234        
# Line 183  void Paso_Pattern_RS(Paso_SparseMatrix* Line 249  void Paso_Pattern_RS(Paso_SparseMatrix*
249      Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");      Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
250      return;      return;
251    }    }
252      /*#pragma omp parallel for private(i,iptr,max_offdiagonal,threshold,j) schedule(static)*/      /*#pragma omp parallel for private(i,iptr,max_offdiagonal,threshold,j,fnorm,bi) schedule(static)*/
253      for (i=0;i<n;++i) {      for (i=0;i<n;++i) {
254        if(mis_marker[i]==IS_AVAILABLE) {        if(mis_marker[i]==IS_AVAILABLE) {
255          max_offdiagonal = DBL_MIN;          max_offdiagonal = DBL_MIN;
256          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) {
257              if(A->pattern->index[iptr] != i){              if(A->pattern->index[iptr] != i){
258                    if(block_size==1) {
259                    max_offdiagonal = MAX(max_offdiagonal,ABS(A->val[iptr]));                    max_offdiagonal = MAX(max_offdiagonal,ABS(A->val[iptr]));
260                    } else {
261                        fnorm=0;
262                        for(bi=0;bi<block_size*block_size;++bi)
263                        {
264                           fnorm+=A->val[iptr*block_size*block_size+bi]*A->val[iptr*block_size*block_size+bi];
265                        }
266                        fnorm=sqrt(fnorm);
267                        max_offdiagonal = MAX(max_offdiagonal,fnorm);
268                    }
269              }              }
270          }          }
271                    
272          threshold = theta*max_offdiagonal;          threshold = theta*max_offdiagonal;
273          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) {
274              j=A->pattern->index[iptr];              j=A->pattern->index[iptr];
275              if(ABS(A->val[iptr])>=threshold && i!=j) {               if(block_size==1) {
276                  Paso_IndexList_insertIndex(&(index_list[i]),j);                      if(ABS(A->val[iptr])>=threshold && i!=j) {
277                  /*Paso_IndexList_insertIndex(&(index_list[j]),i);*/                          Paso_IndexList_insertIndex(&(index_list[i]),j);
278              }                          /*Paso_IndexList_insertIndex(&(index_list[j]),i);*/
279                        }
280                 } else {
281                        if (i!=j) {
282                            fnorm=0;
283                            for(bi=0;bi<block_size*block_size;++bi)
284                            {
285                               fnorm+=A->val[iptr*block_size*block_size+bi]*A->val[iptr*block_size*block_size+bi];
286                            }
287                            fnorm=sqrt(fnorm);
288                            if(fnorm>=threshold) {
289                                Paso_IndexList_insertIndex(&(index_list[i]),j);
290                                /*Paso_IndexList_insertIndex(&(index_list[j]),i);*/
291                            }
292                        }
293                 }
294          }          }
295         }         }
296        }        }
# Line 221  void Paso_Pattern_RS(Paso_SparseMatrix* Line 312  void Paso_Pattern_RS(Paso_SparseMatrix*
312    
313  void Paso_Pattern_Aggregiation(Paso_SparseMatrix* A, index_t* mis_marker, double theta)  void Paso_Pattern_Aggregiation(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
314  {  {
315    dim_t i,j,n;    dim_t i,j,n,bi;
316    index_t iptr;    index_t iptr;
317    double diag,eps_Aii,val;    double diag,eps_Aii,val;
318    double* diags;    double* diags;
319      double fnorm;
320      dim_t block_size=A->row_block_size;
321    
322    
323    Paso_Pattern *out=NULL;    Paso_Pattern *out=NULL;
# Line 248  void Paso_Pattern_Aggregiation(Paso_Spar Line 341  void Paso_Pattern_Aggregiation(Paso_Spar
341    }    }
342    
343    
344      #pragma omp parallel for private(i,iptr,diag) schedule(static)      #pragma omp parallel for private(i,iptr,diag,fnorm,bi) schedule(static)
345        for (i=0;i<n;++i) {        for (i=0;i<n;++i) {
346          diag = 0;          diag = 0;
347          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) {
348              if(A->pattern->index[iptr] == i){              if(A->pattern->index[iptr] == i){
349                  diag+=A->val[iptr];                  if(block_size==1) {
350                        diag+=A->val[iptr];
351                    } else {
352                        fnorm=0;
353                        for(bi=0;bi<block_size*block_size;++bi)
354                        {
355                           fnorm+=A->val[iptr*block_size*block_size+bi]*A->val[iptr*block_size*block_size+bi];
356                        }
357                        fnorm=sqrt(fnorm);
358                        diag+=fnorm;
359                    }
360                    
361              }              }
362          }          }
363          diags[i]=ABS(diag);          diags[i]=ABS(diag);
364        }        }
365    
366    
367      #pragma omp parallel for private(i,iptr,j,val,eps_Aii) schedule(static)      #pragma omp parallel for private(i,iptr,j,val,eps_Aii,fnorm,bi) schedule(static)
368       for (i=0;i<n;++i) {       for (i=0;i<n;++i) {
369         if (mis_marker[i]==IS_AVAILABLE) {         if (mis_marker[i]==IS_AVAILABLE) {
370          eps_Aii = theta*theta*diags[i];          eps_Aii = theta*theta*diags[i];
371          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) {
372              j=A->pattern->index[iptr];              j=A->pattern->index[iptr];
373              val=A->val[iptr];              if(block_size==1) {
374                    val=A->val[iptr];
375                } else {
376                    fnorm=0;
377                    for(bi=0;bi<block_size*block_size;++bi)
378                    {
379                       fnorm+=A->val[iptr*block_size*block_size+bi]*A->val[iptr*block_size*block_size+bi];
380                    }
381                    fnorm=sqrt(fnorm);
382                    val=fnorm;
383                }
384                
385                if((val*val)>=(eps_Aii*diags[j])) {                if((val*val)>=(eps_Aii*diags[j])) {
386                 Paso_IndexList_insertIndex(&(index_list[i]),j);                 Paso_IndexList_insertIndex(&(index_list[i]),j);
387                }                }
# Line 713  void Paso_Pattern_Standard(Paso_SparseMa Line 828  void Paso_Pattern_Standard(Paso_SparseMa
828          threshold = theta*max_offdiagonal;          threshold = theta*max_offdiagonal;
829          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) {
830              j=A->pattern->index[iptr];              j=A->pattern->index[iptr];
831              if(ABS(A->val[iptr])>=threshold && i!=j) {              if(ABS(A->val[iptr])>threshold && i!=j) {
832                  Paso_IndexList_insertIndex(&(index_list[i]),j);                  Paso_IndexList_insertIndex(&(index_list[i]),j);
833              }              }
834          }          }
# Line 799  void Paso_Pattern_Standard(Paso_SparseMa Line 914  void Paso_Pattern_Standard(Paso_SparseMa
914                  }                  }
915              }              }
916          }          }
917          for (iptr=S->ptr[i];iptr<S->ptr[i+1]; ++iptr) {         for (iptr=S->ptr[i];iptr<S->ptr[i+1]; ++iptr) {
918              j=S->index[iptr];              j=S->index[iptr];
919              if(mis_marker[j]==IS_AVAILABLE) {              if(mis_marker[j]==IS_AVAILABLE) {
920                             lambda[j]--;                             lambda[j]--;
921              }              }
922          }          }
               
923      }      }
924            
925     /* Used when transpose of S is not available */     /* Used when transpose of S is not available */
# Line 1047  Paso_Pattern* Paso_Pattern_getTranspose( Line 1161  Paso_Pattern* Paso_Pattern_getTranspose(
1161  }  }
1162    
1163    
1164    /************** BLOCK COARSENENING *********************/
1165    
1166    void Paso_Pattern_Standard_Block(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
1167    {
1168      dim_t i,n,j,k;
1169      index_t iptr,jptr;
1170      /*index_t *index,*where_p;*/
1171      double threshold,max_offdiagonal;
1172      dim_t *lambda;   /*mesure of importance */
1173      dim_t maxlambda=0;
1174      index_t index_maxlambda=0;
1175      double time0=0;
1176      bool_t verbose=0;
1177      dim_t n_block=A->row_block_size;
1178      
1179      double fnorm=0;
1180      dim_t bi;
1181      
1182      Paso_Pattern *S=NULL;
1183      Paso_Pattern *ST=NULL;
1184      Paso_IndexList* index_list=NULL;
1185      
1186      /*dim_t lk;*/
1187    
1188      index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);
1189       if (! Paso_checkPtr(index_list)) {
1190            #pragma omp parallel for private(i) schedule(static)
1191            for(i=0;i<A->pattern->numOutput;++i) {
1192                 index_list[i].extension=NULL;
1193                 index_list[i].n=0;
1194            }
1195        }
1196      
1197      
1198      n=A->numRows;
1199      if (A->pattern->type & PATTERN_FORMAT_SYM) {
1200        Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
1201        return;
1202      }
1203      
1204        time0=Paso_timer();
1205       k=0;
1206       /*Paso_Pattern_getReport(n,mis_marker);*/
1207       /*printf("Blocks %d %d\n",n_block,A->len);*/
1208      
1209        /*S_i={j \in N_i; i strongly coupled to j}*/
1210        #pragma omp parallel for private(i,bi,fnorm,iptr,max_offdiagonal,threshold,j) schedule(static)
1211        for (i=0;i<n;++i) {
1212          if(mis_marker[i]==IS_AVAILABLE) {
1213            max_offdiagonal = DBL_MIN;
1214            for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
1215                j=A->pattern->index[iptr];
1216                if( j != i){
1217                    fnorm=0;
1218                    for(bi=0;bi<n_block*n_block;++bi)
1219                       {
1220                        fnorm+=A->val[iptr*n_block*n_block+bi]*A->val[iptr*n_block*n_block+bi];
1221                       }
1222                    fnorm=sqrt(fnorm);
1223                    max_offdiagonal = MAX(max_offdiagonal,fnorm);
1224                }
1225            }
1226            threshold = theta*max_offdiagonal;
1227            for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
1228                j=A->pattern->index[iptr];
1229                if( j != i){
1230                    fnorm=0;
1231                    for(bi=0;bi<n_block*n_block;++bi)
1232                       {
1233                        fnorm+=A->val[iptr*n_block*n_block+bi]*A->val[iptr*n_block*n_block+bi];
1234                       }
1235                    fnorm=sqrt(fnorm);
1236                if(fnorm>=threshold) {
1237                    Paso_IndexList_insertIndex(&(index_list[i]),j);
1238                }
1239                }
1240            }
1241          }
1242        }
1243      
1244      S=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
1245      ST=Paso_Pattern_getTranspose(S);
1246      
1247      /*printf("Patterns len %d %d\n",S->len,ST->len);*/
1248        
1249      time0=Paso_timer()-time0;
1250      if (verbose) fprintf(stdout,"timing: RS filtering and pattern creation: %e\n",time0);
1251      
1252      lambda=TMPMEMALLOC(n,dim_t);
1253      
1254      #pragma omp parallel for private(i) schedule(static)
1255      for (i=0;i<n;++i) { lambda[i]=IS_NOT_AVAILABLE; }
1256      
1257      /*S_i={j \in N_i; i strongly coupled to j}*/
1258    
1259      /*
1260      #pragma omp parallel for private(i,iptr,lk) schedule(static)
1261      for (i=0;i<n;++i) {
1262            lk=0;
1263            for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
1264                if(ABS(A->val[iptr])>1.e-15 && A->pattern->index[iptr]!=i )
1265                   lk++;
1266            }
1267            #pragma omp critical
1268            k+=lk;
1269            if(k==0) {
1270                mis_marker[i]=IS_IN_F;
1271            }
1272      }
1273      */
1274    
1275    
1276      k=0;
1277      maxlambda=0;
1278      
1279      time0=Paso_timer();
1280      
1281        for (i=0;i<n;++i) {
1282          if(mis_marker[i]==IS_AVAILABLE) {
1283            lambda[i]=how_many(k,ST,FALSE);   /* if every row is available then k and i are the same.*/
1284            /*lambda[i]=how_many(i,S,TRUE);*/
1285            /*printf("lambda[%d]=%d, k=%d \n",i,lambda[i],k);*/
1286            k++;
1287            if(maxlambda<lambda[i]) {
1288                maxlambda=lambda[i];
1289                index_maxlambda=i;
1290            }
1291          }
1292        }
1293    
1294      k=0;
1295      time0=Paso_timer()-time0;
1296      if (verbose) fprintf(stdout,"timing: Lambdas computations at the begining: %e\n",time0);
1297      
1298      time0=Paso_timer();
1299      
1300      /*Paso_Pattern_getReport(n,mis_marker);*/
1301      
1302      while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
1303        
1304        if(index_maxlambda<0) {
1305            break;
1306        }
1307        
1308        i=index_maxlambda;
1309        if(mis_marker[i]==IS_AVAILABLE) {
1310            mis_marker[index_maxlambda]=IS_IN_C;
1311            lambda[index_maxlambda]=IS_NOT_AVAILABLE;
1312            for (iptr=ST->ptr[i];iptr<ST->ptr[i+1]; ++iptr) {
1313                j=ST->index[iptr];
1314                if(mis_marker[j]==IS_AVAILABLE) {
1315                    mis_marker[j]=IS_IN_F;
1316                    lambda[j]=IS_NOT_AVAILABLE;
1317                    for (jptr=S->ptr[j];jptr<S->ptr[j+1]; ++jptr) {
1318                            k=S->index[jptr];
1319                            if(mis_marker[k]==IS_AVAILABLE) {
1320                               lambda[k]++;
1321                            }
1322                    }
1323                }
1324            }
1325           for (iptr=S->ptr[i];iptr<S->ptr[i+1]; ++iptr) {
1326                j=S->index[iptr];
1327                if(mis_marker[j]==IS_AVAILABLE) {
1328                               lambda[j]--;
1329                }
1330            }
1331        }
1332        
1333       /* Used when transpose of S is not available */
1334       /*
1335        for (i=0;i<n;++i) {
1336            if(mis_marker[i]==IS_AVAILABLE) {
1337                if (i==index_maxlambda) {
1338                    mis_marker[index_maxlambda]=IS_IN_C;
1339                    lambda[index_maxlambda]=-1;
1340                    for (j=0;j<n;++j) {
1341                        if(mis_marker[j]==IS_AVAILABLE) {
1342                            index=&(S->index[S->ptr[j]]);
1343                            where_p=(index_t*)bsearch(&i,
1344                                            index,
1345                                            S->ptr[j + 1]-S->ptr[j],
1346                                            sizeof(index_t),
1347                                            Paso_comparIndex);
1348                            if (where_p!=NULL) {
1349                                mis_marker[j]=IS_IN_F;
1350                                lambda[j]=-1;
1351                                for (iptr=S->ptr[j];iptr<S->ptr[j+1]; ++iptr) {
1352                                    k=S->index[iptr];
1353                                    if(mis_marker[k]==IS_AVAILABLE) {
1354                                       lambda[k]++;
1355                                    }
1356                                }
1357                            }
1358                            
1359                        }
1360                    }
1361                    
1362                }
1363            }
1364        }
1365       */
1366       index_maxlambda=arg_max(n,lambda, IS_NOT_AVAILABLE);
1367      }
1368      
1369      time0=Paso_timer()-time0;
1370      if (verbose) fprintf(stdout,"timing: Loop : %e\n",time0);
1371    
1372      /*Paso_Pattern_getReport(n,mis_marker);*/
1373      
1374      #pragma omp parallel for private(i) schedule(static)
1375      for (i=0;i<n;++i)
1376          if(mis_marker[i]==IS_AVAILABLE) {
1377            mis_marker[i]=IS_IN_F;
1378           }
1379          
1380      /*Paso_Pattern_getReport(n,mis_marker);*/
1381    
1382        TMPMEMFREE(lambda);
1383        
1384         /* clean up */
1385        if (index_list!=NULL) {
1386            #pragma omp parallel for private(i) schedule(static)
1387            for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
1388         }
1389    
1390        TMPMEMFREE(index_list);
1391        Paso_Pattern_free(S);
1392    
1393        /* swap to TRUE/FALSE in mis_marker */
1394         #pragma omp parallel for private(i) schedule(static)
1395         for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_F);
1396        
1397    }
1398    
1399    
1400    
1401  #undef IS_AVAILABLE  #undef IS_AVAILABLE
1402  #undef IS_NOT_AVAILABLE  #undef IS_NOT_AVAILABLE
1403  #undef IS_IN_F  #undef IS_IN_F

Legend:
Removed from v.3009  
changed lines
  Added in v.3010

  ViewVC Help
Powered by ViewVC 1.1.26