From: Wolfgang Bangerth <bangerth@colostate.edu>
Date: Fri, 20 Oct 2017 22:32:01 +0000 (-0600)
Subject: Make Vector::values a std::unique_ptr.
X-Git-Tag: v9.0.0-rc1~897^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=7a084ac0ac36f8ba33f111158d7bf36154bd201d;p=dealii.git

Make Vector::values a std::unique_ptr.
---

diff --git a/include/deal.II/lac/vector.h b/include/deal.II/lac/vector.h
index d2de900688..9ee58634f9 100644
--- a/include/deal.II/lac/vector.h
+++ b/include/deal.II/lac/vector.h
@@ -222,7 +222,7 @@ public:
    * Destructor, deallocates memory. Made virtual to allow for derived classes
    * to behave properly.
    */
-  virtual ~Vector ();
+  virtual ~Vector () = default;
 
   /**
    * This function does nothing but exists for compatibility with the parallel
@@ -891,8 +891,12 @@ protected:
 
   /**
    * Pointer to the array of elements of this vector.
+   *
+   * Because we allocate these arrays via Utilities::System::posix_memalign,
+   * we need to use a custom deleter for this object that does not call
+   * <code>delete[]</code>, but instead calls @p free().
    */
-  Number *values;
+  std::unique_ptr<Number[], decltype(&free)> values;
 
   /**
    * For parallel loops with TBB, this member variable stores the affinity
@@ -922,11 +926,6 @@ private:
    * allocated memory is determined by @p max_vec_size .
    */
   void allocate();
-
-  /**
-   * Deallocate @p values.
-   */
-  void deallocate();
 };
 
 /*@}*/
@@ -950,7 +949,7 @@ Vector<Number>::Vector ()
   :
   vec_size(0),
   max_vec_size(0),
-  values(nullptr)
+  values(nullptr, &free)
 {
   reinit(0);
 }
@@ -963,7 +962,7 @@ Vector<Number>::Vector (const InputIterator first, const InputIterator last)
   :
   vec_size (0),
   max_vec_size (0),
