Changeset 8f822c for src/vector.cpp

Timestamp:

Jul 1, 2010, 1:55:08 PM (15 years ago)

Author:

Tillmann Crueger <crueger@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, Candidate_v1.7.0, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable

Children:

b5bf84

Parents:

0d5dce (diff), 014475 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

Message:

Merge branch 'VectorRefactoring' into Shapes

File:

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

src/vector.cpp

@@ -8 +8 @@
 
 #include "vector.hpp"
+#include "Matrix.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;
@@ -34 +38 @@
 {
   content = gsl_vector_alloc(NDIM);
-  gsl_vector_set(content,0,src[0]);
-  gsl_vector_set(content,1,src[1]);
-  gsl_vector_set(content,2,src[2]);
+  gsl_vector_memcpy(content, src.content);
 }
 
@@ -49 +51 @@
 };
 
+Vector::Vector(gsl_vector *_content) :
+  content(_content)
+{}
+
 /**
  * Assignment operator
@@ -55 +61 @@
   // check for self assignment
   if(&src!=this){
-    gsl_vector_set(content,0,src[0]);
-    gsl_vector_set(content,1,src[1]);
-    gsl_vector_set(content,2,src[2]);
+    gsl_vector_memcpy(content, src.content);
   }
   return *this;
@@ -94 +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
@@ -203 +105 @@
 {
   double res = 0.;
-  for (int i=NDIM;i--;)
-    res += at(i)*y[i];
+  gsl_blas_ddot(content, y.content, &res);
   return (res);
 };
@@ -461 +362 @@
 };
 
+Vector &Vector::operator*=(const Matrix &mat){
+  (*this) = mat*(*this);
+  return *this;
+}
+
+Vector operator*(const Matrix &mat,const Vector &vec){
+  gsl_vector *res = gsl_vector_calloc(NDIM);
+  gsl_blas_dgemv( CblasNoTrans, 1.0, mat.content, vec.content, 0.0, res);
+  return Vector(res);
+}
+
+
 /** Factors given vector \a a times \a m.
  * \param a vector
@@ -508 +421 @@
 void Vector::Scale(const double factor)
 {
-  for (int i=NDIM;i--;)
-    at(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++) {
-    //at(i) += 0.5;  // set to center of box
-    while (at(i) >= 1.)
-      at(i) -= 1.;
-    while (at(i) < 0.)
-      at(i) += 1.;
-  }
-  MatrixMultiplication(M);
+  gsl_vector_scale(content,factor);
 };
 
@@ -546 +440 @@
   return make_pair(res,helper);
 }
-
-/** Do a matrix multiplication.
- * \param *matrix NDIM_NDIM array
- */
-void Vector::MatrixMultiplication(const double * const M)
-{
-  Vector tmp;
-  // do the matrix multiplication
-  for(int i=NDIM;i--;)
-    tmp[i] = M[i]*at(0)+M[i+3]*at(1)+M[i+6]*at(2);
-
-  (*this) = tmp;
-};
-
-/** 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
-
-    MatrixMultiplication(B);
-
-    return true;
-  } else {
-    return false;
-  }
-};
-
 
 /** Creates this vector as the b y *factors' components scaled linear combination of the given three.
@@ -679 +529 @@
 void Vector::AddVector(const Vector &y)
 {
-  for(int i=NDIM;i--;)
-    at(i) += y[i];
+  gsl_vector_add(content, y.content);
 }
 
@@ -688 +537 @@
 void Vector::SubtractVector(const Vector &y)
 {
-  for(int i=NDIM;i--;)
-    at(i) -= y[i];
+  gsl_vector_sub(content, y.content);
 }
 
@@ -699 +547 @@
  * @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 Vector::IsInParallelepiped(const Vector &offset, const Matrix& _parallelepiped) const
+{
+  Matrix parallelepiped = _parallelepiped.invert();
+  Vector a = parallelepiped * ((*this)-offset);
   bool isInside = true;