[Getfem-commits] r4838 - in /trunk/getfem/src: getfem_models.cc gmm/gmm_range_basis.h

[Yves . Renard]

Author: renard
Date: Sat Dec 20 13:40:22 2014
New Revision: 4838

URL: http://svn.gna.org/viewcvs/getfem?rev=4838&view=rev
Log:
fix a little bug in range_basis

Modified:
    trunk/getfem/src/getfem_models.cc
    trunk/getfem/src/gmm/gmm_range_basis.h

Modified: trunk/getfem/src/getfem_models.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_models.cc?rev=4838&r1=4837&r2=4838&view=diff
==============================================================================
--- trunk/getfem/src/getfem_models.cc   (original)
+++ trunk/getfem/src/getfem_models.cc   Sat Dec 20 13:40:22 2014
@@ -287,7 +287,7 @@
 
       std::vector<std::string> &mults = multipliers[itbd->first];
       VAR_SET::iterator it2 = variables.find(itbd->first);
-
+      
       gmm::col_matrix< gmm::rsvector<scalar_type> > MGLOB;
       if (mults.size() > 1) {
         size_type s = 0;
@@ -297,13 +297,13 @@
         }
         gmm::resize(MGLOB, it2->second.mf->nb_dof(), s);
       }
-
+      
       size_type s = 0;
       std::set<size_type> glob_columns;
       std::vector<dal::bit_vector> mult_kept_dofs;
       for (size_type k = 0; k < mults.size(); ++k) {
         VAR_SET::iterator it = variables.find(mults[k]);
-
+        
         // This step forces the recomputation of corresponding bricks.
         // A test to check if a modification is really necessary could
         // be done first ... (difficult to coordinate with other multipliers)
@@ -320,64 +320,66 @@
           MM(it2->second.associated_mf().nb_dof(), it->second.mf->nb_dof());   
     
         bool termadded = false;
 
-  if (it->second.filter == VDESCRFILTER_CTERM) {
+        if (it->second.filter == VDESCRFILTER_CTERM) {
 
           for (dal::bv_visitor ib(valid_bricks); !ib.finished(); ++ib) {
-      const brick_description &brick = bricks[ib];
-      bool bupd = false;
-      bool cplx = is_complex() && brick.pbr->is_complex();
-      
-      for (size_type j = 0; j < brick.tlist.size(); ++j) {
-        
-        const term_description &term = brick.tlist[j];
-        
-        if (term.is_matrix_term && !mults[k].compare(term.var1) &&
-      !it2->first.compare(term.var2)) {
-    if (!bupd) {
-      brick.terms_to_be_computed = true;
-      update_brick(ib, BUILD_MATRIX);
-      bupd = true;
-    }
-    if (cplx)
-      gmm::add(gmm::transposed(gmm::real_part(brick.cmatlist[j])),
-         MM);
-    else
-      gmm::add(gmm::transposed(brick.rmatlist[j]), MM);
-    termadded = true;
-    
-        } else if (term.is_matrix_term && !mults[k].compare(term.var2) &&
-       !it2->first.compare(term.var1)) {
-    if (!bupd) {
-      brick.terms_to_be_computed = true;
-      update_brick(ib, BUILD_MATRIX);
-      bupd = true;
-    }
-    if (cplx)
-      gmm::add(gmm::real_part(brick.cmatlist[j]), MM);
-    else
-      gmm::add(brick.rmatlist[j], MM);
-    termadded = true;
-        }
-      }
-    }
-
-    if (!termadded)
-      GMM_WARNING1("No term found to filter multiplier " << it->first
-       << ". Variable is cancelled");
-  } else if (it->second.filter == VDESCRFILTER_INFSUP) {
+            const brick_description &brick = bricks[ib];
+            bool bupd = false;
+            bool cplx = is_complex() && brick.pbr->is_complex();
+            
+            for (size_type j = 0; j < brick.tlist.size(); ++j) {
+              
+              const term_description &term = brick.tlist[j];
+              
+              if (term.is_matrix_term && !mults[k].compare(term.var1) &&
+                  !it2->first.compare(term.var2)) {
+                if (!bupd) {
+                  brick.terms_to_be_computed = true;
+                  update_brick(ib, BUILD_MATRIX);
+                  bupd = true;
+                }
+                if (cplx)
+                  gmm::add(gmm::transposed(gmm::real_part(brick.cmatlist[j])),
+                           MM);
+                else
+                  gmm::add(gmm::transposed(brick.rmatlist[j]), MM);
+                termadded = true;
+                
+              } else if (term.is_matrix_term && !mults[k].compare(term.var2) &&
+                         !it2->first.compare(term.var1)) {
+                if (!bupd) {
+                  brick.terms_to_be_computed = true;
+                  update_brick(ib, BUILD_MATRIX);
+                  bupd = true;
+                }
+                if (cplx)
+                  gmm::add(gmm::real_part(brick.cmatlist[j]), MM);
+                else
+                  gmm::add(brick.rmatlist[j], MM);
+                termadded = true;
+              }
+            }
+          }
+          
+          if (!termadded)
+            GMM_WARNING1("No term found to filter multiplier " << it->first
+                         << ". Variable is cancelled");
+        } else if (it->second.filter == VDESCRFILTER_INFSUP) {
           mesh_region rg(it->second.m_region);
           it->second.mim->linked_mesh().intersect_with_mpi_region(rg);
-    asm_mass_matrix(MM, *(it->second.mim), it2->second.associated_mf(),
-        *(it->second.mf), rg);
-  }
-
+          asm_mass_matrix(MM, *(it->second.mim), it2->second.associated_mf(),
+                          *(it->second.mf), rg);
+        }
+        
         MPI_SUM_SPARSE_MATRIX(MM);
-
+        
         //
         // filtering
         //
         std::set<size_type> columns;
         gmm::range_basis(MM, columns);
+        
+
         if (mults.size() > 1) {
           gmm::copy(MM, gmm::sub_matrix
                     (MGLOB, gmm::sub_interval(0,

Modified: trunk/getfem/src/gmm/gmm_range_basis.h
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/gmm/gmm_range_basis.h?rev=4838&r1=4837&r2=4838&view=diff
==============================================================================
--- trunk/getfem/src/gmm/gmm_range_basis.h      (original)
+++ trunk/getfem/src/gmm/gmm_range_basis.h      Sat Dec 20 13:40:22 2014
@@ -437,7 +437,7 @@
           c1.clear();
         }
       }
-      if (c1.size() > 10)
+      if (c1.size() > 1)
         range_basis_eff_Gram_Schmidt_dense(B, c1, c_ortho, EPS);
       for (std::set<size_type>::iterator it = c1.begin(); it != c1.end(); ++it)
         cres.insert(*it);
@@ -446,7 +446,6 @@
       if (columns.size() == nc_r) break;
       if (sizesm[k] >= 350 && columns.size() > (nc_r*19)/20) break;
     }
-    // cout << "size of columns " << columns.size() << endl;
     if (columns.size() > std::max(size_type(10), actsize))
       range_basis_eff_Lanczos(B, columns, EPS);
     else