-  values (nullptr)
+  values (nullptr, &free)
 {
   // allocate memory. do not initialize it, as we will copy over to it in a
   // second
@@ -979,26 +978,13 @@ Vector<Number>::Vector (const size_type n)
   :
   vec_size(0),
   max_vec_size(0),
-  values(nullptr)
+  values(nullptr, &free)
 {
   reinit (n, false);
 }
 
 
 
-template <typename Number>
-inline
-Vector<Number>::~Vector ()
-{
-  if (values)
-    {
-      deallocate();
-      values = nullptr;
-    }
-}
-
-
-
 template <typename Number>
 inline
 std::size_t Vector<Number>::size () const
@@ -1022,7 +1008,7 @@ inline
 typename Vector<Number>::iterator
 Vector<Number>::begin ()
 {
-  return values;
+  return values.get();
 }
 
 
@@ -1032,7 +1018,7 @@ inline
 typename Vector<Number>::const_iterator
 Vector<Number>::begin () const
 {
-  return values;
+  return values.get();
 }
 
 
@@ -1042,7 +1028,7 @@ inline
 typename Vector<Number>::iterator
 Vector<Number>::end ()
 {
-  return values + vec_size;
+  return values.get() + vec_size;
 }
 
 
@@ -1052,7 +1038,7 @@ inline
 typename Vector<Number>::const_iterator
 Vector<Number>::end () const
 {
-  return values + vec_size;
+  return values.get() + vec_size;
 }
 
 
@@ -1161,7 +1147,7 @@ Vector<Number>::add (const std::vector<size_type> &indices,
 {
   Assert (indices.size() == values.size(),
           ExcDimensionMismatch(indices.size(), values.size()));
-  add (indices.size(), indices.data(), values.values);
+  add (indices.size(), indices.data(), values.values.get());
 }
 
 
@@ -1238,7 +1224,7 @@ Vector<Number>::save (Archive &ar, const unsigned int) const
   ar &static_cast<const Subscriptor &>(*this);
 
   ar &vec_size &max_vec_size ;
-  ar &boost::serialization::make_array(values, max_vec_size);
+  ar &boost::serialization::make_array(values.get(), max_vec_size);
 }
 
 
@@ -1250,14 +1236,14 @@ void
 Vector<Number>::load (Archive &ar, const unsigned int)
 {
   // get rid of previous content
-  deallocate();
+  values.reset();
 
   // the load stuff again from the archive
   ar &static_cast<Subscriptor &>(*this);
   ar &vec_size &max_vec_size ;
 
   allocate();
-  ar &boost::serialization::make_array(values, max_vec_size);
+  ar &boost::serialization::make_array(values.get(), max_vec_size);
 }
 
 #endif
diff --git a/include/deal.II/lac/vector.templates.h b/include/deal.II/lac/vector.templates.h
index a1b19dc83c..c301b955e5 100644
--- a/include/deal.II/lac/vector.templates.h
+++ b/include/deal.II/lac/vector.templates.h
@@ -48,7 +48,7 @@ Vector<Number>::Vector (const Vector<Number> &v)
   Subscriptor(),
   vec_size(v.size()),
   max_vec_size(v.size()),
-  values(nullptr)
+  values(nullptr, &free)
 {
   if (vec_size != 0)
     {
@@ -65,7 +65,7 @@ Vector<Number>::Vector (Vector<Number> &&v)
   Subscriptor(std::move(v)),
   vec_size(v.vec_size),
   max_vec_size(v.max_vec_size),
-  values(v.values),
+  values(std::move(v.values)),
   thread_loop_partitioner(std::move(v.thread_loop_partitioner))
 {
   v.vec_size = 0;
@@ -82,7 +82,7 @@ Vector<Number>::Vector (const Vector<OtherNumber> &v)
   Subscriptor(),
   vec_size(v.size()),
   max_vec_size(v.size()),
-  values(nullptr)
+  values(nullptr, &free)
 {
   if (vec_size != 0)
     {
@@ -145,7 +145,7 @@ Vector<Number>::Vector (const PETScWrappers::VectorBase &v)
   Subscriptor(),
   vec_size(0),
   max_vec_size(0),
-  values(nullptr)
+  values(nullptr, &free)
 {
   if (v.size() != 0)
     {
@@ -163,7 +163,7 @@ Vector<Number>::Vector (const TrilinosWrappers::MPI::Vector &v)
   Subscriptor(),
   vec_size(v.size()),
   max_vec_size(v.size()),
-  values(nullptr)
+  values(nullptr, &free)
 {
   if (vec_size != 0)
     {
@@ -207,7 +207,7 @@ Vector<Number>::operator= (const Vector<Number> &v)
 
   if (vec_size>0)
     {
-      dealii::internal::VectorOperations::Vector_copy<Number,Number> copier(v.values, values);
+      dealii::internal::VectorOperations::Vector_copy<Number,Number> copier(v.values.get(), values.get());
       internal::VectorOperations::parallel_for(copier,0,vec_size,thread_loop_partitioner);
     }
 
@@ -224,16 +224,13 @@ Vector<Number>::operator= (Vector<Number> &&v)
 {
   Subscriptor::operator=(std::move(v));
 
-  if (values) deallocate();
-
   vec_size = v.vec_size;
   max_vec_size = v.max_vec_size;
-  values = v.values;
+  values = std::move(v.values);
   thread_loop_partitioner = std::move(v.thread_loop_partitioner);
 
   v.vec_size = 0;
   v.max_vec_size = 0;
-  v.values = nullptr;
 
   return *this;
 }
@@ -251,7 +248,7 @@ Vector<Number>::operator= (const Vector<Number2> &v)
   if (vec_size != v.vec_size)
     reinit (v, true);
 
-  dealii::internal::VectorOperations::Vector_copy<Number,Number2> copier(v.values, values);
+  dealii::internal::VectorOperations::Vector_copy<Number,Number2> copier(v.values.get(), values.get());
   internal::VectorOperations::parallel_for(copier,0,vec_size,thread_loop_partitioner);
 
   return *this;
@@ -266,19 +263,17 @@ void Vector<Number>::reinit (const size_type n,
 {
   if (n==0)
     {
-      if (values) deallocate();
-      values = nullptr;
+      values.reset();
       max_vec_size = vec_size = 0;
       thread_loop_partitioner.reset(new parallel::internal::TBBPartitioner());
       return;
-    };
+    }
 
   if (n>max_vec_size)
     {
-      if (values) deallocate();
       max_vec_size = n;
       allocate();
-    };
+    }
 
   if (vec_size != n)
     {
@@ -305,18 +300,16 @@ void Vector<Number>::reinit (const Vector<Number2> &v,
 
   if (v.vec_size==0)
     {
-      if (values) deallocate();
-      values = nullptr;
+      values.reset();
       max_vec_size = vec_size = 0;
       return;
-    };
+    }
 
   if (v.vec_size>max_vec_size)
     {
-      if (values) deallocate();
       max_vec_size = v.vec_size;
       allocate();
-    };
+    }
   vec_size = v.vec_size;
   if (omit_zeroing_entries == false)
     *this = Number();
@@ -363,7 +356,7 @@ Vector<Number>::operator= (const Number s)
 
   if (vec_size>0)
     {
-      internal::VectorOperations::Vector_set<Number> setter(s, values);
+      internal::VectorOperations::Vector_set<Number> setter(s, values.get());
       internal::VectorOperations::parallel_for(setter,0,vec_size,thread_loop_partitioner);
     }
 
@@ -379,7 +372,7 @@ Vector<Number> &Vector<Number>::operator *= (const Number factor)
 
   Assert (vec_size!=0, ExcEmptyObject());
 
-  internal::VectorOperations::Vectorization_multiply_factor<Number> vector_multiply(values, factor);
+  internal::VectorOperations::Vectorization_multiply_factor<Number> vector_multiply(values.get(), factor);
 
   internal::VectorOperations::parallel_for(vector_multiply,0,vec_size,thread_loop_partitioner);
 
@@ -398,7 +391,7 @@ Vector<Number>::add (const Number a,
   Assert (vec_size!=0, ExcEmptyObject());
   Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size));
 
-  internal::VectorOperations::Vectorization_add_av<Number> vector_add_av(values, v.values, a);
+  internal::VectorOperations::Vectorization_add_av<Number> vector_add_av(values.get(), v.values.get(), a);
   internal::VectorOperations::parallel_for(vector_add_av,0,vec_size,thread_loop_partitioner);
 }
 
@@ -416,7 +409,7 @@ Vector<Number>::sadd (const Number x,
   Assert (vec_size!=0, ExcEmptyObject());
   Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size));
 
-  internal::VectorOperations::Vectorization_sadd_xav<Number> vector_sadd_xav(values, v.values, a, x);
+  internal::VectorOperations::Vectorization_sadd_xav<Number> vector_sadd_xav(values.get(), v.values.get(), a, x);
   internal::VectorOperations::parallel_for(vector_sadd_xav,0,vec_size,thread_loop_partitioner);
 }
 
@@ -435,7 +428,7 @@ Number Vector<Number>::operator * (const Vector<Number2> &v) const
           ExcDimensionMismatch(vec_size, v.size()));
 
   Number sum;
-  internal::VectorOperations::Dot<Number,Number2> dot(values, v.values);
+  internal::VectorOperations::Dot<Number,Number2> dot(values.get(), v.values.get());
   internal::VectorOperations::parallel_reduce (dot, 0, vec_size, sum, thread_loop_partitioner);
   AssertIsFinite(sum);
 
@@ -451,7 +444,7 @@ Vector<Number>::norm_sqr () const
   Assert (vec_size!=0, ExcEmptyObject());
 
   real_type sum;
-  internal::VectorOperations::Norm2<Number,real_type> norm2(values);
+  internal::VectorOperations::Norm2<Number,real_type> norm2(values.get());
   internal::VectorOperations::parallel_reduce (norm2, 0, vec_size, sum, thread_loop_partitioner);
 
   AssertIsFinite(sum);
@@ -467,7 +460,7 @@ Number Vector<Number>::mean_value () const
   Assert (vec_size!=0, ExcEmptyObject());
 
   Number sum;
-  internal::VectorOperations::MeanValue<Number> mean(values);
+  internal::VectorOperations::MeanValue<Number> mean(values.get());
   internal::VectorOperations::parallel_reduce (mean, 0, vec_size, sum, thread_loop_partitioner);
 
   return sum / real_type(size());
@@ -482,7 +475,7 @@ Vector<Number>::l1_norm () const
   Assert (vec_size!=0, ExcEmptyObject());
 
   real_type sum;
-  internal::VectorOperations::Norm1<Number, real_type> norm1(values);
+  internal::VectorOperations::Norm1<Number, real_type> norm1(values.get());
   internal::VectorOperations::parallel_reduce (norm1, 0, vec_size, sum, thread_loop_partitioner);
 
   return sum;
@@ -502,7 +495,7 @@ Vector<Number>::l2_norm () const
   Assert (vec_size!=0, ExcEmptyObject());
 
   real_type norm_square;
-  internal::VectorOperations::Norm2<Number, real_type> norm2(values);
+  internal::VectorOperations::Norm2<Number, real_type> norm2(values.get());
   internal::VectorOperations::parallel_reduce (norm2, 0, vec_size, norm_square,
                                                thread_loop_partitioner);
   if (numbers::is_finite(norm_square) &&
@@ -546,7 +539,7 @@ Vector<Number>::lp_norm (const real_type p) const
     return l2_norm();
 
   real_type sum;
-  internal::VectorOperations::NormP<Number, real_type> normp(values, p);
+  internal::VectorOperations::NormP<Number, real_type> normp(values.get(), p);
   internal::VectorOperations::parallel_reduce (normp, 0, vec_size, sum, thread_loop_partitioner);
 
   if (numbers::is_finite(sum) && sum >= std::numeric_limits<real_type>::min())
@@ -603,7 +596,7 @@ Vector<Number>::add_and_dot (const Number          a,
   AssertDimension (vec_size, W.size());
 
   Number sum;
-  internal::VectorOperations::AddAndDot<Number> adder(this->values, V.values, W.values, a);
+  internal::VectorOperations::AddAndDot<Number> adder(this->values.get(), V.values.get(), W.values.get(), a);
   internal::VectorOperations::parallel_reduce (adder, 0, vec_size, sum, thread_loop_partitioner);
   AssertIsFinite(sum);
 
@@ -618,7 +611,7 @@ Vector<Number> &Vector<Number>::operator += (const Vector<Number> &v)
   Assert (vec_size!=0, ExcEmptyObject());
   Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size));
 
-  internal::VectorOperations::Vectorization_add_v<Number> vector_add(values, v.values);
+  internal::VectorOperations::Vectorization_add_v<Number> vector_add(values.get(), v.values.get());
   internal::VectorOperations::parallel_for(vector_add,0,vec_size,thread_loop_partitioner);
   return *this;
 }
@@ -631,7 +624,7 @@ Vector<Number> &Vector<Number>::operator -= (const Vector<Number> &v)
   Assert (vec_size!=0, ExcEmptyObject());
   Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size));
 
-  internal::VectorOperations::Vectorization_subtract_v<Number> vector_subtract(values, v.values);
+  internal::VectorOperations::Vectorization_subtract_v<Number> vector_subtract(values.get(), v.values.get());
   internal::VectorOperations::parallel_for(vector_subtract,0,vec_size,thread_loop_partitioner);
 
   return *this;
@@ -644,7 +637,7 @@ void Vector<Number>::add (const Number v)
 {
   Assert (vec_size!=0, ExcEmptyObject());
 
-  internal::VectorOperations::Vectorization_add_factor<Number> vector_add(values, v);
+  internal::VectorOperations::Vectorization_add_factor<Number> vector_add(values.get(), v);
   internal::VectorOperations::parallel_for(vector_add,0,vec_size,thread_loop_partitioner);
 }
 
@@ -661,7 +654,7 @@ void Vector<Number>::add (const Number a, const Vector<Number> &v,
   Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size));
   Assert (vec_size == w.vec_size, ExcDimensionMismatch(vec_size, w.vec_size));
 
-  internal::VectorOperations::Vectorization_add_avpbw<Number> vector_add(values, v.values, w.values, a, b);
+  internal::VectorOperations::Vectorization_add_avpbw<Number> vector_add(values.get(), v.values.get(), w.values.get(), a, b);
   internal::VectorOperations::parallel_for(vector_add,0,vec_size,thread_loop_partitioner);
 }
 
