[Getfem-commits] r5374 - in /trunk/getfem: contrib/opt_assembly/opt_assembly.cc src/gmm/gmm_blas.h src/gmm/gmm_vector.h

[Yves . Renard]

Author: renard
Date: Mon Sep 26 21:16:41 2016
New Revision: 5374

URL: http://svn.gna.org/viewcvs/getfem?rev=5374&view=rev
Log:
small fixes

Modified:
    trunk/getfem/contrib/opt_assembly/opt_assembly.cc
    trunk/getfem/src/gmm/gmm_blas.h
    trunk/getfem/src/gmm/gmm_vector.h

Modified: trunk/getfem/contrib/opt_assembly/opt_assembly.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/contrib/opt_assembly/opt_assembly.cc?rev=5374&r1=5373&r2=5374&view=diff
==============================================================================
--- trunk/getfem/contrib/opt_assembly/opt_assembly.cc   (original)
+++ trunk/getfem/contrib/opt_assembly/opt_assembly.cc   Mon Sep 26 21:16:41 2016
@@ -424,21 +424,19 @@
   GMM_SET_EXCEPTION_DEBUG; // Exceptions make a memory fault, to debug.
   FE_ENABLE_EXCEPT;        // Enable floating point exception for Nan.
 
-  gmm::row_matrix<gmm::dsvector<double>> MM(15, 15);
+  gmm::row_matrix<gmm::dsvector<double>> MM(15, 200000);
   MM(1, 1) = 5.; MM(1, 1) += 5.;
   MM(1, 3) = 2.;
   MM(2, 2) += 5.; MM(2, 2) -= 4.999;
+  MM(2, 150000) = 6.;
 
   cout << "MM(1, 1) = " << MM(1, 1) << endl;
   cout << "MM(2, 2) = " << MM(2, 2) << endl;
   cout << "nnz(MM) = " << gmm::nnz(MM) << endl;
 
-
-  gmm::resize(MM, 20, 20);
-
   cout << "MM = " << MM << endl;
 
-  gmm::row_matrix<gmm::dsvector<double>> MM2(20, 20);
+  gmm::row_matrix<gmm::dsvector<double>> MM2(15, 200000);
   gmm::copy(MM, MM2);
 
   cout << "MM2 = " << MM2 << endl;

