Changes in src/vector.cpp [112b09:9d5ddf]

File:

• src/vector.cpp (modified) (19 diffs)

src/vector.cpp

@@ -8 +8 @@
 
 #include "vector.hpp"
+#include "VectorContent.hpp"
 #include "verbose.hpp"
 #include "World.hpp"
 #include "Helpers/Assert.hpp"
 #include "Helpers/fast_functions.hpp"
+#include "Exceptions/MathException.hpp"
 
 #include <iostream>
+#include <gsl/gsl_blas.h>
+
 
 using namespace std;
@@ -24 +28 @@
 Vector::Vector()
 {
-  x[0] = x[1] = x[2] = 0.;
+  content = new VectorContent();
 };
 
@@ -33 +37 @@
 Vector::Vector(const Vector& src)
 {
-  x[0] = src[0];
-  x[1] = src[1];
-  x[2] = src[2];
+  content = new VectorContent();
+  gsl_vector_memcpy(content->content, src.content->content);
 }
 
@@ -42 +45 @@
 Vector::Vector(const double x1, const double x2, const double x3)
 {
-  x[0] = x1;
-  x[1] = x2;
-  x[2] = x3;
-};
+  content = new VectorContent();
+  gsl_vector_set(content->content,0,x1);
+  gsl_vector_set(content->content,1,x2);
+  gsl_vector_set(content->content,2,x3);
+};
+
+Vector::Vector(VectorContent *_content) :
+  content(_content)
+{}
 
 /**
@@ -53 +61 @@
   // check for self assignment
   if(&src!=this){
-    x[0] = src[0];
-    x[1] = src[1];
-    x[2] = src[2];
+    gsl_vector_memcpy(content->content, src.content->content);
   }
   return *this;
@@ -62 +68 @@
 /** Desctructor of class vector.
  */
