diff --git a/examples/r_KLU_GLU.cpp b/examples/r_KLU_GLU.cpp
index 82c949db..c072eb34 100644
--- a/examples/r_KLU_GLU.cpp
+++ b/examples/r_KLU_GLU.cpp
@@ -145,7 +145,7 @@ int main(int argc, char *argv[])
     vec_r->update(rhs, "cpu", "cuda");
 
 
-    matrix_handler->setValuesChanged(true);
+    matrix_handler->setValuesChanged(true, "cuda");
     matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr", "cuda"); 
 
     printf("\t 2-Norm of the residual: %16.16e\n", sqrt(vector_handler->dot(vec_r, vec_r, "cuda")));
diff --git a/examples/r_KLU_GLU_matrix_values_update.cpp b/examples/r_KLU_GLU_matrix_values_update.cpp
index 19d7ec3c..af797883 100644
--- a/examples/r_KLU_GLU_matrix_values_update.cpp
+++ b/examples/r_KLU_GLU_matrix_values_update.cpp
@@ -155,7 +155,7 @@ int main(int argc, char *argv[])
       vec_r->update(rhs, "cpu", "cuda");
 
 
-      matrix_handler->setValuesChanged(true);
+      matrix_handler->setValuesChanged(true, "cuda");
       matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr", "cuda"); 
 
       printf("\t 2-Norm of the residual: %16.16e\n", sqrt(vector_handler->dot(vec_r, vec_r, "cuda")));
diff --git a/examples/r_KLU_KLU.cpp b/examples/r_KLU_KLU.cpp
index 85c84921..e777833a 100644
--- a/examples/r_KLU_KLU.cpp
+++ b/examples/r_KLU_KLU.cpp
@@ -134,7 +134,7 @@ int main(int argc, char *argv[])
     }
     vec_r->update(rhs, "cpu", "cpu");
 
-    matrix_handler->setValuesChanged(true);
+    matrix_handler->setValuesChanged(true, "cpu");
 
     matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE, "csr", "cpu"); 
     real_type* test = vec_r->getData("cpu");
diff --git a/examples/r_KLU_KLU_standalone.cpp b/examples/r_KLU_KLU_standalone.cpp
index f0b7c3a4..daeae43a 100644
--- a/examples/r_KLU_KLU_standalone.cpp
+++ b/examples/r_KLU_KLU_standalone.cpp
@@ -96,7 +96,7 @@ int main(int argc, char *argv[])
   std::cout << "KLU solve status: " << status << std::endl;      
   vec_r->update(rhs, "cpu", "cpu");
 
-  matrix_handler->setValuesChanged(true);
+  matrix_handler->setValuesChanged(true, "cpu");
 
   matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE, "csr", "cpu"); 
   real_type* test = vec_r->getData("cpu");
diff --git a/examples/r_KLU_rf.cpp b/examples/r_KLU_rf.cpp
index f2e96dfe..8954724d 100644
--- a/examples/r_KLU_rf.cpp
+++ b/examples/r_KLU_rf.cpp
@@ -158,7 +158,7 @@ int main(int argc, char *argv[] )
     }
     vec_r->update(rhs, "cpu", "cuda");
 
-    matrix_handler->setValuesChanged(true);
+    matrix_handler->setValuesChanged(true, "cuda");
 
     matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr", "cuda"); 
 
