Changes in / [c7b39a:f0f1cc]


Ignore:
Location:
src
Files:
4 deleted
10 edited

Legend:

Unmodified
Added
Removed

  • src/Actions/WorldAction/RepeatBoxAction.cpp

    rc7b39a rf0f1cc  
    1313#include "molecule.hpp"
    1414#include "vector.hpp"
    15 #include "Matrix.hpp"
    1615#include "verbose.hpp"
    1716#include "World.hpp"
    … …  
    6160    (cout << "Repeating box " << Repeater << " times for (x,y,z) axis." << endl);
    6261    double * const cell_size = World::getInstance().getDomain();
    63     double *M_double = ReturnFullMatrixforSymmetric(cell_size);
    64     Matrix M = Matrix(M_double);
    65     delete[] M_double;
     62    double *M = ReturnFullMatrixforSymmetric(cell_size);
    6663    Vector x,y;
    6764    int n[NDIM];
    … …  
    121118      }
    122119    }
     120    delete(M);
    123121    delete dialog;
    124122    return Action::success;

  • src/Makefile.am

    rc7b39a rf0f1cc  
    3030  atom_trajectoryparticleinfo.hpp
    3131
    32 EXCEPTIONSOURCE = Exceptions/CustomException.cpp \
    33                                   Exceptions/LinearDependenceException.cpp \
    34                                   Exceptions/MathException.cpp \
    35                                   Exceptions/NotInvertibleException.cpp \
    36                                   Exceptions/SkewException.cpp \
    37                                   Exceptions/ZeroVectorException.cpp
    38                                  
    39 EXCEPTIONHEADER = Exceptions/CustomException.hpp \
    40                                   Exceptions/LinearDependenceException.hpp \
    41                                   Exceptions/MathException.hpp \
    42                                   Exceptions/NotInvertibleException.hpp \
    43                                   Exceptions/SkewException.hpp \
    44                                   Exceptions/ZeroVectorException.hpp
    45 
    46 # TODO: Move Exceptionsource back down, when transition is done
    4732LINALGSOURCE = \
    48   ${EXCEPTIONSOURCE} \
    4933  ${HELPERSOURCE} \
    5034  gslmatrix.cpp \
    5135  gslvector.cpp \
    5236  linearsystemofequations.cpp \
    53   Matrix.cpp \
    5437  Space.cpp \
    5538  vector.cpp
    56  
    5739                           
    58 LINALGHEADER = \
    59   ${EXCEPTIONHEADER} \
    60   gslmatrix.hpp \
     40LINALGHEADER = gslmatrix.hpp \
    6141  gslvector.hpp \
    6242  linearsystemofequations.hpp \
    63   Matrix.hpp \
    6443  Space.hpp \
    6544  vector.hpp
    … …  
    139118  Descriptors/MoleculeNameDescriptor.hpp \
    140119  Descriptors/MoleculePtrDescriptor.hpp
    141 
    142 #TODO: Exceptionsource should go here, when the transisition makes this possible                                   
     120                                   
     121EXCEPTIONSOURCE = Exceptions/CustomException.cpp \
     122                                  Exceptions/LinearDependenceException.cpp \
     123                                  Exceptions/MathException.cpp \
     124                                  Exceptions/SkewException.cpp \
     125                                  Exceptions/ZeroVectorException.cpp
     126                                 
     127EXCEPTIONHEADER = Exceptions/CustomException.hpp \
     128                                  Exceptions/LinearDependenceException.hpp \
     129                                  Exceptions/MathException.hpp \
     130                                  Exceptions/SkewException.hpp \
     131                                  Exceptions/ZeroVectorException.hpp
     132
    143133SOURCE = \
    144134  ${ANALYSISSOURCE} \
    … …  
    149139  ${DESCRIPTORSOURCE} \
    150140  ${HELPERSOURCE} \
     141  ${EXCEPTIONSOURCE} \
    151142  bond.cpp \
    152143  bondgraph.cpp \
    … …  
    194185  ${PATTERNHEADER} \
    195186  ${DESCRIPTORHEADER} \
     187  ${EXCEPTIONHEADER} \
    196188  bond.hpp \
    197189  bondgraph.hpp \

  • src/analysis_correlation.cpp

    rc7b39a rf0f1cc  
    1919#include "triangleintersectionlist.hpp"
    2020#include "vector.hpp"
    21 #include "Matrix.hpp"
    2221#include "verbose.hpp"
    2322#include "World.hpp"
    … …  
    136135  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++){
    137136    if ((*MolWalker)->ActiveFlag) {
    138       double * FullMatrix_double = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
    139       Matrix FullMatrix = Matrix(FullMatrix_double);
    140       Matrix FullInverseMatrix = FullMatrix.invert();
    141       delete[](FullMatrix_double);
     137      double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
     138      double * FullInverseMatrix = InverseMatrix(FullMatrix);
    142139      DoeLog(2) && (eLog()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    143140      eLog() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
    … …  
    179176            }
    180177      }
     178      delete[](FullMatrix);
     179      delete[](FullInverseMatrix);
    181180    }
    182181  }
    … …  
    247246  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    248247    if ((*MolWalker)->ActiveFlag) {
    249       double * FullMatrix_double = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
    250       Matrix FullMatrix = Matrix(FullMatrix_double);
    251       Matrix FullInverseMatrix = FullMatrix.invert();
    252       delete[] FullMatrix_double;
     248      double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
     249      double * FullInverseMatrix = InverseMatrix(FullMatrix);
    253250      DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    254251      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    … …  
    270267          }
    271268      }
     269      delete[](FullMatrix);
     270      delete[](FullInverseMatrix);
    272271    }
    273272
    … …  
    353352  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    354353    if ((*MolWalker)->ActiveFlag) {
    355       double * FullMatrix_double = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
    356       Matrix FullMatrix = Matrix(FullMatrix_double);
    357       Matrix FullInverseMatrix = FullMatrix.invert();
    358       delete[](FullMatrix_double);
     354      double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
     355      double * FullInverseMatrix = InverseMatrix(FullMatrix);
    359356      DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    360357      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    … …  
    384381          }
    385382      }
     383      delete[](FullMatrix);
     384      delete[](FullInverseMatrix);
    386385    }
    387386

  • src/boundary.cpp

    rc7b39a rf0f1cc  
    2222#include "World.hpp"
    2323#include "Plane.hpp"
    24 #include "Matrix.hpp"
    2524
    2625#include<gsl/gsl_poly.h>
    … …  
    765764  int N[NDIM];
    766765  int n[NDIM];
    767   double *M_double =  ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
    768   Matrix M = Matrix(M_double);
    769   delete[](M_double);
    770   Matrix Rotations;
    771   Matrix MInverse = M.invert();
     766  double *M =  ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
     767  double Rotations[NDIM*NDIM];
     768  double *MInverse = InverseMatrix(M);
    772769  Vector AtomTranslations;
    773770  Vector FillerTranslations;
    … …  
    803800  // calculate filler grid in [0,1]^3
    804801  FillerDistance = Vector(distance[0], distance[1], distance[2]);
    805   FillerDistance.MatrixMultiplication(MInverse);
     802  FillerDistance.InverseMatrixMultiplication(M);
    806803  for(int i=0;i<NDIM;i++)
    807804    N[i] = (int) ceil(1./FillerDistance[i]);
    … …  
    838835            }
    839836
    840             Rotations.at(0,0) =   cos(phi[0])            *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2]));
    841             Rotations.at(0,1) =   sin(phi[0])            *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2]));
    842             Rotations.at(0,2) =               cos(phi[1])*sin(phi[2])                                     ;
    843             Rotations.at(1,0) = - sin(phi[0])*cos(phi[1])                                                ;
    844             Rotations.at(1,1) =   cos(phi[0])*cos(phi[1])                                                ;
    845             Rotations.at(1,2) =               sin(phi[1])                                                ;
    846             Rotations.at(2,0) = - cos(phi[0])            *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2]));
    847             Rotations.at(2,1) = - sin(phi[0])            *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2]));
    848             Rotations.at(2,2) =               cos(phi[1])*cos(phi[2])                                     ;
     837            Rotations[0] =   cos(phi[0])            *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2]));
     838            Rotations[3] =   sin(phi[0])            *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2]));
     839            Rotations[6] =               cos(phi[1])*sin(phi[2])                                     ;
     840            Rotations[1] = - sin(phi[0])*cos(phi[1])                                                ;
     841            Rotations[4] =   cos(phi[0])*cos(phi[1])                                                ;
     842            Rotations[7] =               sin(phi[1])                                                ;
     843            Rotations[3] = - cos(phi[0])            *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2]));
     844            Rotations[5] = - sin(phi[0])            *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2]));
     845            Rotations[8] =               cos(phi[1])*cos(phi[2])                                     ;
    849846          }
    850847
    … …  
    895892            }
    896893      }
     894  delete[](M);
     895  delete[](MInverse);
    897896
    898897  return Filling;

  • src/ellipsoid.cpp

    rc7b39a rf0f1cc  
    2121#include "tesselation.hpp"
    2222#include "vector.hpp"
    23 #include "Matrix.hpp"
    2423#include "verbose.hpp"
    2524
    … …  
    3534  Vector helper, RefPoint;
    3635  double distance = -1.;
    37   Matrix Matrix;
     36  double Matrix[NDIM*NDIM];
    3837  double InverseLength[3];
    3938  double psi,theta,phi; // euler angles in ZX'Z'' convention
    … …  
    5251  theta = EllipsoidAngle[1];
    5352  phi = EllipsoidAngle[2];
    54   Matrix.at(0,0) = cos(psi)*cos(phi) - sin(psi)*cos(theta)*sin(phi);
    55   Matrix.at(1,0) = -cos(psi)*sin(phi) - sin(psi)*cos(theta)*cos(phi);
    56   Matrix.at(2,0) = sin(psi)*sin(theta);
    57   Matrix.at(0,1) = sin(psi)*cos(phi) + cos(psi)*cos(theta)*sin(phi);
    58   Matrix.at(1,1) = cos(psi)*cos(theta)*cos(phi) - sin(psi)*sin(phi);
    59   Matrix.at(2,1) = -cos(psi)*sin(theta);
    60   Matrix.at(0,2) = sin(theta)*sin(phi);
    61   Matrix.at(1,2) = sin(theta)*cos(phi);
    62   Matrix.at(2,2) = cos(theta);
     53  Matrix[0] = cos(psi)*cos(phi) - sin(psi)*cos(theta)*sin(phi);
     54  Matrix[1] = -cos(psi)*sin(phi) - sin(psi)*cos(theta)*cos(phi);
     55  Matrix[2] = sin(psi)*sin(theta);
     56  Matrix[3] = sin(psi)*cos(phi) + cos(psi)*cos(theta)*sin(phi);
     57  Matrix[4] = cos(psi)*cos(theta)*cos(phi) - sin(psi)*sin(phi);
     58  Matrix[5] = -cos(psi)*sin(theta);
     59  Matrix[6] = sin(theta)*sin(phi);
     60  Matrix[7] = sin(theta)*cos(phi);
     61  Matrix[8] = cos(theta);
    6362  helper.MatrixMultiplication(Matrix);
    6463  helper.ScaleAll(InverseLength);
    … …  
    7473  phi = -EllipsoidAngle[2];
    7574  helper.ScaleAll(EllipsoidLength);
    76   Matrix.at(0,0) = cos(psi)*cos(phi) - sin(psi)*cos(theta)*sin(phi);
    77   Matrix.at(1,0) = -cos(psi)*sin(phi) - sin(psi)*cos(theta)*cos(phi);
    78   Matrix.at(2,0) = sin(psi)*sin(theta);
    79   Matrix.at(0,1) = sin(psi)*cos(phi) + cos(psi)*cos(theta)*sin(phi);
    80   Matrix.at(1,1) = cos(psi)*cos(theta)*cos(phi) - sin(psi)*sin(phi);
    81   Matrix.at(2,1) = -cos(psi)*sin(theta);
    82   Matrix.at(0,2) = sin(theta)*sin(phi);
    83   Matrix.at(1,2) = sin(theta)*cos(phi);
    84   Matrix.at(2,2) = cos(theta);
     75  Matrix[0] = cos(psi)*cos(phi) - sin(psi)*cos(theta)*sin(phi);
     76  Matrix[1] = -cos(psi)*sin(phi) - sin(psi)*cos(theta)*cos(phi);
     77  Matrix[2] = sin(psi)*sin(theta);
     78  Matrix[3] = sin(psi)*cos(phi) + cos(psi)*cos(theta)*sin(phi);
     79  Matrix[4] = cos(psi)*cos(theta)*cos(phi) - sin(psi)*sin(phi);
     80  Matrix[5] = -cos(psi)*sin(theta);
     81  Matrix[6] = sin(theta)*sin(phi);
     82  Matrix[7] = sin(theta)*cos(phi);
     83  Matrix[8] = cos(theta);
    8584  helper.MatrixMultiplication(Matrix);
    8685  //Log() << Verbose(4) << "Intersection is at " << helper << "." << endl;

  • src/molecule.cpp

    rc7b39a rf0f1cc  
    2727#include "tesselation.hpp"
    2828#include "vector.hpp"
    29 #include "Matrix.hpp"
    3029#include "World.hpp"
    3130#include "Plane.hpp"
    … …  
    285284  Vector Orthovector1, Orthovector2;  // temporary vectors in coordination construction
    286285  Vector InBondvector;    // vector in direction of *Bond
    287   Matrix matrix;
     286  double *matrix = NULL;
    288287  bond *Binder = NULL;
    289288  double * const cell_size = World::getInstance().getDomain();
    … …  
    308307      } // (signs are correct, was tested!)
    309308    }
    310     double *matrix_double = ReturnFullMatrixforSymmetric(cell_size);
    311     matrix = Matrix(matrix_double);
    312     delete[](matrix_double);
     309    matrix = ReturnFullMatrixforSymmetric(cell_size);
    313310    Orthovector1.MatrixMultiplication(matrix);
    314311    InBondvector -= Orthovector1; // subtract just the additional translation
     312    delete[](matrix);
    315313    bondlength = InBondvector.Norm();
    316314//    Log() << Verbose(4) << "Corrected InBondvector is now: ";
    … …  
    543541      break;
    544542  }
     543  delete[](matrix);
    545544
    546545//  Log() << Verbose(3) << "End of AddHydrogenReplacementAtom." << endl;

  • src/molecule_fragmentation.cpp

    rc7b39a rf0f1cc  
    2222#include "periodentafel.hpp"
    2323#include "World.hpp"
    24 #include "Matrix.hpp"
    2524
    2625/************************************* Functions for class molecule *********************************/
    … …  
    17101709  atom *OtherWalker = NULL;
    17111710  double * const cell_size = World::getInstance().getDomain();
    1712   double *matrix_double = ReturnFullMatrixforSymmetric(cell_size);
    1713   Matrix matrix = Matrix(matrix_double);
    1714   delete[](matrix_double);
     1711  double *matrix = ReturnFullMatrixforSymmetric(cell_size);
    17151712  enum Shading *ColorList = NULL;
    17161713  double tmp;
    … …  
    17831780  delete(AtomStack);
    17841781  delete[](ColorList);
     1782  delete[](matrix);
    17851783  DoLog(2) && (Log() << Verbose(2) << "End of ScanForPeriodicCorrection." << endl);
    17861784};

  • src/molecule_geometry.cpp

    rc7b39a rf0f1cc  
    1919#include "World.hpp"
    2020#include "Plane.hpp"
    21 #include "Matrix.hpp"
    2221#include <boost/foreach.hpp>
    2322
    … …  
    155154
    156155  const double *cell_size = World::getInstance().getDomain();
    157   double *M_double = ReturnFullMatrixforSymmetric(cell_size);
    158   Matrix M = Matrix(M_double);
    159   delete[](M_double);
     156  double *M = ReturnFullMatrixforSymmetric(cell_size);
    160157  a->MatrixMultiplication(M);
     158  delete[](M);
    161159
    162160  return a;
    … …  
    276274{
    277275  double * const cell_size = World::getInstance().getDomain();
    278   double *matrix_double = ReturnFullMatrixforSymmetric(cell_size);
    279   Matrix matrix = Matrix(matrix_double);
    280   delete[](matrix_double);
    281   Matrix inversematrix = matrix.invert();
     276  double *matrix = ReturnFullMatrixforSymmetric(cell_size);
     277  double *inversematrix = InverseMatrix(matrix);
    282278  double tmp;
    283279  bool flag;
    … …  
    328324    }
    329325  } while (!flag);
     326  delete[](matrix);
     327  delete[](inversematrix);
    330328
    331329  Center.Scale(1./static_cast<double>(getAtomCount()));
    … …  
    389387    DoLog(1) && (Log() << Verbose(1) << "Transforming molecule into PAS ... ");
    390388    // the eigenvectors specify the transformation matrix
    391     Matrix M = Matrix(evec->data);
    392     ActOnAllVectors( &Vector::MatrixMultiplication, static_cast<const Matrix>(M));
     389    ActOnAllVectors( &Vector::MatrixMultiplication, (const double *) evec->data );
    393390    DoLog(0) && (Log() << Verbose(0) << "done." << endl);
    394391

  • src/vector.cpp

    rc7b39a rf0f1cc  
    88
    99#include "vector.hpp"
    10 #include "Matrix.hpp"
    1110#include "verbose.hpp"
    1211#include "World.hpp"
    1312#include "Helpers/Assert.hpp"
    1413#include "Helpers/fast_functions.hpp"
    15 #include "Exceptions/MathException.hpp"
    1614
    1715#include <iostream>
    18 #include <gsl/gsl_blas.h>
    19 
    2016
    2117using namespace std;
    … …  
    3834{
    3935  content = gsl_vector_alloc(NDIM);
    40   gsl_vector_memcpy(content, src.content);
     36  gsl_vector_set(content,0,src[0]);
     37  gsl_vector_set(content,1,src[1]);
     38  gsl_vector_set(content,2,src[2]);
    4139}
    4240
    … …  
    5149};
    5250
    53 Vector::Vector(gsl_vector *_content) :
    54   content(_content)
    55 {}
    56 
    5751/**
    5852 * Assignment operator
    … …  
    6155  // check for self assignment
    6256  if(&src!=this){
    63     gsl_vector_memcpy(content, src.content);
     57    gsl_vector_set(content,0,src[0]);
     58    gsl_vector_set(content,1,src[1]);
     59    gsl_vector_set(content,2,src[2]);
    6460  }
    6561  return *this;
    … …  
    105101double Vector::PeriodicDistance(const Vector &y, const double * const cell_size) const
    106102{
    107   double res = distance(y), tmp;
    108   Matrix matrix;
     103  double res = distance(y), tmp, matrix[NDIM*NDIM];
    109104    Vector Shiftedy, TranslationVector;
    110105    int N[NDIM];
    111     matrix.at(0,0) = cell_size[0];
    112     matrix.at(1,0) = cell_size[1];
    113     matrix.at(2,0) = cell_size[3];
    114     matrix.at(0,1) = cell_size[1];
    115     matrix.at(1,1) = cell_size[2];
    116     matrix.at(2,1) = cell_size[4];
    117     matrix.at(0,2) = cell_size[3];
    118     matrix.at(1,2) = cell_size[4];
    119     matrix.at(2,2) = cell_size[5];
     106    matrix[0] = cell_size[0];
     107    matrix[1] = cell_size[1];
     108    matrix[2] = cell_size[3];
     109    matrix[3] = cell_size[1];
     110    matrix[4] = cell_size[2];
     111    matrix[5] = cell_size[4];
     112    matrix[6] = cell_size[3];
     113    matrix[7] = cell_size[4];
     114    matrix[8] = cell_size[5];
    120115    // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
    121116    for (N[0]=-1;N[0]<=1;N[0]++)
    … …  
    143138double Vector::PeriodicDistanceSquared(const Vector &y, const double * const cell_size) const
    144139{
    145   double res = DistanceSquared(y), tmp;
    146   Matrix matrix;
     140  double res = DistanceSquared(y), tmp, matrix[NDIM*NDIM];
    147141    Vector Shiftedy, TranslationVector;
    148142    int N[NDIM];
    149     matrix.at(0,0) = cell_size[0];
    150     matrix.at(1,0) = cell_size[1];
    151     matrix.at(2,0) = cell_size[3];
    152     matrix.at(0,1) = cell_size[1];
    153     matrix.at(1,1) = cell_size[2];
    154     matrix.at(2,1) = cell_size[4];
    155     matrix.at(0,2) = cell_size[3];
    156     matrix.at(1,2) = cell_size[4];
    157     matrix.at(2,2) = cell_size[5];
     143    matrix[0] = cell_size[0];
     144    matrix[1] = cell_size[1];
     145    matrix[2] = cell_size[3];
     146    matrix[3] = cell_size[1];
     147    matrix[4] = cell_size[2];
     148    matrix[5] = cell_size[4];
     149    matrix[6] = cell_size[3];
     150    matrix[7] = cell_size[4];
     151    matrix[8] = cell_size[5];
    158152    // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
    159153    for (N[0]=-1;N[0]<=1;N[0]++)
    … …  
    178172 * Tries to translate a vector into each adjacent neighbouring cell.
    179173 */
    180 void Vector::KeepPeriodic(const double * const _matrix)
    181 {
    182   Matrix matrix = Matrix(_matrix);
     174void Vector::KeepPeriodic(const double * const matrix)
     175{
    183176  //  int N[NDIM];
    184177  //  bool flag = false;
    … …  
    188181  //  Output(out);
    189182  //  Log() << Verbose(0) << endl;
    190     MatrixMultiplication(matrix.invert());
     183    InverseMatrixMultiplication(matrix);
    191184    for(int i=NDIM;i--;) { // correct periodically
    192185      if (at(i) < 0) {  // get every coefficient into the interval [0,1)
    … …  
    210203{
    211204  double res = 0.;
    212   gsl_blas_ddot(content, y.content, &res);
     205  for (int i=NDIM;i--;)
     206    res += at(i)*y[i];
    213207  return (res);
    214208};
    … …  
    467461};
    468462
    469 Vector &Vector::operator*=(const Matrix &mat){
    470   (*this) = mat*(*this);
    471   return *this;
    472 }
    473 
    474 Vector operator*(const Matrix &mat,const Vector &vec){
    475   gsl_vector *res = gsl_vector_calloc(NDIM);
    476   gsl_blas_dgemv( CblasNoTrans, 1.0, mat.content, vec.content, 0.0, res);
    477   return Vector(res);
    478 }
    479 
    480 
    481463/** Factors given vector \a a times \a m.
    482464 * \param a vector
    … …  
    526508void Vector::Scale(const double factor)
    527509{
    528   gsl_vector_scale(content,factor);
     510  for (int i=NDIM;i--;)
     511    at(i) *= factor;
    529512};
    530513
    … …  
    533516 * \param *Minv inverse matrix
    534517 */
    535 void Vector::WrapPeriodically(const double * const _M, const double * const _Minv)
    536 {
    537   Matrix M = Matrix(_M);
    538   Matrix Minv = Matrix(_Minv);
     518void Vector::WrapPeriodically(const double * const M, const double * const Minv)
     519{
    539520  MatrixMultiplication(Minv);
    540521  // truncate to [0,1] for each axis
    … …  
    569550 * \param *matrix NDIM_NDIM array
    570551 */
    571 void Vector::MatrixMultiplication(const Matrix &M)
    572 {
    573   (*this) *= M;
     552void Vector::MatrixMultiplication(const double * const M)
     553{
     554  Vector tmp;
     555  // do the matrix multiplication
     556  for(int i=NDIM;i--;)
     557    tmp[i] = M[i]*at(0)+M[i+3]*at(1)+M[i+6]*at(2);
     558
     559  (*this) = tmp;
    574560};
    575561
    … …  
    579565bool Vector::InverseMatrixMultiplication(const double * const A)
    580566{
    581   /*
    582567  double B[NDIM*NDIM];
    583568  double detA = RDET3(A);
    … …  
    601586    return true;
    602587  } else {
    603     return false;
    604   }
    605   */
    606   Matrix mat = Matrix(A);
    607   try{
    608     (*this) *= mat.invert();
    609     return true;
    610   }
    611   catch(MathException &excpt){
    612588    return false;
    613589  }
    … …  
    703679void Vector::AddVector(const Vector &y)
    704680{
    705   gsl_vector_add(content, y.content);
     681  for(int i=NDIM;i--;)
     682    at(i) += y[i];
    706683}
    707684
    … …  
    711688void Vector::SubtractVector(const Vector &y)
    712689{
    713   gsl_vector_sub(content, y.content);
     690  for(int i=NDIM;i--;)
     691    at(i) -= y[i];
    714692}
    715693

  • src/vector.hpp

    rc7b39a rf0f1cc  
    2424
    2525class Vector;
    26 class Matrix;
    2726
    2827typedef std::vector<Vector> pointset;
    … …  
    3231 */
    3332class Vector : public Space{
    34   friend Vector operator*(const Matrix&,const Vector&);
    3533public:
     34
    3635  Vector();
    3736  Vector(const double x1, const double x2, const double x3);
    … …  
    6059  void ScaleAll(const double *factor);
    6160  void Scale(const double factor);
    62   void MatrixMultiplication(const Matrix &M);
     61  void MatrixMultiplication(const double * const M);
    6362  bool InverseMatrixMultiplication(const double * const M);
    6463  void KeepPeriodic(const double * const matrix);
    … …  
    9796  Vector const operator+(const Vector& b) const;
    9897  Vector const operator-(const Vector& b) const;
    99   Vector& operator*=(const Matrix&);
    10098
    10199  // Methods inherited from Space
    … …  
    106104
    107105private:
    108   Vector(gsl_vector *);
    109106  gsl_vector *content;
    110107
    … …  
    121118Vector const operator*(const Vector& a, const double m);
    122119Vector const operator*(const double m, const Vector& a);
    123 Vector operator*(const Matrix&,const Vector&);
    124 
    125120
    126121#endif /*VECTOR_HPP_*/
Note: See TracChangeset for help on using the changeset viewer.