Modified: trunk/getfem/src/gmm/gmm_blas.h
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/gmm/gmm_blas.h?rev=5374&r1=5373&r2=5374&view=diff
==============================================================================
--- trunk/getfem/src/gmm/gmm_blas.h     (original)
+++ trunk/getfem/src/gmm/gmm_blas.h     Mon Sep 26 21:16:41 2016
@@ -691,7 +691,7 @@
     
     std::vector<R> aux(mat_ncols(m));
     for (size_type i = 0; i < mat_nrows(m); ++i) {
-      auto row = mat_const_row(m, i);
+      typename linalg_traits<M>::const_sub_row_type row = mat_const_row(m, i);
       auto it = vect_const_begin(row), ite = vect_const_end(row);
       for (size_type k = 0; it != ite; ++it, ++k)
        aux[index_of_it(it, k, store_type())] += gmm::abs(*it);
@@ -745,7 +745,7 @@
     
     std::vector<R> aux(mat_nrows(m));
     for (size_type i = 0; i < mat_ncols(m); ++i) {
-      auto col = mat_const_col(m, i);
+      typename linalg_traits<M>::const_sub_col_type col = mat_const_col(m, i);
       auto it = vect_const_begin(col), ite = vect_const_end(col);
       for (size_type k = 0; it != ite; ++it, ++k)
        aux[index_of_it(it, k, store_type())] += gmm::abs(*it);
@@ -1156,9 +1156,13 @@
   void copy_vect(const L1& l1, L2& l2, abstract_sparse, abstract_sparse) {
     auto  it  = vect_const_begin(l1), ite = vect_const_end(l1);
     clear(l2);
-    for (; it != ite; ++it)
+    // cout << "copy " << l1 << " of size " << vect_size(l1) << " nnz = " << 
nnz(l1) << endl;
+    for (; it != ite; ++it) {
+      // cout << "*it = " << *it << endl;
+      // cout << "it.index() = " << it.index() << endl;
       if (*it != (typename linalg_traits<L1>::value_type)(0))
        l2[it.index()] = *it;
+    }
   }
   
   template <typename L1, typename L2>
@@ -1950,7 +1954,7 @@
       size_type i,j, k = mat_nrows(l1);
       
       for (i = 0; i < k; ++i) {
-       auto r1=mat_const_row(l1, i);
+       typename linalg_traits<L1>::const_sub_row_type r1=mat_const_row(l1, i);
        for (it2 = it2b, j = 0; it2 != ite; ++it2, ++j)
          l3(i,j) = vect_sp(r1, linalg_traits<L2>::col(it2));
       }
@@ -1982,7 +1986,7 @@
     clear(l3);
     size_type nn = mat_nrows(l3);
     for (size_type i = 0; i < nn; ++i) {
-      auto rl1 = mat_const_row(l1, i);
+      typename linalg_traits<L1>::const_sub_row_type rl1=mat_const_row(l1, i);
       auto it = vect_const_begin(rl1), ite = vect_const_end(rl1);
       for (; it != ite; ++it)
        add(scaled(mat_const_row(l2, it.index()), *it), mat_row(l3, i));
@@ -2024,7 +2028,7 @@
     clear(l3);
     size_type nn = mat_ncols(l3);
     for (size_type i = 0; i < nn; ++i) {
-      auto rc2 = mat_const_col(l2, i);
+      typename linalg_traits<L2>::const_sub_col_type rc2 = mat_const_col(l2, 
i);
       auto it = vect_const_begin(rc2), ite = vect_const_end(rc2);
       for (; it != ite; ++it)
        add(scaled(mat_const_col(l1, it.index()), *it), mat_col(l3, i));
@@ -2075,7 +2079,7 @@
     clear(l3);
     size_type nn = mat_ncols(l1);
     for (size_type i = 0; i < nn; ++i) {
-      auto rc1 = mat_const_col(l1, i);
+      typename linalg_traits<L1>::const_sub_col_type rc1 = mat_const_col(l1, 
i);
       auto it = vect_const_begin(rc1), ite = vect_const_end(rc1);
       for (; it != ite; ++it)
        add(scaled(mat_const_row(l2, i), *it), mat_row(l3, it.index()));
@@ -2127,7 +2131,7 @@
   bool is_symmetric(const MAT &A, magnitude_of_linalg(MAT) tol, 
                    row_major) {
     for (size_type i = 0; i < mat_nrows(A); ++i) {
-      auto row = mat_const_row(A, i);
+      typename linalg_traits<MAT>::const_sub_row_type row = mat_const_row(A, 
i);
       auto it = vect_const_begin(row), ite = vect_const_end(row);
       for (; it != ite; ++it)
        if (gmm::abs(*it - A(it.index(), i)) > tol) return false;
@@ -2139,7 +2143,7 @@
   bool is_symmetric(const MAT &A, magnitude_of_linalg(MAT) tol, 
                    col_major) {
     for (size_type i = 0; i < mat_ncols(A); ++i) {
-      auto col = mat_const_col(A, i);
+      typename linalg_traits<MAT>::const_sub_col_type col = mat_const_col(A, 
i);
       auto it = vect_const_begin(col), ite = vect_const_end(col);
       for (; it != ite; ++it)
        if (gmm::abs(*it - A(i, it.index())) > tol) return false;
@@ -2188,7 +2192,7 @@
   bool is_hermitian(const MAT &A, magnitude_of_linalg(MAT) tol, 
                    row_major) {
     for (size_type i = 0; i < mat_nrows(A); ++i) {
-      auto row = mat_const_row(A, i);
+      typename linalg_traits<MAT>::const_sub_row_type row = mat_const_row(A, 
i);
       auto it = vect_const_begin(row), ite = vect_const_end(row);
       for (; it != ite; ++it)
        if (gmm::abs(gmm::conj(*it) - A(it.index(), i)) > tol) return false;
@@ -2200,7 +2204,7 @@
   bool is_hermitian(const MAT &A, magnitude_of_linalg(MAT) tol, 
                    col_major) {
     for (size_type i = 0; i < mat_ncols(A); ++i) {
-      auto col = mat_const_col(A, i);
+      typename linalg_traits<MAT>::const_sub_col_type col = mat_const_col(A, 
i);
       auto it = vect_const_begin(col), ite = vect_const_end(col);
       for (; it != ite; ++it)
        if (gmm::abs(gmm::conj(*it) - A(i, it.index())) > tol) return false;

Modified: trunk/getfem/src/gmm/gmm_vector.h
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/gmm/gmm_vector.h?rev=5374&r1=5373&r2=5374&view=diff
==============================================================================
--- trunk/getfem/src/gmm/gmm_vector.h   (original)
+++ trunk/getfem/src/gmm/gmm_vector.h   Mon Sep 26 21:16:41 2016
@@ -365,7 +365,7 @@
        size_type j = (i & my_mask) >> my_shift;
        void_pointer q = ((void_pointer *)(p))[j];
        if (!q) {
-         if (k+1 == depth) {
+         if (k+1 != depth) {
            q = new void_pointer[16];
            std::memset(q, 0, 16*sizeof(void_pointer));
          } else {
@@ -417,7 +417,7 @@
       if (my_depth) {
        for (size_type k = 0; k < 16; ++k)
          if (((void_pointer *)(p))[k] && (base + (k+1)*(mask+1)) >= i)
-           rec_clean_i(((void_pointer *)(p))[k], my_depth-1, (my_mask << 4),
+           rec_clean_i(((void_pointer *)(p))[k], my_depth-1, (my_mask >> 4),
                        i, base + k*(mask+1));
       } else {
        for (size_type k = 0; k < 16; ++k)
@@ -441,8 +441,8 @@
 
     void copy_rec(void_pointer &p, const_void_pointer q, size_type my_depth) {
       if (depth) {
-       p = new void_pointer [16];
-       std::memset(root_ptr, 0, 16*sizeof(void_pointer));
+       p = new void_pointer[16];
+       std::memset(p, 0, 16*sizeof(void_pointer));
        for (size_type l = 0; l < 16; ++l)
          if (((const const_void_pointer *)(q))[l])
            copy_rec(((void_pointer *)(p))[l],
@@ -461,12 +461,15 @@
 
     void next_pos_rec(void_pointer p, size_type my_depth, size_type my_mask,
                      const_pointer &pp, size_type &i, size_type base) const {
+      // cout << "base = " << base << endl;
+      // cout << "mask+1 = " << my_mask+1 << " (mask >> 4) = " << (my_mask >> 
4) << endl;
       size_type ii = i;
       if (my_depth) {
+       my_mask = (my_mask >> 4);
        for (size_type k = 0; k < 16; ++k)
-         if (((void_pointer *)(p))[k] && (base + (k+1)*(mask+1)) >= i) {
-           next_pos_rec(((void_pointer *)(p))[k], my_depth-1, (my_mask << 4),
-                        pp, i, base + k*(mask+1));
+         if (((void_pointer *)(p))[k] && (base + (k+1)*(my_mask+1)) >= i) {
+           next_pos_rec(((void_pointer *)(p))[k], my_depth-1, my_mask,
+                        pp, i, base + k*(my_mask+1));
            if (i != size_type(-1)) return; else i = ii;
        }
        i = size_type(-1); pp = 0;
@@ -483,10 +486,11 @@
                          size_type base) const {
       size_type ii = i;
       if (my_depth) {
+       my_mask = (my_mask >> 4);
        for (size_type k = 15; k != size_type(-1); --k)
-         if (((void_pointer *)(p))[k] && ((base + k*(mask+1)) < i)) {
-           next_pos_rec(((void_pointer *)(p))[k], my_depth-1, (my_mask << 4),
-                        pp, i, base + k*(mask+1));
+         if (((void_pointer *)(p))[k] && ((base + k*(my_mask+1)) < i)) {
+           previous_pos_rec(((void_pointer *)(p))[k], my_depth-1,
+                            my_mask, pp, i, base + k*(my_mask+1));
            if (i != size_type(-1)) return; else i = ii;
        }
        i = size_type(-1); pp = 0;