diff --git a/examples/r_KLU_rf_FGMRES.cpp b/examples/r_KLU_rf_FGMRES.cpp
index d9543895..52f54a0d 100644
--- a/examples/r_KLU_rf_FGMRES.cpp
+++ b/examples/r_KLU_rf_FGMRES.cpp
@@ -125,7 +125,7 @@ int main(int argc, char *argv[])
     real_type norm_b;
     if (i < 2){
       KLU->setup(A);
-      matrix_handler->setValuesChanged(true);
+      matrix_handler->setValuesChanged(true, "cuda");
       status = KLU->analyze();
       std::cout<<"KLU analysis status: "<<status<<std::endl;
       status = KLU->factorize();
@@ -135,7 +135,7 @@ int main(int argc, char *argv[])
       vec_r->update(rhs, "cpu", "cuda");
       norm_b = vector_handler->dot(vec_r, vec_r, "cuda");
       norm_b = sqrt(norm_b);
-      matrix_handler->setValuesChanged(true);
+      matrix_handler->setValuesChanged(true, "cuda");
       matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr", "cuda"); 
       printf("\t 2-Norm of the residual : %16.16e\n", sqrt(vector_handler->dot(vec_r, vec_r, "cuda"))/norm_b);
       if (i == 1) {
@@ -165,7 +165,7 @@ int main(int argc, char *argv[])
        norm_b = vector_handler->dot(vec_r, vec_r, "cuda");
       norm_b = sqrt(norm_b);
 
-      //matrix_handler->setValuesChanged(true);
+      //matrix_handler->setValuesChanged(true, "cuda");
       FGMRES->resetMatrix(A);
       FGMRES->setupPreconditioner("CuSolverRf", Rf);
       
diff --git a/examples/r_KLU_rf_FGMRES_reuse_factorization.cpp b/examples/r_KLU_rf_FGMRES_reuse_factorization.cpp
index d50bdb4a..2c484926 100644
--- a/examples/r_KLU_rf_FGMRES_reuse_factorization.cpp
+++ b/examples/r_KLU_rf_FGMRES_reuse_factorization.cpp
@@ -126,7 +126,7 @@ int main(int argc, char *argv[])
     real_type norm_b;
     if (i < 2){
       KLU->setup(A);
-      matrix_handler->setValuesChanged(true);
+      matrix_handler->setValuesChanged(true, "cuda");
       status = KLU->analyze();
       std::cout<<"KLU analysis status: "<<status<<std::endl;
       status = KLU->factorize();
@@ -136,7 +136,7 @@ int main(int argc, char *argv[])
       vec_r->update(rhs, "cpu", "cuda");
       norm_b = vector_handler->dot(vec_r, vec_r, "cuda");
       norm_b = sqrt(norm_b);
-      matrix_handler->setValuesChanged(true);
+      matrix_handler->setValuesChanged(true, "cuda");
       matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr", "cuda"); 
       printf("\t 2-Norm of the residual : %16.16e\n", sqrt(vector_handler->dot(vec_r, vec_r, "cuda"))/norm_b);
       if (i == 1) {
@@ -178,7 +178,7 @@ int main(int argc, char *argv[])
        norm_b = vector_handler->dot(vec_r, vec_r, "cuda");
       norm_b = sqrt(norm_b);
 
-      matrix_handler->setValuesChanged(true);
+      matrix_handler->setValuesChanged(true, "cuda");
       FGMRES->resetMatrix(A);
       
       matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr", "cuda"); 
diff --git a/resolve/LinSolverIterativeFGMRES.cpp b/resolve/LinSolverIterativeFGMRES.cpp
index 2a6da732..bd7dca42 100644
--- a/resolve/LinSolverIterativeFGMRES.cpp
+++ b/resolve/LinSolverIterativeFGMRES.cpp
@@ -308,7 +308,7 @@ namespace ReSolve
   int  LinSolverIterativeFGMRES::resetMatrix(matrix::Sparse* new_matrix)
   {
     A_ = new_matrix;
-    matrix_handler_->setValuesChanged(true);
+    matrix_handler_->setValuesChanged(true, "cuda");
     return 0;
   }
 
diff --git a/resolve/matrix/MatrixHandler.cpp b/resolve/matrix/MatrixHandler.cpp
index 7005046d..6ebadfa7 100644
--- a/resolve/matrix/MatrixHandler.cpp
+++ b/resolve/matrix/MatrixHandler.cpp
@@ -50,7 +50,6 @@ namespace ReSolve {
   MatrixHandler::MatrixHandler()
   {
     this->new_matrix_ = true;
-    this->values_changed_ = true;
     cpuImpl_  = new MatrixHandlerCpu();
     cudaImpl_ = new MatrixHandlerCuda();
   }
@@ -66,10 +65,15 @@ namespace ReSolve {
     cudaImpl_ = new MatrixHandlerCuda(new_workspace);
   }
 
-  void MatrixHandler::setValuesChanged(bool toWhat)
+  void MatrixHandler::setValuesChanged(bool isValuesChanged, std::string memspace)
   {
-    this->values_changed_ = toWhat;
-    cpuImpl_->setValuesChanged(values_changed_);
+    if (memspace == "cpu") {
+      cpuImpl_->setValuesChanged(isValuesChanged);
+    } else if (memspace == "cuda") {
+      cudaImpl_->setValuesChanged(isValuesChanged);
+    } else {
+      out::error() << "Unsupported device " << memspace << "\n";
+    }
   }
 
   int MatrixHandler::coo2csr(matrix::Coo* A_coo, matrix::Csr* A_csr, std::string memspace)
diff --git a/resolve/matrix/MatrixHandler.hpp b/resolve/matrix/MatrixHandler.hpp
index af444abd..59ace249 100644
--- a/resolve/matrix/MatrixHandler.hpp
+++ b/resolve/matrix/MatrixHandler.hpp
@@ -71,15 +71,14 @@ namespace ReSolve {
                  std::string matrix_type,
                  std::string memspace);
       int Matrix1Norm(matrix::Sparse *A, real_type* norm);
-      void setValuesChanged(bool toWhat); 
+      void setValuesChanged(bool toWhat, std::string memspace); 
     
     private: 
       LinAlgWorkspace* workspace_{nullptr};
       bool new_matrix_{true};     ///< if the structure changed, you need a new handler.
-      bool values_changed_{true}; ///< needed for matvec
 
       MemoryHandler mem_; ///< Device memory manager object
-      MatrixHandlerImpl* cpuImpl_{nullptr};
+      MatrixHandlerImpl*  cpuImpl_{nullptr};
       MatrixHandlerImpl* cudaImpl_{nullptr};
   };
 
diff --git a/resolve/matrix/MatrixHandlerCpu.cpp b/resolve/matrix/MatrixHandlerCpu.cpp
index d06a4ce5..522d2a1b 100644
--- a/resolve/matrix/MatrixHandlerCpu.cpp
+++ b/resolve/matrix/MatrixHandlerCpu.cpp
@@ -65,10 +65,10 @@ namespace ReSolve {
   //   return this->values_changed_;
   // }
 
-  // void MatrixHandlerCpu::setValuesChanged(bool toWhat)
-  // {
-  //   this->values_changed_ = toWhat;
-  // }
+  void MatrixHandlerCpu::setValuesChanged(bool values_changed)
+  {
+    values_changed_ = values_changed;
+  }
 
 //   int MatrixHandlerCpu::coo2csr(matrix::Coo* A_coo, matrix::Csr* A_csr, std::string memspace)
 //   {
diff --git a/resolve/matrix/MatrixHandlerCpu.hpp b/resolve/matrix/MatrixHandlerCpu.hpp
index 60079a89..fb760750 100644
--- a/resolve/matrix/MatrixHandlerCpu.hpp
+++ b/resolve/matrix/MatrixHandlerCpu.hpp
@@ -51,12 +51,12 @@ namespace ReSolve {
                  const real_type* beta,
                  std::string matrix_type);
       virtual int Matrix1Norm(matrix::Sparse *A, real_type* norm);
-      // void setValuesChanged(bool toWhat); 
+      void setValuesChanged(bool isValuesChanged); 
     
     private: 
       LinAlgWorkspace* workspace_{nullptr};
       // bool new_matrix_{true};     ///< if the structure changed, you need a new handler.
-      // bool values_changed_{true}; ///< needed for matvec
+      bool values_changed_{true}; ///< needed for matvec
 
       MemoryHandler mem_; ///< Device memory manager object
   };
diff --git a/resolve/matrix/MatrixHandlerCuda.cpp b/resolve/matrix/MatrixHandlerCuda.cpp
index 31515398..f40de0d8 100644
--- a/resolve/matrix/MatrixHandlerCuda.cpp
+++ b/resolve/matrix/MatrixHandlerCuda.cpp
@@ -65,10 +65,10 @@ namespace ReSolve {
   //   return this->values_changed_;
   // }
 
-  // void MatrixHandlerCuda::setValuesChanged(bool toWhat)
-  // {
-  //   this->values_changed_ = toWhat;
-  // }
+  void MatrixHandlerCuda::setValuesChanged(bool values_changed)
+  {
+    values_changed_ = values_changed;
+  }
 
 //   int MatrixHandlerCuda::coo2csr(matrix::Coo* A_coo, matrix::Csr* A_csr, std::string memspace)
 //   {
diff --git a/resolve/matrix/MatrixHandlerCuda.hpp b/resolve/matrix/MatrixHandlerCuda.hpp
index ae0c3709..f2250008 100644
--- a/resolve/matrix/MatrixHandlerCuda.hpp
+++ b/resolve/matrix/MatrixHandlerCuda.hpp
@@ -51,12 +51,12 @@ namespace ReSolve {
                  const real_type* beta,
                  std::string matrix_type);
       virtual int Matrix1Norm(matrix::Sparse *A, real_type* norm);
-      // void setValuesChanged(bool toWhat); 
+      void setValuesChanged(bool isValuesChanged); 
     
     private: 
       LinAlgWorkspace* workspace_{nullptr};
       // bool new_matrix_{true};     ///< if the structure changed, you need a new handler.
-      // bool values_changed_{true}; ///< needed for matvec
+      bool values_changed_{true}; ///< needed for matvec
 
       MemoryHandler mem_; ///< Device memory manager object
   };
diff --git a/resolve/matrix/MatrixHandlerImpl.hpp b/resolve/matrix/MatrixHandlerImpl.hpp
index ba24d498..015412d0 100644
--- a/resolve/matrix/MatrixHandlerImpl.hpp
+++ b/resolve/matrix/MatrixHandlerImpl.hpp
@@ -46,27 +46,19 @@ namespace ReSolve {
 
       /// Should compute vec_result := alpha*A*vec_x + beta*vec_result, but at least on cpu alpha and beta are flipped
       virtual int matvec(matrix::Sparse* A,
-                 vector_type* vec_x,
-                 vector_type* vec_result,
-                 const real_type* alpha,
-                 const real_type* beta,
-                 std::string matrix_type) = 0;
+                         vector_type* vec_x,
+                         vector_type* vec_result,
+                         const real_type* alpha,
+                         const real_type* beta,
+                         std::string matrix_type) = 0;
       virtual int Matrix1Norm(matrix::Sparse* A, real_type* norm) = 0;
 
-      bool isValuesChanged()
-      {
-        return values_changed_;
-      }
-
-      void setValuesChanged(bool isValuesChanged)
-      {
-        values_changed_ = isValuesChanged;
-      } 
+      virtual void setValuesChanged(bool isValuesChanged) = 0;
     
     protected: 
       // LinAlgWorkspace* workspace_{nullptr};
       // bool new_matrix_{true};     ///< if the structure changed, you need a new handler.
-      bool values_changed_{true}; ///< needed for matvec
+      // bool values_changed_{true}; ///< needed for matvec
 
       // MemoryHandler mem_; ///< Device memory manager object
   };
diff --git a/tests/functionality/testKLU.cpp b/tests/functionality/testKLU.cpp
index a60083cd..70ee83b0 100644
--- a/tests/functionality/testKLU.cpp
+++ b/tests/functionality/testKLU.cpp
@@ -104,7 +104,7 @@ int main(int argc, char *argv[])
   vec_diff->update(x_data, "cpu", "cpu");
 
   // real_type normXmatrix1 = sqrt(vector_handler->dot(vec_test, vec_test, "cpu"));
-  matrix_handler->setValuesChanged(true);
+  matrix_handler->setValuesChanged(true, "cpu");
   status = matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr","cpu"); 
   error_sum += status;
   
@@ -174,7 +174,7 @@ int main(int argc, char *argv[])
   error_sum += status;
 
   vec_r->update(rhs, "cpu", "cpu");
-  matrix_handler->setValuesChanged(true);
+  matrix_handler->setValuesChanged(true, "cpu");
 
   status = matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE, "csr", "cpu"); 
   error_sum += status;
diff --git a/tests/functionality/testKLU_GLU.cpp b/tests/functionality/testKLU_GLU.cpp
index 60dcce50..1c66fda1 100644
--- a/tests/functionality/testKLU_GLU.cpp
+++ b/tests/functionality/testKLU_GLU.cpp
@@ -127,7 +127,7 @@ int main(int argc, char *argv[])
   vec_diff->update(x_data, "cpu", "cuda");
 
   // real_type normXmatrix1 = sqrt(vector_handler->dot(vec_test, vec_test, "cuda"));
-  matrix_handler->setValuesChanged(true);
+  matrix_handler->setValuesChanged(true, "cuda");
   status = matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr","cuda"); 
   error_sum += status;
   
@@ -198,7 +198,7 @@ int main(int argc, char *argv[])
   error_sum += status;
 
    vec_r->update(rhs, "cpu", "cuda");
-   matrix_handler->setValuesChanged(true);
+   matrix_handler->setValuesChanged(true, "cuda");
 
   status = matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE, "csr", "cuda"); 
   error_sum += status;
diff --git a/tests/functionality/testKLU_Rf.cpp b/tests/functionality/testKLU_Rf.cpp
index 50add7cd..656aa724 100644
--- a/tests/functionality/testKLU_Rf.cpp
+++ b/tests/functionality/testKLU_Rf.cpp
@@ -111,7 +111,7 @@ int main(int argc, char *argv[])
   vec_diff->update(x_data, "cpu", "cuda");
 
   // real_type normXmatrix1 = sqrt(vector_handler->dot(vec_test, vec_test, "cuda"));
-  matrix_handler->setValuesChanged(true);
+  matrix_handler->setValuesChanged(true, "cuda");
   status = matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr","cuda"); 
   error_sum += status;
   
@@ -194,8 +194,8 @@ int main(int argc, char *argv[])
   status = Rf->solve(vec_rhs, vec_x);
   error_sum += status;
 
-   vec_r->update(rhs, "cpu", "cuda");
-  matrix_handler->setValuesChanged(true);
+  vec_r->update(rhs, "cpu", "cuda");
+  matrix_handler->setValuesChanged(true, "cuda");
 
   status = matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE, "csr", "cuda"); 
   error_sum += status;
