Changes in / [0c5eeb:33d774]


Ignore:
Location:
src
Files:
4 added
10 edited

Legend:

Unmodified
Added
Removed

  • src/Actions/WorldAction/RepeatBoxAction.cpp

    r0c5eeb r33d774  
    1313#include "molecule.hpp"
    1414#include "vector.hpp"
     15#include "Matrix.hpp"
    1516#include "verbose.hpp"
    1617#include "World.hpp"
    … …  
    6061    (cout << "Repeating box " << Repeater << " times for (x,y,z) axis." << endl);
    6162    double * const cell_size = World::getInstance().getDomain();
    62     double *M = ReturnFullMatrixforSymmetric(cell_size);
     63    double *M_double = ReturnFullMatrixforSymmetric(cell_size);
     64    Matrix M = Matrix(M_double);
     65    delete[] M_double;
    6366    Vector x,y;
    6467    int n[NDIM];
    … …  
    118121      }
    119122    }
    120     delete(M);
    121123    delete dialog;
    122124    return Action::success;

  • src/Makefile.am

    r0c5eeb r33d774  
    3030  atom_trajectoryparticleinfo.hpp
    3131
     32EXCEPTIONSOURCE = Exceptions/CustomException.cpp \
     33                                  Exceptions/LinearDependenceException.cpp \
     34                                  Exceptions/MathException.cpp \
     35                                  Exceptions/NotInvertibleException.cpp \
     36                                  Exceptions/SkewException.cpp \
     37                                  Exceptions/ZeroVectorException.cpp
     38                                 
     39EXCEPTIONHEADER = 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
    3247LINALGSOURCE = \
     48  ${EXCEPTIONSOURCE} \
    3349  ${HELPERSOURCE} \
    3450  gslmatrix.cpp \
    3551  gslvector.cpp \
    3652  linearsystemofequations.cpp \
     53  Matrix.cpp \
    3754  Space.cpp \
    3855  vector.cpp
     56 
    3957                           
    40 LINALGHEADER = gslmatrix.hpp \
     58LINALGHEADER = \
     59  ${EXCEPTIONHEADER} \
     60  gslmatrix.hpp \
    4161  gslvector.hpp \
    4262  linearsystemofequations.hpp \
     63  Matrix.hpp \
    4364  Space.hpp \
    4465  vector.hpp
    … …  
    118139  Descriptors/MoleculeNameDescriptor.hpp \
    119140  Descriptors/MoleculePtrDescriptor.hpp
    120                                    
    121 EXCEPTIONSOURCE = Exceptions/CustomException.cpp \
    122                                   Exceptions/LinearDependenceException.cpp \
    123                                   Exceptions/MathException.cpp \
    124                                   Exceptions/SkewException.cpp \
    125                                   Exceptions/ZeroVectorException.cpp
    126                                  
    127 EXCEPTIONHEADER = Exceptions/CustomException.hpp \
    128                                   Exceptions/LinearDependenceException.hpp \
    129                                   Exceptions/MathException.hpp \
    130                                   Exceptions/SkewException.hpp \
    131                                   Exceptions/ZeroVectorException.hpp
    132 
     141
     142#TODO: Exceptionsource should go here, when the transisition makes this possible                                   
    133143SOURCE = \
    134144  ${ANALYSISSOURCE} \
    … …  
    139149  ${DESCRIPTORSOURCE} \
    140150  ${HELPERSOURCE} \
    141   ${EXCEPTIONSOURCE} \
    142151  bond.cpp \
    143152  bondgraph.cpp \
    … …  
    185194  ${PATTERNHEADER} \
    186195  ${DESCRIPTORHEADER} \
    187   ${EXCEPTIONHEADER} \
    188196  bond.hpp \
    189197  bondgraph.hpp \

  • src/analysis_correlation.cpp

    r0c5eeb r33d774  
    1919#include "triangleintersectionlist.hpp"
    2020#include "vector.hpp"
     21#include "Matrix.hpp"
    2122#include "verbose.hpp"
    2223#include "World.hpp"
    … …  
    135136  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++){
    136137    if ((*MolWalker)->ActiveFlag) {
    137       double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
    138       double * FullInverseMatrix = InverseMatrix(FullMatrix);
     138      double * FullMatrix_double = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
     139      Matrix FullMatrix = Matrix(FullMatrix_double);
     140      Matrix FullInverseMatrix = FullMatrix.invert();
     141      delete[](FullMatrix_double);
    139142      DoeLog(2) && (eLog()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    140143      eLog() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
    … …  
    176179            }
    177180      }
    178       delete[](FullMatrix);
    179       delete[](FullInverseMatrix);
    180181    }
    181182  }
    … …  
    246247  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    247248    if ((*MolWalker)->ActiveFlag) {
    248       double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
    249       double * FullInverseMatrix = InverseMatrix(FullMatrix);
     249      double * FullMatrix_double = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
     250      Matrix FullMatrix = Matrix(FullMatrix_double);
     251      Matrix FullInverseMatrix = FullMatrix.invert();
     252      delete[] FullMatrix_double;
    250253      DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    251254      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    … …  
    267270          }
    268271      }
    269       delete[](FullMatrix);
    270       delete[](FullInverseMatrix);
    271272    }
    272273
    … …  
    352353  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    353354    if ((*MolWalker)->ActiveFlag) {
    354       double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
    355       double * FullInverseMatrix = InverseMatrix(FullMatrix);
     355      double * FullMatrix_double = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
     356      Matrix FullMatrix = Matrix(FullMatrix_double);
     357      Matrix FullInverseMatrix = FullMatrix.invert();
     358      delete[](FullMatrix_double);
    356359      DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    357360      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    … …  
    381384          }
    382385      }
    383       delete[](FullMatrix);
    384       delete[](FullInverseMatrix);
    385386    }
    386387

  • src/boundary.cpp

    r0c5eeb r33d774  
    2222#include "World.hpp"
    2323#include "Plane.hpp"
     24#include "Matrix.hpp"
    2425
    2526#include<gsl/gsl_poly.h>
    … …  
    764765  int N[NDIM];
    765766  int n[NDIM];
    766   double *M =  ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
    767   double Rotations[NDIM*NDIM];
    768   double *MInverse = InverseMatrix(M);
     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();
    769772  Vector AtomTranslations;
    770773  Vector FillerTranslations;
    … …  
    800803  // calculate filler grid in [0,1]^3
    801804  FillerDistance = Vector(distance[0], distance[1], distance[2]);
    802   FillerDistance.InverseMatrixMultiplication(M);
     805  FillerDistance.MatrixMultiplication(MInverse);
    803806  for(int i=0;i<NDIM;i++)
    804807    N[i] = (int) ceil(1./FillerDistance[i]);
    … …  
    835838            }
    836839
    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])                                     ;
     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])                                     ;
    846849          }
    847850
    … …  
    892895            }
    893896      }
    894   delete[](M);
    895   delete[](MInverse);
    896897
    897898  return Filling;

  • src/ellipsoid.cpp

    r0c5eeb r33d774  
    2121#include "tesselation.hpp"
    2222#include "vector.hpp"
     23#include "Matrix.hpp"
    2324#include "verbose.hpp"
    2425
    … …  
    3435  Vector helper, RefPoint;
    3536  double distance = -1.;
    36   double Matrix[NDIM*NDIM];
     37  Matrix Matrix;
    3738  double InverseLength[3];
    3839  double psi,theta,phi; // euler angles in ZX'Z'' convention
    … …  
    5152  theta = EllipsoidAngle[1];
    5253  phi = EllipsoidAngle[2];
    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);
     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);
    6263  helper.MatrixMultiplication(Matrix);
    6364  helper.ScaleAll(InverseLength);
    … …  
    7374  phi = -EllipsoidAngle[2];
    7475  helper.ScaleAll(EllipsoidLength);
    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);
     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);
    8485  helper.MatrixMultiplication(Matrix);
    8586  //Log() << Verbose(4) << "Intersection is at " << helper << "." << endl;

  • src/molecule.cpp

    r0c5eeb r33d774  
    2727#include "tesselation.hpp"
    2828#include "vector.hpp"
     29#include "Matrix.hpp"
    2930#include "World.hpp"
    3031#include "Plane.hpp"
    … …  
    284285  Vector Orthovector1, Orthovector2;  // temporary vectors in coordination construction
    285286  Vector InBondvector;    // vector in direction of *Bond
    286   double *matrix = NULL;
     287  Matrix matrix;
    287288  bond *Binder = NULL;
    288289  double * const cell_size = World::getInstance().getDomain();
    … …  
    307308      } // (signs are correct, was tested!)
    308309    }
    309     matrix = ReturnFullMatrixforSymmetric(cell_size);
     310    double *matrix_double = ReturnFullMatrixforSymmetric(cell_size);
     311    matrix = Matrix(matrix_double);
     312    delete[](matrix_double);
    310313    Orthovector1.MatrixMultiplication(matrix);
    311314    InBondvector -= Orthovector1; // subtract just the additional translation
    312     delete[](matrix);
    313315    bondlength = InBondvector.Norm();
    314316//    Log() << Verbose(4) << "Corrected InBondvector is now: ";
    … …  
    541543      break;
    542544  }
    543   delete[](matrix);
    544545
    545546//  Log() << Verbose(3) << "End of AddHydrogenReplacementAtom." << endl;

  • src/molecule_fragmentation.cpp

    r0c5eeb r33d774  
    2222#include "periodentafel.hpp"
    2323#include "World.hpp"
     24#include "Matrix.hpp"
    2425
    2526/************************************* Functions for class molecule *********************************/
    … …  
    17091710  atom *OtherWalker = NULL;
    17101711  double * const cell_size = World::getInstance().getDomain();
    1711   double *matrix = ReturnFullMatrixforSymmetric(cell_size);
     1712  double *matrix_double = ReturnFullMatrixforSymmetric(cell_size);
     1713  Matrix matrix = Matrix(matrix_double);
     1714  delete[](matrix_double);
    17121715  enum Shading *ColorList = NULL;
    17131716  double tmp;
    … …  
    17801783  delete(AtomStack);
    17811784  delete[](ColorList);
    1782   delete[](matrix);
    17831785  DoLog(2) && (Log() << Verbose(2) << "End of ScanForPeriodicCorrection." << endl);
    17841786};

  • src/molecule_geometry.cpp

    r0c5eeb r33d774  
    1919#include "World.hpp"
    2020#include "Plane.hpp"
     21#include "Matrix.hpp"
    2122#include <boost/foreach.hpp>
    2223
    … …  
    154155
    155156  const double *cell_size = World::getInstance().getDomain();
    156   double *M = ReturnFullMatrixforSymmetric(cell_size);
     157  double *M_double = ReturnFullMatrixforSymmetric(cell_size);
     158  Matrix M = Matrix(M_double);
     159  delete[](M_double);
    157160  a->MatrixMultiplication(M);
    158   delete[](M);
    159161
    160162  return a;
    … …  
    274276{
    275277  double * const cell_size = World::getInstance().getDomain();
    276   double *matrix = ReturnFullMatrixforSymmetric(cell_size);
    277   double *inversematrix = InverseMatrix(matrix);
     278  double *matrix_double = ReturnFullMatrixforSymmetric(cell_size);
     279  Matrix matrix = Matrix(matrix_double);
     280  delete[](matrix_double);
     281  Matrix inversematrix = matrix.invert();
    278282  double tmp;
    279283  bool flag;
    … …  
    324328    }
    325329  } while (!flag);
    326   delete[](matrix);
    327   delete[](inversematrix);
    328330
    329331  Center.Scale(1./static_cast<double>(getAtomCount()));
    … …  
    387389    DoLog(1) && (Log() << Verbose(1) << "Transforming molecule into PAS ... ");
    388390    // the eigenvectors specify the transformation matrix
    389     ActOnAllVectors( &Vector::MatrixMultiplication, (const double *) evec->data );
     391    Matrix M = Matrix(evec->data);
     392    ActOnAllVectors( &Vector::MatrixMultiplication, static_cast<const Matrix>(M));
    390393    DoLog(0) && (Log() << Verbose(0) << "done." << endl);
    391394

  • src/vector.cpp

    r0c5eeb r33d774  
    88
    99#include "vector.hpp"
     10#include "Matrix.hpp"
    1011#include "verbose.hpp"
    1112#include "World.hpp"
    1213#include "Helpers/Assert.hpp"
    1314#include "Helpers/fast_functions.hpp"
     15#include "Exceptions/MathException.hpp"
    1416
    1517#include <iostream>
     18#include <gsl/gsl_blas.h>
     19
    1620
    1721using namespace std;
    … …  
    3438{
    3539  content = gsl_vector_alloc(NDIM);
    36   gsl_vector_set(content,0,src[0]);
    37   gsl_vector_set(content,1,src[1]);
    38   gsl_vector_set(content,2,src[2]);
     40  gsl_vector_memcpy(content, src.content);
    3941}
    4042
    … …  
    4951};
    5052
     53Vector::Vector(gsl_vector *_content) :
     54  content(_content)
     55{}
     56
    5157/**
    5258 * Assignment operator
    … …  
    5561  // check for self assignment
    5662  if(&src!=this){
    57     gsl_vector_set(content,0,src[0]);
    58     gsl_vector_set(content,1,src[1]);
    59     gsl_vector_set(content,2,src[2]);
     63    gsl_vector_memcpy(content, src.content);
    6064  }
    6165  return *this;
    … …  
    101105double Vector::PeriodicDistance(const Vector &y, const double * const cell_size) const
    102106{
    103   double res = distance(y), tmp, matrix[NDIM*NDIM];
     107  double res = distance(y), tmp;
     108  Matrix matrix;
    104109    Vector Shiftedy, TranslationVector;
    105110    int N[NDIM];
    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];
     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];
    115120    // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
    116121    for (N[0]=-1;N[0]<=1;N[0]++)
    … …  
    138143double Vector::PeriodicDistanceSquared(const Vector &y, const double * const cell_size) const
    139144{
    140   double res = DistanceSquared(y), tmp, matrix[NDIM*NDIM];
     145  double res = DistanceSquared(y), tmp;
     146  Matrix matrix;
    141147    Vector Shiftedy, TranslationVector;
    142148    int N[NDIM];
    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];
     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];
    152158    // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
    153159    for (N[0]=-1;N[0]<=1;N[0]++)
    … …  
    172178 * Tries to translate a vector into each adjacent neighbouring cell.
    173179 */
    174 void Vector::KeepPeriodic(const double * const matrix)
    175 {
     180void Vector::KeepPeriodic(const double * const _matrix)
     181{
     182  Matrix matrix = Matrix(_matrix);
    176183  //  int N[NDIM];
    177184  //  bool flag = false;
    … …  
    181188  //  Output(out);
    182189  //  Log() << Verbose(0) << endl;
    183     InverseMatrixMultiplication(matrix);
     190    MatrixMultiplication(matrix.invert());
    184191    for(int i=NDIM;i--;) { // correct periodically
    185192      if (at(i) < 0) {  // get every coefficient into the interval [0,1)
    … …  
    203210{
    204211  double res = 0.;
    205   for (int i=NDIM;i--;)
    206     res += at(i)*y[i];
     212  gsl_blas_ddot(content, y.content, &res);
    207213  return (res);
    208214};
    … …  
    461467};
    462468
     469Vector &Vector::operator*=(const Matrix &mat){
     470  (*this) = mat*(*this);
     471  return *this;
     472}
     473
     474Vector 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
    463481/** Factors given vector \a a times \a m.
    464482 * \param a vector
    … …  
    508526void Vector::Scale(const double factor)
    509527{
    510   for (int i=NDIM;i--;)
    511     at(i) *= factor;
     528  gsl_vector_scale(content,factor);
    512529};
    513530
    … …  
    516533 * \param *Minv inverse matrix
    517534 */
    518 void Vector::WrapPeriodically(const double * const M, const double * const Minv)
    519 {
     535void Vector::WrapPeriodically(const double * const _M, const double * const _Minv)
     536{
     537  Matrix M = Matrix(_M);
     538  Matrix Minv = Matrix(_Minv);
    520539  MatrixMultiplication(Minv);
    521540  // truncate to [0,1] for each axis
    … …  
    550569 * \param *matrix NDIM_NDIM array
    551570 */
    552 void 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;
     571void Vector::MatrixMultiplication(const Matrix &M)
     572{
     573  (*this) *= M;
    560574};
    561575
    … …  
    565579bool Vector::InverseMatrixMultiplication(const double * const A)
    566580{
     581  /*
    567582  double B[NDIM*NDIM];
    568583  double detA = RDET3(A);
    … …  
    586601    return true;
    587602  } 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){
    588612    return false;
    589613  }
    … …  
    679703void Vector::AddVector(const Vector &y)
    680704{
    681   for(int i=NDIM;i--;)
    682     at(i) += y[i];
     705  gsl_vector_add(content, y.content);
    683706}
    684707
    … …  
    688711void Vector::SubtractVector(const Vector &y)
    689712{
    690   for(int i=NDIM;i--;)
    691     at(i) -= y[i];
     713  gsl_vector_sub(content, y.content);
    692714}
    693715

  • src/vector.hpp

    r0c5eeb r33d774  
    2424
    2525class Vector;
     26class Matrix;
    2627
    2728typedef std::vector<Vector> pointset;
    … …  
    3132 */
    3233class Vector : public Space{
     34  friend Vector operator*(const Matrix&,const Vector&);
    3335public:
    34 
    3536  Vector();
    3637  Vector(const double x1, const double x2, const double x3);
    … …  
    5960  void ScaleAll(const double *factor);
    6061  void Scale(const double factor);
    61   void MatrixMultiplication(const double * const M);
     62  void MatrixMultiplication(const Matrix &M);
    6263  bool InverseMatrixMultiplication(const double * const M);
    6364  void KeepPeriodic(const double * const matrix);
    … …  
    9697  Vector const operator+(const Vector& b) const;
    9798  Vector const operator-(const Vector& b) const;
     99  Vector& operator*=(const Matrix&);
    98100
    99101  // Methods inherited from Space
    … …  
    104106
    105107private:
     108  Vector(gsl_vector *);
    106109  gsl_vector *content;
    107110
    … …  
    118121Vector const operator*(const Vector& a, const double m);
    119122Vector const operator*(const double m, const Vector& a);
     123Vector operator*(const Matrix&,const Vector&);
     124
    120125
    121126#endif /*VECTOR_HPP_*/
Note: See TracChangeset for help on using the changeset viewer.