-Vector::~Vector() {};
+Vector::~Vector() {
+  delete content;
+};
 
 /** Calculates square of distance between this and another vector.
@@ -72 +80 @@
   double res = 0.;
   for (int i=NDIM;i--;)
-    res += (x[i]-y[i])*(x[i]-y[i]);
+    res += (at(i)-y[i])*(at(i)-y[i]);
   return (res);
 };
@@ -90 +98 @@
 }
 
-/** Calculates distance between this and another vector in a periodic cell.
- * \param *y array to second vector
- * \param *cell_size 6-dimensional array with (xx, xy, yy, xz, yz, zz) entries specifying the periodic cell
- * \return \f$| x - y |\f$
- */
-double Vector::PeriodicDistance(const Vector &y, const double * const cell_size) const
-{
-  double res = distance(y), tmp, matrix[NDIM*NDIM];
-    Vector Shiftedy, TranslationVector;
-    int N[NDIM];
-    matrix[0] = cell_size[0];
-    matrix[1] = cell_size[1];
-    matrix[2] = cell_size[3];
-    matrix[3] = cell_size[1];
-    matrix[4] = cell_size[2];
-    matrix[5] = cell_size[4];
-    matrix[6] = cell_size[3];
-    matrix[7] = cell_size[4];
-    matrix[8] = cell_size[5];
-    // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
-    for (N[0]=-1;N[0]<=1;N[0]++)
-      for (N[1]=-1;N[1]<=1;N[1]++)
-        for (N[2]=-1;N[2]<=1;N[2]++) {
-          // create the translation vector
-          TranslationVector.Zero();
-          for (int i=NDIM;i--;)
-            TranslationVector[i] = (double)N[i];
-          TranslationVector.MatrixMultiplication(matrix);
-          // add onto the original vector to compare with
-          Shiftedy = y + TranslationVector;
-          // get distance and compare with minimum so far
-          tmp = distance(Shiftedy);
-          if (tmp < res) res = tmp;
-        }
-    return (res);
-};
-
-/** Calculates distance between this and another vector in a periodic cell.
- * \param *y array to second vector
- * \param *cell_size 6-dimensional array with (xx, xy, yy, xz, yz, zz) entries specifying the periodic cell
- * \return \f$| x - y |^2\f$
- */
-double Vector::PeriodicDistanceSquared(const Vector &y, const double * const cell_size) const
-{
-  double res = DistanceSquared(y), tmp, matrix[NDIM*NDIM];
-    Vector Shiftedy, TranslationVector;
-    int N[NDIM];
-    matrix[0] = cell_size[0];
-    matrix[1] = cell_size[1];
-    matrix[2] = cell_size[3];
-    matrix[3] = cell_size[1];
-    matrix[4] = cell_size[2];
-    matrix[5] = cell_size[4];
-    matrix[6] = cell_size[3];
-    matrix[7] = cell_size[4];
-    matrix[8] = cell_size[5];
-    // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
-    for (N[0]=-1;N[0]<=1;N[0]++)
-      for (N[1]=-1;N[1]<=1;N[1]++)
-        for (N[2]=-1;N[2]<=1;N[2]++) {
-          // create the translation vector
-          TranslationVector.Zero();
-          for (int i=NDIM;i--;)
-            TranslationVector[i] = (double)N[i];
-          TranslationVector.MatrixMultiplication(matrix);
-          // add onto the original vector to compare with
-          Shiftedy = y + TranslationVector;
-          // get distance and compare with minimum so far
-          tmp = DistanceSquared(Shiftedy);
-          if (tmp < res) res = tmp;
-        }
-    return (res);
-};
-
-/** Keeps the vector in a periodic cell, defined by the symmetric \a *matrix.
- * \param *out ofstream for debugging messages
- * Tries to translate a vector into each adjacent neighbouring cell.
- */
-void Vector::KeepPeriodic(const double * const matrix)
-{
-  //  int N[NDIM];
-  //  bool flag = false;
-    //vector Shifted, TranslationVector;
-  //  Log() << Verbose(1) << "Begin of KeepPeriodic." << endl;
-  //  Log() << Verbose(2) << "Vector is: ";
-  //  Output(out);
-  //  Log() << Verbose(0) << endl;
-    InverseMatrixMultiplication(matrix);
-    for(int i=NDIM;i--;) { // correct periodically
-      if (at(i) < 0) {  // get every coefficient into the interval [0,1)
-        at(i) += ceil(at(i));
-      } else {
-        at(i) -= floor(at(i));
-      }
-    }
-    MatrixMultiplication(matrix);
-  //  Log() << Verbose(2) << "New corrected vector is: ";
-  //  Output(out);
-  //  Log() << Verbose(0) << endl;
-  //  Log() << Verbose(1) << "End of KeepPeriodic." << endl;
-};
-
 /** Calculates scalar product between this and another vector.
  * \param *y array to second vector
@@ -199 +105 @@
 {
   double res = 0.;
-  for (int i=NDIM;i--;)
-    res += x[i]*y[i];
+  gsl_blas_ddot(content->content, y.content->content, &res);
   return (res);
 };
@@ -214 +119 @@
 {
   Vector tmp;
-  tmp[0] = x[1]* y[2] - x[2]* y[1];
-  tmp[1] = x[2]* y[0] - x[0]* y[2];
-  tmp[2] = x[0]* y[1] - x[1]* y[0];
+  for(int i=NDIM;i--;)
+    tmp[i] = at((i+1)%NDIM)*y[(i+2)%NDIM] - at((i+2)%NDIM)*y[(i+1)%NDIM];
   (*this) = tmp;
 };
@@ -309 +213 @@
 bool Vector::IsZero() const
 {
-  return (fabs(x[0])+fabs(x[1])+fabs(x[2]) < MYEPSILON);
+  return (fabs(at(0))+fabs(at(1))+fabs(at(2)) < MYEPSILON);
 };
 
@@ -338 +242 @@
   bool status = true;
   for (int i=0;i<NDIM;i++) {
-    if (fabs(x[i] - a[i]) > MYEPSILON)
+    if (fabs(at(i) - a[i]) > MYEPSILON)
       status = false;
   }
@@ -366 +270 @@
 double& Vector::operator[](size_t i){
   ASSERT(i<=NDIM && i>=0,"Vector Index out of Range");
-  return x[i];
+  return *gsl_vector_ptr (content->content, i);
 }
 
 const double& Vector::operator[](size_t i) const{
   ASSERT(i<=NDIM && i>=0,"Vector Index out of Range");
-  return x[i];
+  return *gsl_vector_ptr (content->content, i);
 }
 
@@ -382 +286 @@
 }
 
-double* Vector::get(){
-  return x;
+VectorContent* Vector::get(){
+  return content;
 }
 
@@ -498 +402 @@
 {
   for (int i=NDIM;i--;)
-    x[i] *= factor[i];
-};
-
+    at(i) *= factor[i];
+};
+
+void Vector::ScaleAll(const Vector &factor){
+  gsl_vector_mul(content->content, factor.content->content);
+}
 
 
 void Vector::Scale(const double factor)
 {
-  for (int i=NDIM;i--;)
-    x[i] *= factor;
-};
-
-/** Given a box by its matrix \a *M and its inverse *Minv the vector is made to point within that box.
- * \param *M matrix of box
- * \param *Minv inverse matrix
- */
-void Vector::WrapPeriodically(const double * const M, const double * const Minv)
-{
-  MatrixMultiplication(Minv);
-  // truncate to [0,1] for each axis
-  for (int i=0;i<NDIM;i++) {
-    //x[i] += 0.5;  // set to center of box
-    while (x[i] >= 1.)
-      x[i] -= 1.;
-    while (x[i] < 0.)
-      x[i] += 1.;
-  }
-  MatrixMultiplication(M);
+  gsl_vector_scale(content->content,factor);
 };
 
