.

include/sctl/fft_wrapper.hpp

@@ -300,9 +300,15 @@ template <> class FFT<double> : public FFT_Generic<double, FFT<double>> {
 
   void Setup(FFT_Type fft_type_, Long howmany_, const Vector<Long>& dim_vec) {
     if (Dim(0) && Dim(1)) fftw_destroy_plan(plan);
-    Long rank = dim_vec.Dim();
     this->fft_type = fft_type_;
     this->howmany = howmany_;
+
+    Long rank = dim_vec.Dim();
+    Vector<int> dim_vec_(rank);
+    for (Integer i = 0; i < rank; i++) {
+      dim_vec_[i] = dim_vec[i];
+    }
+
     Long N0 = 0, N1 = 0;
     { // Set N0, N1
       Long N = howmany;
@@ -324,21 +330,18 @@ template <> class FFT<double> : public FFT_Generic<double, FFT<double>> {
       this->dim[1] = N1;
     }
     if (!N0 || !N1) return;
-
-    in .ReInit(N0);
-    out.ReInit(N1);
-    Vector<int> dim_vec_(rank);
-    for (Integer i = 0; i < rank; i++) dim_vec_[i] = dim_vec[i];
+    Vector<ValueType> in(N0), out(N1);
 
     if (fft_type == FFT_Type::R2C) {
-      plan = fftw_plan_many_dft_r2c(rank, &dim_vec_[0], howmany_, (double*)&in[0], NULL, 1, N0 / howmany, (fftw_complex*)&out[0], NULL, 1, N1 / 2 / howmany, FFTW_ESTIMATE);
+      plan = fftw_plan_many_dft_r2c(rank, &dim_vec_[0], howmany_, &in[0], NULL, 1, N0 / howmany, (fftw_complex*)&out[0], NULL, 1, N1 / 2 / howmany, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT);
     } else if (fft_type == FFT_Type::C2C) {
-      plan = fftw_plan_many_dft(rank, &dim_vec_[0], howmany_, (fftw_complex*)&in[0], NULL, 1, N0 / 2 / howmany, (fftw_complex*)&out[0], NULL, 1, N1 / 2 / howmany, FFTW_FORWARD, FFTW_ESTIMATE);
+      plan = fftw_plan_many_dft(rank, &dim_vec_[0], howmany_, (fftw_complex*)&in[0], NULL, 1, N0 / 2 / howmany, (fftw_complex*)&out[0], NULL, 1, N1 / 2 / howmany, FFTW_FORWARD, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT);
     } else if (fft_type == FFT_Type::C2C_INV) {
-      plan = fftw_plan_many_dft(rank, &dim_vec_[0], howmany_, (fftw_complex*)&in[0], NULL, 1, N0 / 2 / howmany, (fftw_complex*)&out[0], NULL, 1, N1 / 2 / howmany, FFTW_BACKWARD, FFTW_ESTIMATE);
+      plan = fftw_plan_many_dft(rank, &dim_vec_[0], howmany_, (fftw_complex*)&in[0], NULL, 1, N0 / 2 / howmany, (fftw_complex*)&out[0], NULL, 1, N1 / 2 / howmany, FFTW_BACKWARD, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT);
     } else if (fft_type == FFT_Type::C2R) {
-      plan = fftw_plan_many_dft_c2r(rank, &dim_vec_[0], howmany_, (fftw_complex*)&in[0], NULL, 1, N0 / 2 / howmany, (double*)&out[0], NULL, 1, N1 / howmany, FFTW_ESTIMATE);
+      plan = fftw_plan_many_dft_c2r(rank, &dim_vec_[0], howmany_, (fftw_complex*)&in[0], NULL, 1, N0 / 2 / howmany, &out[0], NULL, 1, N1 / howmany, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT);
     }
+    SCTL_ASSERT(plan);
   }
 
   void Execute(const Vector<ValueType>& in, Vector<ValueType>& out) const {
@@ -367,7 +370,6 @@ template <> class FFT<double> : public FFT_Generic<double, FFT<double>> {
 
  private:
 
-  Vector<ValueType> in, out;
   fftw_plan plan;
 };
 #endif
@@ -383,9 +385,15 @@ template <> class FFT<float> : public FFT_Generic<float, FFT<float>> {
 
   void Setup(FFT_Type fft_type_, Long howmany_, const Vector<Long>& dim_vec) {
     if (Dim(0) && Dim(1)) fftwf_destroy_plan(plan);
-    Long rank = dim_vec.Dim();
     this->fft_type = fft_type_;
     this->howmany = howmany_;
+
+    Long rank = dim_vec.Dim();
+    Vector<int> dim_vec_(rank);
+    for (Integer i = 0; i < rank; i++) {
+      dim_vec_[i] = dim_vec[i];
+    }
+
     Long N0, N1;
     { // Set N0, N1
       Long N = howmany;
@@ -407,21 +415,18 @@ template <> class FFT<float> : public FFT_Generic<float, FFT<float>> {
       this->dim[1] = N1;
     }
     if (!N0 || !N1) return;
-
-    in .ReInit(N0);
-    out.ReInit(N1);
-    Vector<int> dim_vec_(rank);
-    for (Integer i = 0; i < rank; i++) dim_vec_[i] = dim_vec[i];
+    Vector<ValueType> in (N0), out(N1);
 
     if (fft_type == FFT_Type::R2C) {
-      plan = fftwf_plan_many_dft_r2c(rank, &dim_vec_[0], howmany_, (float*)&in[0], NULL, 1, N0 / howmany, (fftwf_complex*)&out[0], NULL, 1, N1 / 2 / howmany, FFTW_ESTIMATE);
+      plan = fftwf_plan_many_dft_r2c(rank, &dim_vec_[0], howmany_, &in[0], NULL, 1, N0 / howmany, (fftwf_complex*)&out[0], NULL, 1, N1 / 2 / howmany, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT);
     } else if (fft_type == FFT_Type::C2C) {
-      plan = fftwf_plan_many_dft(rank, &dim_vec_[0], howmany_, (fftwf_complex*)&in[0], NULL, 1, N0 / 2 / howmany, (fftwf_complex*)&out[0], NULL, 1, N1 / 2 / howmany, FFTW_FORWARD, FFTW_ESTIMATE);
+      plan = fftwf_plan_many_dft(rank, &dim_vec_[0], howmany_, (fftwf_complex*)&in[0], NULL, 1, N0 / 2 / howmany, (fftwf_complex*)&out[0], NULL, 1, N1 / 2 / howmany, FFTW_FORWARD, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT);
     } else if (fft_type == FFT_Type::C2C_INV) {
-      plan = fftwf_plan_many_dft(rank, &dim_vec_[0], howmany_, (fftwf_complex*)&in[0], NULL, 1, N0 / 2 / howmany, (fftwf_complex*)&out[0], NULL, 1, N1 / 2 / howmany, FFTW_BACKWARD, FFTW_ESTIMATE);
+      plan = fftwf_plan_many_dft(rank, &dim_vec_[0], howmany_, (fftwf_complex*)&in[0], NULL, 1, N0 / 2 / howmany, (fftwf_complex*)&out[0], NULL, 1, N1 / 2 / howmany, FFTW_BACKWARD, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT);
     } else if (fft_type == FFT_Type::C2R) {
-      plan = fftwf_plan_many_dft_c2r(rank, &dim_vec_[0], howmany_, (fftwf_complex*)&in[0], NULL, 1, N0 / 2 / howmany, (float*)&out[0], NULL, 1, N1 / howmany, FFTW_ESTIMATE);
+      plan = fftwf_plan_many_dft_c2r(rank, &dim_vec_[0], howmany_, (fftwf_complex*)&in[0], NULL, 1, N0 / 2 / howmany, &out[0], NULL, 1, N1 / howmany, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT);
     }
+    SCTL_ASSERT(plan);
   }
 
   void Execute(const Vector<ValueType>& in, Vector<ValueType>& out) const {
@@ -450,7 +455,6 @@ template <> class FFT<float> : public FFT_Generic<float, FFT<float>> {
 
  private:
 
-  Vector<ValueType> in, out;
   fftwf_plan plan;
 };
 #endif
@@ -466,9 +470,13 @@ template <> class FFT<long double> : public FFT_Generic<long double, FFT<long do
 
   void Setup(FFT_Type fft_type_, Long howmany_, const Vector<Long>& dim_vec) {
     if (Dim(0) && Dim(1)) fftwl_destroy_plan(plan);
-    Long rank = dim_vec.Dim();
     this->fft_type = fft_type_;
     this->howmany = howmany_;
+
+    Long rank = dim_vec.Dim();
+    Vector<int> dim_vec_(rank);
+    for (Integer i = 0; i < rank; i++) dim_vec_[i] = dim_vec[i];
+
     Long N0, N1;
     { // Set N0, N1
       Long N = howmany;
@@ -490,21 +498,18 @@ template <> class FFT<long double> : public FFT_Generic<long double, FFT<long do
       this->dim[1] = N1;
     }
     if (!N0 || !N1) return;
-
-    in .ReInit(N0);
-    out.ReInit(N1);
-    Vector<int> dim_vec_(rank);
-    for (Integer i = 0; i < rank; i++) dim_vec_[i] = dim_vec[i];
+    Vector<ValueType> in (N0), out(N1);
 
     if (fft_type == FFT_Type::R2C) {
-      plan = fftwl_plan_many_dft_r2c(rank, &dim_vec_[0], howmany_, (long double*)&in[0], NULL, 1, N0 / howmany, (fftwl_complex*)&out[0], NULL, 1, N1 / 2 / howmany, FFTW_ESTIMATE);
+      plan = fftwl_plan_many_dft_r2c(rank, &dim_vec_[0], howmany_, &in[0], NULL, 1, N0 / howmany, (fftwl_complex*)&out[0], NULL, 1, N1 / 2 / howmany, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT);
     } else if (fft_type == FFT_Type::C2C) {
-      plan = fftwl_plan_many_dft(rank, &dim_vec_[0], howmany_, (fftwl_complex*)&in[0], NULL, 1, N0 / 2 / howmany, (fftwl_complex*)&out[0], NULL, 1, N1 / 2 / howmany, FFTW_FORWARD, FFTW_ESTIMATE);
+      plan = fftwl_plan_many_dft(rank, &dim_vec_[0], howmany_, (fftwl_complex*)&in[0], NULL, 1, N0 / 2 / howmany, (fftwl_complex*)&out[0], NULL, 1, N1 / 2 / howmany, FFTW_FORWARD, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT);
     } else if (fft_type == FFT_Type::C2C_INV) {
-      plan = fftwl_plan_many_dft(rank, &dim_vec_[0], howmany_, (fftwl_complex*)&in[0], NULL, 1, N0 / 2 / howmany, (fftwl_complex*)&out[0], NULL, 1, N1 / 2 / howmany, FFTW_BACKWARD, FFTW_ESTIMATE);
+      plan = fftwl_plan_many_dft(rank, &dim_vec_[0], howmany_, (fftwl_complex*)&in[0], NULL, 1, N0 / 2 / howmany, (fftwl_complex*)&out[0], NULL, 1, N1 / 2 / howmany, FFTW_BACKWARD, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT);
     } else if (fft_type == FFT_Type::C2R) {
-      plan = fftwl_plan_many_dft_c2r(rank, &dim_vec_[0], howmany_, (fftwl_complex*)&in[0], NULL, 1, N0 / 2 / howmany, (long double*)&out[0], NULL, 1, N1 / howmany, FFTW_ESTIMATE);
+      plan = fftwl_plan_many_dft_c2r(rank, &dim_vec_[0], howmany_, (fftwl_complex*)&in[0], NULL, 1, N0 / 2 / howmany, &out[0], NULL, 1, N1 / howmany, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT);
     }
+    SCTL_ASSERT(plan);
   }
 
   void Execute(const Vector<ValueType>& in, Vector<ValueType>& out) const {
@@ -533,7 +538,6 @@ template <> class FFT<long double> : public FFT_Generic<long double, FFT<long do
 
  private:
 
-  Vector<ValueType> in, out;
   fftwl_plan plan;
 };
 #endif

include/sctl/sph_harm.hpp

@@ -1,7 +1,7 @@
 #ifndef _SCTL_SPH_HARM_HPP_
 #define _SCTL_SPH_HARM_HPP_
 
-#define SCTL_SHMAXDEG 256
+#define SCTL_SHMAXDEG 1024
 
 #include SCTL_INCLUDE(matrix.hpp)
 #include SCTL_INCLUDE(fft_wrapper.hpp)
@@ -9,55 +9,71 @@
 
 namespace SCTL_NAMESPACE {
 
+enum class SHCArrange {
+  // (p+1) x (p+1) complex elements in row-major order.
+  // A : { A(0,0), A(0,1), ... A(0,p), A(1,0), ... A(p,p) }
+  // where, A(n,m) = { Ar(n,m), Ai(n,m) } (real and imaginary parts)
+  ALL,
+
+  // (p+1)(p+2)/2  complex elements in row-major order (lower triangular part)
+  // A : { A(0,0), A(1,0), A(1,1), A(2,0), A(2,1), A(2,2), ... A(p,p) }
+  // where, A(n,m) = { Ar(n,m), Ai(n,m) } (real and imaginary parts)
+  ROW_MAJOR,
+
+  // (p+1)(p+1) real elements in col-major order (non-zero lower triangular part)
+  // A : { Ar(0,0), Ar(1,0), ... Ar(p,0), Ar(1,1), ... Ar(p,1), Ai(1,1), ... Ai(p,1), ..., Ar(p,p), Ai(p,p)
+  // where, A(n,m) = { Ar(n,m), Ai(n,m) } (real and imaginary parts)
+  COL_MAJOR_NONZERO
+};
+
 template <class Real> class SphericalHarmonics{
   static constexpr Integer COORD_DIM = 3;
 
   public:
 
-    // TODO: Ynm *= sqrt(2)*(m==0?1:2);
+    static void Grid2SHC(const Vector<Real>& X_in, Long Nt_in, Long Np_in, Long p_out, Vector<Real>& S_out, SHCArrange arrange_out);
 
-    static void SHC2Grid(const Vector<Real>& S, Long p0, Long p1, Vector<Real>& X, Vector<Real>* X_theta=nullptr, Vector<Real>* X_phi=nullptr);
+    static void SHC2Grid(const Vector<Real>& S_in, SHCArrange arrange_in, Long p_in, Long Nt_out, Long Np_out, Vector<Real>* X_out, Vector<Real>* X_theta_out=nullptr, Vector<Real>* X_phi_out=nullptr);
 
-    static void Grid2SHC(const Vector<Real>& X, Long Nt, Long Np, Long p1, Vector<Real>& S);
-    static void Grid2SHC(const Vector<Real>& X, Long          p0, Long p1, Vector<Real>& S);
+    static void SHC2Pole(const Vector<Real>& S_in, SHCArrange arrange_in, Long p_in, Vector<Real>& P_out);
 
-    static void SHC2GridTranspose(const Vector<Real>& X, Long p0, Long p1, Vector<Real>& S);
+    static void WriteVTK(const char* fname, const Vector<Real>* S, const Vector<Real>* f_val, SHCArrange arrange, Long p_in, Long p_out, Real period=0, const Comm& comm = Comm::World());
 
-    static void SHC2Pole(const Vector<Real>& S, Long p0, Vector<Real>& P);
+    static void test() {
+      int p = 3;
+      int dof = 2;
 
-    static void RotateAll(const Vector<Real>& S, Long p0, Long dof, Vector<Real>& S_);
+      int Ncoeff = (p + 1) * (p + 1);
+      Vector<Real> Xcoeff(dof * Ncoeff);
+      for (int i=0;i<Xcoeff.Dim();i++) Xcoeff[i]=i;
 
-    static void RotateTranspose(const Vector<Real>& S_, Long p0, Long dof, Vector<Real>& S);
+      Vector<Real> Xgrid;
+      int Nt = p+1, Np = 2*p+1;
+      SHC2Grid(Xcoeff, sctl::SHCArrange::COL_MAJOR_NONZERO, p, Nt, Np, &Xgrid);
+      Grid2SHC(Xgrid, Nt, Np, p, Xcoeff, sctl::SHCArrange::ROW_MAJOR);
 
-    static void StokesSingularInteg(const Vector<Real>& S, Long p0, Long p1, Vector<Real>* SLMatrix=nullptr, Vector<Real>* DLMatrix=nullptr);
+      int indx=0;
+      for (int i=0;i<dof;i++) {
+        for (int n=0;n<=p;n++){
+          std::cout<<Vector<Real>(2*n+2,Xcoeff.begin()+indx);
+          indx+=2*n+2;
+        }
+      }
 
-    static void WriteVTK(const char* fname, long p0, long p1, Real period=0, const Vector<Real>* S=nullptr, const Vector<Real>* f_val=nullptr, MPI_Comm comm=MPI_COMM_WORLD);
+      SphericalHarmonics<Real>::WriteVTK("test", nullptr, &Xcoeff, sctl::SHCArrange::ROW_MAJOR, p, 32);
+    }
 
   private:
 
-    static Vector<Real>& LegendreNodes(Long p1);
-
-    static Vector<Real>& LegendreWeights(Long p1);
-
-    static Vector<Real>& SingularWeights(Long p1);
-
-    static Matrix<Real>& MatFourier(Long p0, Long p1);
-
-    static Matrix<Real>& MatFourierInv(Long p0, Long p1);
-
-    static Matrix<Real>& MatFourierGrad(Long p0, Long p1);
-
-    static std::vector<Matrix<Real>>& MatLegendre(Long p0, Long p1);
-
-    static std::vector<Matrix<Real>>& MatLegendreInv(Long p0, Long p1);
-
-    static std::vector<Matrix<Real>>& MatLegendreGrad(Long p0, Long p1);
-
-    static std::vector<Matrix<Real>>& MatRotate(Long p0);
+    // Probably don't work anymore, need to be updated :(
+    static void SHC2GridTranspose(const Vector<Real>& X, Long p0, Long p1, Vector<Real>& S);
+    static void RotateAll(const Vector<Real>& S, Long p0, Long dof, Vector<Real>& S_);
+    static void RotateTranspose(const Vector<Real>& S_, Long p0, Long dof, Vector<Real>& S);
+    static void StokesSingularInteg(const Vector<Real>& S, Long p0, Long p1, Vector<Real>* SLMatrix=nullptr, Vector<Real>* DLMatrix=nullptr);
 
     /**
-     * \brief Computes all the Associated Legendre Polynomials (normalized) upto the specified degree.
-     * \param[in] degree The degree upto which the legendre polynomials have to be computed.
+     * \brief Computes all the Associated Legendre Polynomials (normalized) up to the specified degree.
+     * \param[in] degree The degree up to which the Legendre polynomials have to be computed.
      * \param[in] X The input values for which the polynomials have to be computed.
      * \param[in] N The number of input points.
      * \param[out] poly_val The output array of size (degree+1)*(degree+2)*N/2 containing the computed polynomial values.
@@ -65,16 +81,31 @@ template <class Real> class SphericalHarmonics{
      * P(n,m)[i] => {P(0,0)[0], P(0,0)[1], ..., P(0,0)[N-1], P(1,0)[0], ..., P(1,0)[N-1],
      * P(2,0)[0], ..., P(degree,0)[N-1], P(1,1)[0], ...,P(2,1)[0], ..., P(degree,degree)[N-1]}
      */
-    //static void LegPoly(Vector<Real>& poly_val, const Vector<Real>& X, Long degree);
-    //static void LegPolyDeriv(Vector<Real>& poly_val, const Vector<Real>& X, Long degree);
+    static void LegPoly(Vector<Real>& poly_val, const Vector<Real>& X, Long degree);
+    static void LegPolyDeriv(Vector<Real>& poly_val, const Vector<Real>& X, Long degree);
 
-    static void LegPoly(Real* poly_val, const Real* X, Long N, Long degree);
-    static void LegPolyDeriv(Real* poly_val, const Real* X, Long N, Long degree);
+    static const Vector<Real>& LegendreNodes(Long p1);
+    static const Vector<Real>& LegendreWeights(Long p1);
+    static const Vector<Real>& SingularWeights(Long p1);
+
+    static const Matrix<Real>& MatFourier(Long p0, Long p1);
+    static const Matrix<Real>& MatFourierInv(Long p0, Long p1);
+    static const Matrix<Real>& MatFourierGrad(Long p0, Long p1);
+
+    static const FFT<Real>& OpFourier(Long Np);
+    static const FFT<Real>& OpFourierInv(Long Np);
+
+    static const std::vector<Matrix<Real>>& MatLegendre(Long p0, Long p1);
+    static const std::vector<Matrix<Real>>& MatLegendreInv(Long p0, Long p1);
+    static const std::vector<Matrix<Real>>& MatLegendreGrad(Long p0, Long p1);
+
+    static const std::vector<Matrix<Real>>& MatRotate(Long p0);
 
     template <bool SLayer, bool DLayer> static void StokesSingularInteg_(const Vector<Real>& X0, Long p0, Long p1, Vector<Real>& SL, Vector<Real>& DL);
 
     struct MatrixStorage{
-      MatrixStorage(Long size){
+      MatrixStorage(){
+        const Long size = SCTL_SHMAXDEG;
         Qx_ .resize(size);
         Qw_ .resize(size);
         Sw_ .resize(size);
@@ -96,15 +127,16 @@ template <class Real> class SphericalHarmonics{
       std::vector<std::vector<Matrix<Real>>> Mr_;
       std::vector<Matrix<Real>> Mfinv_ ;
       std::vector<std::vector<Matrix<Real>>> Mlinv_ ;
+
+      StaticArray<FFT<Real>, SCTL_SHMAXDEG> Mfft_ ;
+      StaticArray<FFT<Real>, SCTL_SHMAXDEG> Mfftinv_ ;
     };
     static MatrixStorage& MatrixStore(){
-      static MatrixStorage storage(SCTL_SHMAXDEG);
+      static MatrixStorage storage;
       return storage;
     }
 };
 
-template class SphericalHarmonics<double>;
-
 }  // end namespace
 
 #include SCTL_INCLUDE(sph_harm.txx)

include/sctl/vector.txx

@@ -68,6 +68,10 @@ template <class ValueType> void Vector<ValueType>::Swap(Vector<ValueType>& v1) {
 }
 
 template <class ValueType> void Vector<ValueType>::ReInit(Long dim_, Iterator<ValueType> data_, bool own_data_) {
+#ifdef SCTL_MEMDEBUG
+  Vector<ValueType> tmp(dim_, data_, own_data_);
+  this->Swap(tmp);
+#else
   if (own_data_ && own_data && dim_ <= capacity) {
     dim = dim_;
     if (data_ != NullIterator<ValueType>()) {
@@ -77,6 +81,7 @@ template <class ValueType> void Vector<ValueType>::ReInit(Long dim_, Iterator<Va
     Vector<ValueType> tmp(dim_, data_, own_data_);
     this->Swap(tmp);
   }
+#endif
 }
 
 template <class ValueType> void Vector<ValueType>::Write(const char* fname) const {
@@ -111,7 +116,7 @@ template <class ValueType> void Vector<ValueType>::Read(const char* fname) {
 
 template <class ValueType> inline Long Vector<ValueType>::Dim() const { return dim; }
 
-template <class ValueType> inline Long Vector<ValueType>::Capacity() const { return capacity; }
+//template <class ValueType> inline Long Vector<ValueType>::Capacity() const { return capacity; }
 
 template <class ValueType> void Vector<ValueType>::SetZero() {
   if (dim > 0) memset<ValueType>(data_ptr, 0, dim);
@@ -156,18 +161,14 @@ template <class ValueType> inline const ValueType& Vector<ValueType>::operator[]
 // Vector-Vector operations
 
 template <class ValueType> Vector<ValueType>& Vector<ValueType>::operator=(const std::vector<ValueType>& V) {
-  {
-    if (capacity < V.size()) ReInit(V.size());
-    dim = V.size();
-    memcopy(data_ptr, Ptr2ConstItr<ValueType>(&V[0], V.size()), dim);
-  }
+  if (dim != V.size()) ReInit(V.size());
+  memcopy(data_ptr, Ptr2ConstItr<ValueType>(&V[0], V.size()), dim);
   return *this;
 }
 
 template <class ValueType> Vector<ValueType>& Vector<ValueType>::operator=(const Vector<ValueType>& V) {
   if (this != &V) {
-    if (capacity < V.dim) ReInit(V.dim);
-    dim = V.dim;
+    if (dim != V.dim) ReInit(V.dim);
     memcopy(data_ptr, V.data_ptr, dim);
   }
   return *this;