@@ -676,7 +669,7 @@ void Vector<Number>::sadd (const Number x,
   Assert (vec_size!=0, ExcEmptyObject());
   Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size));
 
-  internal::VectorOperations::Vectorization_sadd_xv<Number> vector_sadd(values, v.values, x);
+  internal::VectorOperations::Vectorization_sadd_xv<Number> vector_sadd(values.get(), v.values.get(), x);
   internal::VectorOperations::parallel_for(vector_sadd,0,vec_size,thread_loop_partitioner);
 }
 
@@ -688,7 +681,7 @@ void Vector<Number>::scale (const Vector<Number> &s)
   Assert (vec_size!=0, ExcEmptyObject());
   Assert (vec_size == s.vec_size, ExcDimensionMismatch(vec_size, s.vec_size));
 
-  internal::VectorOperations::Vectorization_scale<Number> vector_scale(values, s.values);
+  internal::VectorOperations::Vectorization_scale<Number> vector_scale(values.get(), s.values.get());
   internal::VectorOperations::parallel_for(vector_scale,0,vec_size,thread_loop_partitioner);
 }
 
@@ -716,7 +709,7 @@ void Vector<Number>::equ (const Number a,
   Assert (vec_size!=0, ExcEmptyObject());
   Assert (vec_size == u.vec_size, ExcDimensionMismatch(vec_size, u.vec_size));
 
-  internal::VectorOperations::Vectorization_equ_au<Number> vector_equ(values, u.values, a);
+  internal::VectorOperations::Vectorization_equ_au<Number> vector_equ(values.get(), u.values.get(), a);
   internal::VectorOperations::parallel_for(vector_equ,0,vec_size,thread_loop_partitioner);
 }
 