@@ -544 +432 @@
 }
 
-/** Do a matrix multiplication.
- * \param *matrix NDIM_NDIM array
- */
-void Vector::MatrixMultiplication(const double * const M)
-{
-  // do the matrix multiplication
-  at(0) = M[0]*x[0]+M[3]*x[1]+M[6]*x[2];
-  at(1) = M[1]*x[0]+M[4]*x[1]+M[7]*x[2];
-  at(2) = M[2]*x[0]+M[5]*x[1]+M[8]*x[2];
-};
-
-/** Do a matrix multiplication with the \a *A' inverse.
- * \param *matrix NDIM_NDIM array
- */
-bool Vector::InverseMatrixMultiplication(const double * const A)
-{
-  double B[NDIM*NDIM];
-  double detA = RDET3(A);
-  double detAReci;
-
-  // calculate the inverse B
-  if (fabs(detA) > MYEPSILON) {;  // RDET3(A) yields precisely zero if A irregular
-    detAReci = 1./detA;
-    B[0] =  detAReci*RDET2(A[4],A[5],A[7],A[8]);    // A_11
-    B[1] = -detAReci*RDET2(A[1],A[2],A[7],A[8]);    // A_12
-    B[2] =  detAReci*RDET2(A[1],A[2],A[4],A[5]);    // A_13
-    B[3] = -detAReci*RDET2(A[3],A[5],A[6],A[8]);    // A_21
-    B[4] =  detAReci*RDET2(A[0],A[2],A[6],A[8]);    // A_22
-    B[5] = -detAReci*RDET2(A[0],A[2],A[3],A[5]);    // A_23
-    B[6] =  detAReci*RDET2(A[3],A[4],A[6],A[7]);    // A_31
-    B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]);    // A_32
-    B[8] =  detAReci*RDET2(A[0],A[1],A[3],A[4]);    // A_33
-
-    // do the matrix multiplication
-    at(0) = B[0]*x[0]+B[3]*x[1]+B[6]*x[2];
-    at(1) = B[1]*x[0]+B[4]*x[1]+B[7]*x[2];
-    at(2) = B[2]*x[0]+B[5]*x[1]+B[8]*x[2];
-
-    return true;
-  } else {
-    return false;
-  }
-};
-
-
 /** Creates this vector as the b y *factors' components scaled linear combination of the given three.
  * this vector = x1*factors[0] + x2* factors[1] + x3*factors[2]
@@ -616 +459 @@
   SubtractVector(x1);
   for (int i=NDIM;i--;)
-    result = result || (fabs(x[i]) > MYEPSILON);
+    result = result || (fabs(at(i)) > MYEPSILON);
 
   return result;
@@ -677 +520 @@
 void Vector::AddVector(const Vector &y)
 {
-  for(int i=NDIM;i--;)
-    x[i] += y[i];
+  gsl_vector_add(content->content, y.content->content);
 }
 
@@ -686 +528 @@
 void Vector::SubtractVector(const Vector &y)
 {
-  for(int i=NDIM;i--;)
-    x[i] -= y[i];
-}
-
-/**
- * Checks whether this vector is within the parallelepiped defined by the given three vectors and
- * their offset.
- *
- * @param offest for the origin of the parallelepiped
- * @param three vectors forming the matrix that defines the shape of the parallelpiped
- */
-bool Vector::IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const
-{
-  Vector a = (*this)-offset;
-  a.InverseMatrixMultiplication(parallelepiped);
-  bool isInside = true;
-
-  for (int i=NDIM;i--;)
-    isInside = isInside && ((a[i] <= 1) && (a[i] >= 0));
-
-  return isInside;
+  gsl_vector_sub(content->content, y.content->content);
 }