diff --git a/tests/functionality/testKLU_Rf_FGMRES.cpp b/tests/functionality/testKLU_Rf_FGMRES.cpp
index 2dae74b5..14bb0f33 100644
--- a/tests/functionality/testKLU_Rf_FGMRES.cpp
+++ b/tests/functionality/testKLU_Rf_FGMRES.cpp
@@ -117,7 +117,7 @@ int main(int argc, char *argv[])
   vec_diff->update(x_data, "cpu", "cuda");
 
   // real_type normXmatrix1 = sqrt(vector_handler->dot(vec_test, vec_test, "cuda"));
-  matrix_handler->setValuesChanged(true);
+  matrix_handler->setValuesChanged(true, "cuda");
   //evaluate the residual ||b-Ax||
   status = matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr","cuda"); 
   error_sum += status;
@@ -221,7 +221,7 @@ int main(int argc, char *argv[])
   error_sum += status;
 
   vec_r->update(rhs, "cpu", "cuda");
-  matrix_handler->setValuesChanged(true);
+  matrix_handler->setValuesChanged(true, "cuda");
 
   //evaluate final residual
   status = matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE, "csr", "cuda"); 
diff --git a/tests/unit/matrix/MatrixHandlerTests.hpp b/tests/unit/matrix/MatrixHandlerTests.hpp
index 41df9f4b..16fd3215 100644
--- a/tests/unit/matrix/MatrixHandlerTests.hpp
+++ b/tests/unit/matrix/MatrixHandlerTests.hpp
@@ -57,7 +57,7 @@ class MatrixHandlerTests : TestBase
 
     real_type alpha = 2.0/30.0;
     real_type beta  = 2.0;
-    handler.setValuesChanged(true);
+    handler.setValuesChanged(true, memspace_);
     handler.matvec(A, &x, &y, &alpha, &beta, "csr", memspace_);
 
     status *= verifyAnswer(y, 4.0, memspace_);