@@ -756,7 +749,7 @@ void Vector<Number>::ratio (const Vector<Number> &a,
   // we overwrite them anyway
   reinit (a.size(), true);
 
-  internal::VectorOperations::Vectorization_ratio<Number> vector_ratio(values, a.values, b.values);
+  internal::VectorOperations::Vectorization_ratio<Number> vector_ratio(values.get(), a.values.get(), b.values.get());
   internal::VectorOperations::parallel_for(vector_ratio,0,vec_size,thread_loop_partitioner);
 }
 
@@ -997,23 +990,14 @@ template <typename Number>
 void
 Vector<Number>::allocate()
 {
-  // make sure that we don't create a memory leak
-  Assert (values == nullptr, ExcInternalError());
-
-  // then allocate memory with the proper alignment requirements of 64 bytes
-  Utilities::System::posix_memalign ((void **)&values, 64, sizeof(Number)*max_vec_size);
+  // allocate memory with the proper alignment requirements of 64 bytes
+  Number *new_values;
+  Utilities::System::posix_memalign ((void **)&new_values, 64, sizeof(Number)*max_vec_size);
+  values.reset (new_values);
 }
 
 
 
-template <typename Number>
-void
-Vector<Number>::deallocate()
-{
-  free(values);
-  values = nullptr;
-}
-
 DEAL_II_NAMESPACE_CLOSE
 
 #endif
