Changeset bbf1bd for src/LinearAlgebra


Ignore:
Timestamp:
Dec 4, 2010, 11:54:32 PM (14 years ago)
Author:
Frederik Heber <heber@…>
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, 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:
6e06bd, b4cf2b
Parents:
5ada94
git-author:
Frederik Heber <heber@…> (11/15/10 20:00:30)
git-committer:
Frederik Heber <heber@…> (12/04/10 23:54:32)
Message:

Extended VectorContent class.

  • VectorContent formerly has been just a structure to contain the gsl_vector due to forward declaration reasons.
  • now VectorContent is a true wrapper to gsl_vector, i.e. all functionality that is now specific to 3 dimensions has been shifted from GSLVector over to VectorContent.
  • VectorContentView is preserved to allow for VectorContent as a view on a row or column of a matrix.
  • changed and renamed unit test gslvectorunittest -> VectorContentUnitTest
  • GSLVector is not used anymore anywhere
  • one long-sough error was a missing assignment operator filled-in in a wrong manner by the compiler for VectorContent.

Note that:

  • gsl_vector is still used at many places
Location:
src/LinearAlgebra
Files:
1 added
9 edited

Legend:

Unmodified
Added
Removed

  • src/LinearAlgebra/Makefile.am

    r5ada94 rbbf1bd  
    2323  vector_ops.cpp \
    2424  Vector.cpp \
     25  VectorContent.cpp \
    2526  VectorInterface.cpp
    2627                           

  • src/LinearAlgebra/MatrixContent.cpp

    r5ada94 rbbf1bd  
    475475Vector operator*(const MatrixContent &mat,const Vector &vec)
    476476{
    477   Vector res;
    478   gsl_blas_dgemv( CblasNoTrans, 1.0, mat.content, vec.content->content, 0.0, res.content->content);
    479   return res;
    480 }
     477  Vector result;
     478  gsl_blas_dgemv( CblasNoTrans, 1.0, mat.content, vec.content->content, 0.0, result.content->content);
     479  return result;
     480}

  • src/LinearAlgebra/MatrixContent.hpp

    r5ada94 rbbf1bd  
    5757  const double at(size_t i, size_t j) const;
    5858  void set(size_t i, size_t j, const double value);
     59
     60  // Initializing
    5961  void setFromDoubleArray(double * x);
     62  void setIdentity();
     63  void setValue(double _value);
     64  void setZero();
    6065
    6166  // Exchanging elements
    … …  
    7681  bool IsPositiveDefinite() const;
    7782  double Determinant() const;
    78 
    79   // Initializing
    80   void setIdentity();
    81   void setZero();
    82   void setValue(double _value);
    8383
    8484  // Operators

  • src/LinearAlgebra/RealSpaceMatrix.cpp

    r5ada94 rbbf1bd  
    8787    VectorContent *rowContent = new VectorViewContent(gsl_matrix_row(content->content,i));
    8888    rows_ptr[i] = new Vector(rowContent);
     89    ASSERT(rowContent == NULL, "RealSpaceMatrix::createViews() - rowContent was not taken over as supposed to happen!");
    8990  }
    9091  // create column views
    … …  
    9293    VectorContent *columnContent = new VectorViewContent(gsl_matrix_column(content->content,i));
    9394    columns_ptr[i] = new Vector(columnContent);
     95    ASSERT(columnContent == NULL, "RealSpaceMatrix::createViews() - columnContent was not taken over as supposed to happen!");
    9496  }
    9597  // create diagonal view
    9698  VectorContent *diagonalContent = new VectorViewContent(gsl_matrix_diagonal(content->content));
    9799  diagonal_ptr = new Vector(diagonalContent);
     100  ASSERT(diagonalContent == NULL, "RealSpaceMatrix::createViews() - diagonalContent was not taken over as supposed to happen!");
    98101}
    99102

  • src/LinearAlgebra/Vector.cpp

    r5ada94 rbbf1bd  
    4141Vector::Vector()
    4242{
    43   content = new VectorContent();
    44 };
    45 
    46 /**
    47  * Copy constructor
    48  */
    49 
     43  content = new VectorContent((size_t) NDIM);
     44};
     45
     46/** Copy constructor.
     47 * \param &src source Vector reference
     48 */
    5049Vector::Vector(const Vector& src)
    5150{
    52   content = new VectorContent();
    53   gsl_vector_memcpy(content->content, src.content->content);
     51  content = new VectorContent(*(src.content));
    5452}
    5553
    5654/** Constructor of class vector.
     55 * \param x1 first component
     56 * \param x2 second component
     57 * \param x3 third component
    5758 */
    5859Vector::Vector(const double x1, const double x2, const double x3)
    5960{
    60   content = new VectorContent();
    61   gsl_vector_set(content->content,0,x1);
    62   gsl_vector_set(content->content,1,x2);
    63   gsl_vector_set(content->content,2,x3);
     61  content = new VectorContent((size_t) NDIM);
     62  content->at(0) = x1;
     63  content->at(1) = x2;
     64  content->at(2) = x3;
    6465};
    6566
    6667/** Constructor of class vector.
     68 * \param x[3] three values to initialize Vector with
    6769 */
    6870Vector::Vector(const double x[3])
    6971{
    70   content = new VectorContent();
    71   gsl_vector_set(content->content,0,x[0]);
    72   gsl_vector_set(content->content,1,x[1]);
    73   gsl_vector_set(content->content,2,x[2]);
    74 };
    75 
    76 Vector::Vector(VectorContent *_content) :
    77   content(_content)
    78 {}
    79 
    80 /**
    81  * Assignment operator
     72  content = new VectorContent((size_t) NDIM);
     73  for (size_t i = NDIM; i--; )
     74    content->at(i) = x[i];
     75};
     76
     77/** Copy constructor of class vector from VectorContent.
     78 * \note This is destructive, i.e. we take over _content.
     79 */
     80Vector::Vector(VectorContent *&_content) :
     81    content(_content)
     82{
     83  _content = NULL;
     84}
     85
     86/** Copy constructor of class vector from VectorContent.
     87 * \note This is non-destructive, i.e. _content is copied.
     88 */
     89Vector::Vector(VectorContent &_content)
     90{
     91  content = new VectorContent(_content);
     92}
     93
     94/** Assignment operator.
     95 * \param &src source vector to assign \a *this to
     96 * \return reference to \a *this
    8297 */
    8398Vector& Vector::operator=(const Vector& src){
    8499  // check for self assignment
    85100  if(&src!=this){
    86     gsl_vector_memcpy(content->content, src.content->content);
     101    *content = *(src.content);
    87102  }
    88103  return *this;
    … …  
    90105
    91106/** Desctructor of class vector.
     107 * Vector::content is deleted.
    92108 */
    93109Vector::~Vector() {

  • src/LinearAlgebra/Vector.hpp

    r5ada94 rbbf1bd  
    3333class Vector : public Space{
    3434  friend Vector operator*(const MatrixContent&,const Vector&);
    35   friend class RealSpaceMatrix;
     35  friend class RealSpaceMatrix; // needs access to Vector(VectorContent*&) due to views
    3636public:
    3737  Vector();
    … …  
    104104
    105105private:
    106   Vector(VectorContent *);
     106  explicit Vector(VectorContent *&);
     107  explicit Vector(VectorContent &_content);
    107108  VectorContent *content;
    108109

  • src/LinearAlgebra/VectorContent.hpp

    r5ada94 rbbf1bd  
    1818 */
    1919
     20#include <iosfwd>
     21
    2022#include <gsl/gsl_vector.h>
    2123
     24class Vector;
     25
     26/** Dummy structure to create a unique signature.
     27 *
     28 */
    2229struct BaseCase{};
    2330
    24 struct VectorContent{
    25   VectorContent(){
    26     content =  gsl_vector_calloc (NDIM);
    27   }
    28   VectorContent(BaseCase){
     31class VectorContent{
     32  friend std::ostream & operator<< (std::ostream& ost, const VectorContent &m);
     33  friend VectorContent const operator*(const VectorContent& a, const double m);
     34  friend VectorContent const operator*(const double m, const VectorContent& a);
     35  friend VectorContent const operator+(const VectorContent& a, const VectorContent& b);
     36  friend VectorContent const operator-(const VectorContent& a, const VectorContent& b);
    2937
    30   }
    31   virtual ~VectorContent(){
    32     if(content){
    33       gsl_vector_free(content);
    34       content = 0;
    35     }
    36   }
     38public:
     39  explicit VectorContent(size_t _dim);
     40  VectorContent(BaseCase);
     41  VectorContent(const VectorContent * const src);
     42  VectorContent(const VectorContent & src);
     43  virtual ~VectorContent();
     44
     45  // Accessing
     46  double &at(size_t m);
     47  const double at(size_t m) const;
     48  double & operator[](size_t i);
     49  const double operator[](size_t i) const;
     50  double *Pointer(size_t m) const;
     51  const double *const_Pointer(size_t m) const;
     52
     53  // Assignment operator
     54  VectorContent &operator=(const VectorContent& src);
     55
     56  // Initializing
     57  void setFromDoubleArray(double * x);
     58  void setFromVector(Vector &v);
     59  void setValue(double x);
     60  void setZero();
     61  int setBasis(size_t m);
     62
     63  // Exchanging elements
     64  int SwapElements(size_t i, size_t j);
     65  int Reverse();
     66
     67  // checking state
     68  bool IsZero() const;
     69  bool IsOne() const;
     70
     71  // properties
     72  bool IsNormalTo(const Vector &normal) const;
     73  bool IsEqualTo(const Vector &a) const;
     74
     75  // properties relative to another VectorContent
     76  double DistanceSquared(const VectorContent &y) const;
     77  double ScalarProduct(const VectorContent &y) const;
     78  double Angle(const VectorContent &y) const;
     79
     80  // operators
     81  bool operator==(const VectorContent& b) const;
     82  const VectorContent& operator+=(const VectorContent& b);
     83  const VectorContent& operator-=(const VectorContent& b);
     84  const VectorContent& operator*=(const double m);
     85
     86  size_t dimension;
    3787  gsl_vector *content;
     88private:
    3889};
    3990
    40 struct VectorViewContent : public VectorContent{
     91std::ostream & operator << (std::ostream& ost, const VectorContent &m);
     92VectorContent const operator*(const VectorContent& a, const double m);
     93VectorContent const operator*(const double m, const VectorContent& a);
     94VectorContent const operator+(const VectorContent& a, const VectorContent& b);
     95VectorContent const operator-(const VectorContent& a, const VectorContent& b);
     96
     97/** Vector view.
     98 * Extension of VectorContent to contain not a gsl_vector but only a view on a
     99 * gsl_vector (or row/column in a gsl_matrix).
     100 *
     101 * We need the above BaseCase here:
     102 * content, i.e. the gsl_vector, must not be allocated, as it is just a view
     103 * with an internal gsl_vector_view. Hence, BaseCase is just dummy class to
     104 * give the constructor a unique signature.
     105 */
     106struct VectorViewContent : public VectorContent
     107{
    41108  VectorViewContent(gsl_vector_view _view) :
    42109    VectorContent(BaseCase()),
    43110    view(_view)
    44111  {
     112    dimension = _view.vector.size;
    45113    content=&view.vector;
    46114  }
    47   virtual ~VectorViewContent(){
     115  ~VectorViewContent(){
    48116    content=0;
    49117  }

  • src/LinearAlgebra/linearsystemofequations.cpp

    r5ada94 rbbf1bd  
    2222#include "defs.hpp"
    2323#include "LinearAlgebra/MatrixContent.hpp"
    24 #include "LinearAlgebra/gslvector.hpp"
     24#include "LinearAlgebra/VectorContent.hpp"
    2525#include "LinearAlgebra/linearsystemofequations.hpp"
     26#include "Helpers/Assert.hpp"
    2627#include "Helpers/logger.hpp"
    2728#include "LinearAlgebra/Vector.hpp"
    2829
    29 #include <cassert>
    3030#include <gsl/gsl_permutation.h>
    3131
    … …  
    3838{
    3939  A = new MatrixContent(m, n);
    40   b = new GSLVector(m);
    41   x = new GSLVector(n);
     40  b = new VectorContent(m);
     41  x = new VectorContent(n);
    4242};
    4343
    … …  
    5858bool LinearSystemOfEquations::SetSymmetric(bool symmetric)
    5959{
    60   assert (rows == columns && "Rows and columns don't have equal size! Setting symmetric not possible.");
     60  ASSERT (rows == columns, "Rows and columns don't have equal size! Setting symmetric not possible.");
    6161  return (IsSymmetric = symmetric);
    6262};
    … …  
    6767void LinearSystemOfEquations::Setb(Vector *x)
    6868{
    69   assert ( columns == NDIM && "Vector class is always three-dimensional, unlike this LEqS!");
    70   b->SetFromVector(*x);
     69  ASSERT  ( columns == NDIM, "Vector class is always three-dimensional, unlike this LEqS!");
     70  b->setFromVector(*x);
    7171};
    7272
    … …  
    7676void LinearSystemOfEquations::Setb(double *x)
    7777{
    78   b->SetFromDoubleArray(x);
     78  b->setFromDoubleArray(x);
    7979};
    8080
    … …  
    9898  // copy solution
    9999  for (size_t i=0;i<x->dimension;i++) {
    100     array[i] = x->Get(i);
     100    array[i] = x->at(i);
    101101  }
    102102  return status;
    … …  
    109109bool LinearSystemOfEquations::GetSolutionAsVector(Vector &v)
    110110{
    111   assert(rows == NDIM && "Solution can only be returned as vector if number of columns is NDIM.");
     111  ASSERT (rows == NDIM, "Solution can only be returned as vector if number of columns is NDIM.");
    112112  bool status = Solve();
    113113
    114114  // copy solution
    115115  for (size_t i=0;i<x->dimension;i++)
    116     v[i] = x->Get(i);
     116    v[i] = x->at(i);
    117117  return status;
    118118};
    … …  
    130130    gsl_permutation * p = gsl_permutation_alloc (x->dimension);
    131131    gsl_linalg_LU_decomp (A->content, p, &s);
    132     gsl_linalg_LU_solve (A->content, p, b->vector, x->vector);
     132    gsl_linalg_LU_solve (A->content, p, b->content, x->content);
    133133    gsl_permutation_free (p);
    134134  } else {  // use Householder
    135135    //MatrixContent *backup = new MatrixContent(rows,columns);
    136136    //*backup = *A;
    137     gsl_linalg_HH_solve (A->content, b->vector, x->vector);
     137    gsl_linalg_HH_solve (A->content, b->content, x->content);
    138138    //*A = *backup;
    139139    //delete(backup);

  • src/LinearAlgebra/linearsystemofequations.hpp

    r5ada94 rbbf1bd  
    2222class Vector;
    2323class MatrixContent;
    24 class GSLVector;
     24class VectorContent;
    2525
    2626/********************************************** declarations *******************************/
    … …  
    4444  bool GetSolutionAsVector(Vector &v);
    4545
    46   GSLVector *b;
    47   GSLVector *x;
     46  VectorContent *b;
     47  VectorContent *x;
    4848  MatrixContent *A;
    4949
Note: See TracChangeset for help on using the changeset viewer.