diff --git a/include/deal.II/lac/vector_view.h b/include/deal.II/lac/vector_view.h
index ed22d53ed3..282db2e6f7 100644
--- a/include/deal.II/lac/vector_view.h
+++ b/include/deal.II/lac/vector_view.h
@@ -238,7 +238,9 @@ VectorView<Number>::VectorView(const size_type new_size, Number *ptr)
 {
   this->vec_size      = new_size;
   this->max_vec_size  = new_size;
-  this->values        = ptr;
+  // release the pointer, but do not delete the object pointed to
+  this->values.release();
+  this->values.reset (ptr);
 }
 
 
@@ -249,7 +251,7 @@ VectorView<Number>::VectorView(const size_type new_size, const Number *ptr)
 {
   this->vec_size      = new_size;
   this->max_vec_size  = new_size;
-  this->values        = const_cast<Number *>(ptr);
+  this->values.reset (const_cast<Number *>(ptr));
 }
 
 
@@ -262,7 +264,9 @@ VectorView<Number>::~VectorView()
   // memory it doesn't own
   this->vec_size = 0;
   this->max_vec_size = 0;
-  this->values = nullptr;
+
+  // release the pointer, but do not delete the object pointed to
+  this->values.release();
 }
 
 
@@ -284,7 +288,9 @@ void VectorView<Number>::reinit(const size_type new_size, Number *ptr)
 {
   this->vec_size      = new_size;
   this->max_vec_size  = new_size;
-  this->values        = ptr;
+  // release the pointer, but do not delete the object pointed to
+  this->values.release();
+  this->values.reset (ptr);
 }
 
 
@@ -294,7 +300,9 @@ void VectorView<Number>::reinit(const size_type new_size, const Number *ptr)
 {
   this->vec_size      = new_size;
   this->max_vec_size  = new_size;
-  this->values        = const_cast<Number *>(ptr);
+  // release the pointer, but do not delete the object pointed to
+  this->values.release();
+  this->values.reset (const_cast<Number *>(ptr));
 }
 
 
diff --git a/source/lac/lapack_full_matrix.cc b/source/lac/lapack_full_matrix.cc
index 570b1a343f..879c906561 100644
--- a/source/lac/lapack_full_matrix.cc
+++ b/source/lac/lapack_full_matrix.cc
@@ -192,7 +192,7 @@ LAPACKFullMatrix<number>::vmult (
       AssertDimension(v.size(), this->n_cols());
       AssertDimension(w.size(), this->n_rows());
 
-      gemv("N", &mm, &nn, &alpha, &this->values[0], &mm, v.values, &one, &beta, w.values, &one);
+      gemv("N", &mm, &nn, &alpha, &this->values[0], &mm, v.values.get(), &one, &beta, w.values.get(), &one);
       break;
     }
     case svd:
@@ -202,12 +202,12 @@ LAPACKFullMatrix<number>::vmult (
       AssertDimension(w.size(), this->n_rows());
       // Compute V^T v
       work.resize(std::max(mm,nn));
-      gemv("N", &nn, &nn, &alpha, &svd_vt->values[0], &nn, v.values, &one, &null, work.data(), &one);
+      gemv("N", &nn, &nn, &alpha, &svd_vt->values[0], &nn, v.values.get(), &one, &null, work.data(), &one);
       // Multiply by singular values
       for (size_type i=0; i<wr.size(); ++i)
         work[i] *= wr[i];
       // Multiply with U
-      gemv("N", &mm, &mm, &alpha, &svd_u->values[0], &mm, work.data(), &one, &beta, w.values, &one);
+      gemv("N", &mm, &mm, &alpha, &svd_u->values[0], &mm, work.data(), &one, &beta, w.values.get(), &one);
       break;
     }
     case inverse_svd:
@@ -217,12 +217,12 @@ LAPACKFullMatrix<number>::vmult (
       AssertDimension(v.size(), this->n_rows());
       // Compute U^T v
       work.resize(std::max(mm,nn));
-      gemv("T", &mm, &mm, &alpha, &svd_u->values[0], &mm, v.values, &one, &null, work.data(), &one);
+      gemv("T", &mm, &mm, &alpha, &svd_u->values[0], &mm, v.values.get(), &one, &null, work.data(), &one);
       // Multiply by singular values
       for (size_type i=0; i<wr.size(); ++i)
         work[i] *= wr[i];
       // Multiply with V
-      gemv("T", &nn, &nn, &alpha, &svd_vt->values[0], &nn, work.data(), &one, &beta, w.values, &one);
+      gemv("T", &nn, &nn, &alpha, &svd_vt->values[0], &nn, work.data(), &one, &beta, w.values.get(), &one);
       break;
     }
     default:
@@ -252,7 +252,7 @@ LAPACKFullMatrix<number>::Tvmult (
       AssertDimension(w.size(), this->n_cols());
       AssertDimension(v.size(), this->n_rows());
 
-      gemv("T", &mm, &nn, &alpha, &this->values[0], &mm, v.values, &one, &beta, w.values, &one);
+      gemv("T", &mm, &nn, &alpha, &this->values[0], &mm, v.values.get(), &one, &beta, w.values.get(), &one);
       break;
     }
     case svd:
@@ -263,12 +263,12 @@ LAPACKFullMatrix<number>::Tvmult (
 
       // Compute U^T v
       work.resize(std::max(mm,nn));
-      gemv("T", &mm, &mm, &alpha, &svd_u->values[0], &mm, v.values, &one, &null, work.data(), &one);
+      gemv("T", &mm, &mm, &alpha, &svd_u->values[0], &mm, v.values.get(), &one, &null, work.data(), &one);
       // Multiply by singular values
       for (size_type i=0; i<wr.size(); ++i)
         work[i] *= wr[i];
       // Multiply with V
-      gemv("T", &nn, &nn, &alpha, &svd_vt->values[0], &nn, work.data(), &one, &beta, w.values, &one);
+      gemv("T", &nn, &nn, &alpha, &svd_vt->values[0], &nn, work.data(), &one, &beta, w.values.get(), &one);
       break;
     }
     case inverse_svd:
@@ -279,12 +279,12 @@ LAPACKFullMatrix<number>::Tvmult (
 
       // Compute V^T v
       work.resize(std::max(mm,nn));
-      gemv("N", &nn, &nn, &alpha, &svd_vt->values[0], &nn, v.values, &one, &null, work.data(), &one);
+      gemv("N", &nn, &nn, &alpha, &svd_vt->values[0], &nn, v.values.get(), &one, &null, work.data(), &one);
       // Multiply by singular values
       for (size_type i=0; i<wr.size(); ++i)
         work[i] *= wr[i];
       // Multiply with U
-      gemv("N", &mm, &mm, &alpha, &svd_u->values[0], &mm, work.data(), &one, &beta, w.values, &one);
+      gemv("N", &mm, &mm, &alpha, &svd_u->values[0], &mm, work.data(), &one, &beta, w.values.get(), &one);
       break;
     }
     default: