Changeset 6e06bd for src/LinearAlgebra


Ignore:
Timestamp:
Dec 5, 2010, 12:19:19 AM (15 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, 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:
f03705
Parents:
589112 (diff), bbf1bd (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 'MatrixVectorContentRefactoring' into stable

Location:
src/LinearAlgebra
Files:
3 added
3 deleted
9 edited
1 moved

Legend:

Unmodified
Added
Removed

  • src/LinearAlgebra/Line.cpp

    r589112 r6e06bd  
    2828#include "Helpers/Log.hpp"
    2929#include "Helpers/Verbose.hpp"
    30 #include "LinearAlgebra/gslmatrix.hpp"
     30#include "LinearAlgebra/MatrixContent.hpp"
    3131#include "Helpers/Info.hpp"
    3232#include "Exceptions/LinearDependenceException.hpp"
    … …  
    116116  Vector res;
    117117
    118   auto_ptr<GSLMatrix> M = auto_ptr<GSLMatrix>(new GSLMatrix(4,4));
    119 
    120   M->SetAll(1.);
     118  auto_ptr<MatrixContent> M = auto_ptr<MatrixContent>(new MatrixContent(4,4));
     119
     120  M->setValue(1.);
    121121  for (int i=0;i<3;i++) {
    122     M->Set(0, i, Line1a[i]);
    123     M->Set(1, i, Line1b[i]);
    124     M->Set(2, i, Line2a[i]);
    125     M->Set(3, i, Line2b[i]);
     122    M->set(0, i, Line1a[i]);
     123    M->set(1, i, Line1b[i]);
     124    M->set(2, i, Line2a[i]);
     125    M->set(3, i, Line2b[i]);
    126126  }
    127127

  • src/LinearAlgebra/Makefile.am

    r589112 r6e06bd  
    1212  ${HELPERSOURCE} \
    1313  BoxVector.cpp \
    14   gslmatrix.cpp \
    1514  gslvector.cpp \
    1615  Line.cpp \
    … …  
    1817  LineSegmentSet.cpp \
    1918  linearsystemofequations.cpp \
    20   Matrix.cpp \
     19  MatrixContent.cpp \
    2120  Plane.cpp \
     21  RealSpaceMatrix.cpp \
    2222  Space.cpp \
    2323  vector_ops.cpp \
    2424  Vector.cpp \
     25  VectorContent.cpp \
    2526  VectorInterface.cpp
    2627                           
    2728LINALGHEADER = \
    2829  BoxVector.hpp \
    29   gslmatrix.hpp \
    3030  gslvector.hpp \
    3131  Line.hpp \
    … …  
    3333  LineSegmentSet.hpp \
    3434  linearsystemofequations.hpp \
    35   Matrix.hpp \
     35  MatrixContent.hpp \
     36  RealSpaceMatrix.hpp \
    3637  Plane.hpp \
    3738  Space.hpp \
    3839  vector_ops.hpp \
    3940  Vector.hpp \
     41  VectorContent.hpp \
    4042  VectorInterface.hpp
    4143

  • src/LinearAlgebra/MatrixContent.hpp

    r589112 r6e06bd  
    99#define MATRIXCONTENT_HPP_
    1010
    11 /**
    12  * !file
    13  * The way GSL works does not allow for forward definitions of the structures.
    14  * Because of this the pointer to the gsl_matrix struct is wrapped inside another
    15  * (dumb) object that allows for forward definitions. All methods of this class
    16  * are empty, to allow for a maximum of flexibility.
     11/** MatrixContent is a wrapper for gsl_matrix.
     12 *
     13 * The GNU Scientific Library contaisn some very well written routines for
     14 * linear algebra problems. However, it's syntax is C and hence it does not
     15 * lend itself to readable written code, i.e. C = A * B, where A, B, and C
     16 * are all matrices. Writing code this way is very convenient, concise and
     17 * also in the same manner as one would type in MatLab.
     18 * However, with C++ and its feature to overload functions and its operator
     19 * functions such syntax is easy to create.
     20 *
     21 * Hence, this class is a C++ wrapper to gsl_matrix. There already some other
     22 * libraries, however they fall short for the following reasons:
     23 * -# gslwrap: not very well written and uses floats instead of doubles
     24 * -# gslcpp: last change is from 2007 and only very few commits
     25 * -# o2scl: basically, the same functionality as gsl only written in C++,
     26 *    however it uses GPLv3 license which is inconvenient for MoleCuilder.
     27 *
     28 * <h1>Howto use MatrixContent</h1>
     29 *
     30 * One should not use MatrixContent directly but either have it as a member
     31 * variable in a specialized class or inherit from it.
     32 *
    1733 */
    1834
    1935#include <gsl/gsl_matrix.h>
    2036
    21 struct MatrixContent{
     37class Vector;
     38
     39class MatrixContent
     40{
     41  friend Vector operator*(const MatrixContent &mat,const Vector &vec);
     42  friend class RealSpaceMatrix;
     43  friend class LinearSystemOfEquations;
     44
     45  friend class MatrixContentTest;
     46  friend class MatrixContentSymmetricTest;
     47public:
     48  MatrixContent(size_t rows, size_t columns);
     49  MatrixContent(size_t rows, size_t columns, const double *src);
     50  MatrixContent(gsl_matrix *&src);
     51  MatrixContent(const MatrixContent &src);
     52  MatrixContent(const MatrixContent *src);
     53  ~MatrixContent();
     54
     55  // Accessing
     56  double &at(size_t i, size_t j);
     57  const double at(size_t i, size_t j) const;
     58  void set(size_t i, size_t j, const double value);
     59
     60  // Initializing
     61  void setFromDoubleArray(double * x);
     62  void setIdentity();
     63  void setValue(double _value);
     64  void setZero();
     65
     66  // Exchanging elements
     67  bool SwapRows(size_t i, size_t j);
     68  bool SwapColumns(size_t i, size_t j);
     69  bool SwapRowColumn(size_t i, size_t j);
     70
     71  // Transformations
     72  MatrixContent transpose() const;
     73  MatrixContent &transpose();
     74  gsl_vector* transformToEigenbasis();
     75
     76  // Properties
     77  bool IsNull() const;
     78  bool IsPositive() const;
     79  bool IsNegative() const;
     80  bool IsNonNegative() const;
     81  bool IsPositiveDefinite() const;
     82  double Determinant() const;
     83
     84  // Operators
     85  MatrixContent &operator=(const MatrixContent &src);
     86  const MatrixContent &operator+=(const MatrixContent &rhs);
     87  const MatrixContent &operator-=(const MatrixContent &rhs);
     88  const MatrixContent &operator*=(const MatrixContent &rhs);
     89  const MatrixContent operator*(const MatrixContent &rhs) const;
     90  const MatrixContent &operator*=(const double factor);
     91  bool operator==(const MatrixContent &rhs) const;
     92
     93protected:
     94  double *Pointer(size_t m, size_t n);
     95  const double *const_Pointer(size_t m, size_t n) const;
     96
     97private:
    2298  gsl_matrix * content;
     99  const size_t rows;      // store for internal purposes
     100  const size_t columns;  // store for internal purposes
    23101};
    24102
     103const MatrixContent operator*(const double factor,const MatrixContent& mat);
     104const MatrixContent operator*(const MatrixContent &mat,const double factor);
     105Vector operator*(const MatrixContent &mat,const Vector &vec);
     106
    25107#endif /* MATRIXCONTENT_HPP_ */

  • src/LinearAlgebra/RealSpaceMatrix.hpp

    r589112 r6e06bd  
    11/*
    2  * Matrix.hpp
     2 * RealSpaceMatrix.hpp
    33 *
    44 *  Created on: Jun 25, 2010
    … …  
    66 */
    77
    8 #ifndef MATRIX_HPP_
    9 #define MATRIX_HPP_
     8#ifndef REALSPACEMATRIX_HPP_
     9#define REALSPACEMATRIX_HPP_
    1010
    1111#include <iosfwd>
    … …  
    1919class MatrixContent;
    2020
    21 class Matrix
     21/** 3x3 Matrix class.
     22 * This class has some specific features for 3D space such as rotations or
     23 * hard-coded determinants.
     24 */
     25class RealSpaceMatrix
    2226{
    23   friend Vector operator*(const Matrix&,const Vector&);
     27  friend Vector operator*(const RealSpaceMatrix&,const Vector&);
    2428public:
    25   Matrix();
     29  RealSpaceMatrix();
    2630
    2731  /**
    … …  
    4145   *
    4246   */
    43   Matrix(const double*);
    44   Matrix(const Matrix&);
    45   virtual ~Matrix();
     47  RealSpaceMatrix(const double*);
     48  RealSpaceMatrix(const RealSpaceMatrix&);
     49  RealSpaceMatrix(const MatrixContent&);
     50  virtual ~RealSpaceMatrix();
    4651
    4752  /**
    4853   * Set the matrix to a unit matrix.
    4954   */
    50   void one();
     55  void setIdentity();
    5156
    5257  /**
    5358   * Set all matrix entries to zero.
    5459   */
    55   void zero();
     60  void setZero();
    5661
    5762  /**
    … …  
    6065   * with given angles.
    6166   */
    62   void rotation(const double x, const double y, const double z);
     67  void setRotation(const double x, const double y, const double z);
    6368
    6469  /**
    … …  
    106111   * Rather costly, so use precomputation as often as possible.
    107112   */
    108   Matrix invert() const;
     113  RealSpaceMatrix invert() const;
    109114
    110115  /**
    … …  
    119124   * Calculate the transpose of the matrix.
    120125   */
    121   Matrix transpose() const;
    122   Matrix &transpose();
     126  RealSpaceMatrix transpose() const;
     127  RealSpaceMatrix &transpose();
    123128
    124129  // operators
    125   Matrix &operator=(const Matrix&);
     130  RealSpaceMatrix &operator=(const RealSpaceMatrix&);
    126131
    127   const Matrix &operator+=(const Matrix&);
    128   const Matrix &operator-=(const Matrix&);
    129   const Matrix &operator*=(const Matrix&);
     132  const RealSpaceMatrix &operator+=(const RealSpaceMatrix&);
     133  const RealSpaceMatrix &operator-=(const RealSpaceMatrix&);
     134  const RealSpaceMatrix &operator*=(const RealSpaceMatrix&);
    130135
    131   const Matrix &operator*=(const double);
     136  const RealSpaceMatrix &operator*=(const double);
    132137
    133   const Matrix operator+(const Matrix&) const;
    134   const Matrix operator-(const Matrix&) const;
    135   const Matrix operator*(const Matrix&) const;
     138  const RealSpaceMatrix operator+(const RealSpaceMatrix&) const;
     139  const RealSpaceMatrix operator-(const RealSpaceMatrix&) const;
     140  const RealSpaceMatrix operator*(const RealSpaceMatrix&) const;
    136141
    137   bool operator==(const Matrix&) const;
     142  bool operator==(const RealSpaceMatrix&) const;
    138143
    139144private:
    140   Matrix(MatrixContent*);
     145  RealSpaceMatrix(MatrixContent*);
    141146  void createViews();
    142147  MatrixContent *content;
    … …  
    147152};
    148153
    149 const Matrix operator*(const double,const Matrix&);
    150 const Matrix operator*(const Matrix&,const double);
     154const RealSpaceMatrix operator*(const double,const RealSpaceMatrix&);
     155const RealSpaceMatrix operator*(const RealSpaceMatrix&,const double);
    151156
    152157/**
    … …  
    164169 * 5 -> (2,2)
    165170 */
    166 Matrix ReturnFullMatrixforSymmetric(const double * const cell_size);
     171RealSpaceMatrix ReturnFullMatrixforSymmetric(const double * const cell_size);
    167172
    168 std::ostream &operator<<(std::ostream&,const Matrix&);
    169 Vector operator*(const Matrix&,const Vector&);
    170 Vector& operator*=(Vector&,const Matrix&);
     173std::ostream &operator<<(std::ostream&,const RealSpaceMatrix&);
     174Vector operator*(const RealSpaceMatrix&,const Vector&);
     175Vector& operator*=(Vector&,const RealSpaceMatrix&);
    171176
    172 #endif /* MATRIX_HPP_ */
     177#endif /* REALSPACEMATRIX_HPP_ */

  • src/LinearAlgebra/Vector.cpp

    r589112 r6e06bd  
    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

    r589112 r6e06bd  
    2222
    2323class Vector;
    24 class Matrix;
     24class RealSpaceMatrix;
     25class MatrixContent;
    2526struct VectorContent;
    2627
    … …  
    3132 */
    3233class Vector : public Space{
    33   friend Vector operator*(const Matrix&,const Vector&);
    34   friend class Matrix;
     34  friend Vector operator*(const MatrixContent&,const Vector&);
     35  friend class RealSpaceMatrix; // needs access to Vector(VectorContent*&) due to views
    3536public:
    3637  Vector();
    … …  
    103104
    104105private:
    105   Vector(VectorContent *);
     106  explicit Vector(VectorContent *&);
     107  explicit Vector(VectorContent &_content);
    106108  VectorContent *content;
    107109

  • src/LinearAlgebra/VectorContent.hpp

    r589112 r6e06bd  
    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

    r589112 r6e06bd  
    2121
    2222#include "defs.hpp"
    23 #include "LinearAlgebra/gslmatrix.hpp"
    24 #include "LinearAlgebra/gslvector.hpp"
     23#include "LinearAlgebra/MatrixContent.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
    … …  
    3737LinearSystemOfEquations::LinearSystemOfEquations(int m, int n) : rows(m), columns(n), IsSymmetric(false)
    3838{
    39   A = new GSLMatrix(m, n);
    40   b = new GSLVector(m);
    41   x = new GSLVector(n);
     39  A = new MatrixContent(m, 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
    … …  
    8585void LinearSystemOfEquations::SetA(double *x)
    8686{
    87   A->SetFromDoubleArray(x);
     87  A->setFromDoubleArray(x);
    8888};
    8989
    … …  
    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};
    … …  
    129129  if (IsSymmetric) { // use LU
    130130    gsl_permutation * p = gsl_permutation_alloc (x->dimension);
    131     gsl_linalg_LU_decomp (A->matrix, p, &s);
    132     gsl_linalg_LU_solve (A->matrix, p, b->vector, x->vector);
     131    gsl_linalg_LU_decomp (A->content, p, &s);
     132    gsl_linalg_LU_solve (A->content, p, b->content, x->content);
    133133    gsl_permutation_free (p);
    134134  } else {  // use Householder
    135     //GSLMatrix *backup = new GSLMatrix(rows,columns);
     135    //MatrixContent *backup = new MatrixContent(rows,columns);
    136136    //*backup = *A;
    137     gsl_linalg_HH_solve (A->matrix, 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

    r589112 r6e06bd  
    2121
    2222class Vector;
    23 class GSLMatrix;
    24 class GSLVector;
     23class MatrixContent;
     24class VectorContent;
    2525
    2626/********************************************** declarations *******************************/
    … …  
    4444  bool GetSolutionAsVector(Vector &v);
    4545
    46   GSLVector *b;
    47   GSLVector *x;
    48   GSLMatrix *A;
     46  VectorContent *b;
     47  VectorContent *x;
     48  MatrixContent *A;
    4949
    5050private:

  • src/LinearAlgebra/vector_ops.cpp

    r589112 r6e06bd  
    2424#include "Helpers/Log.hpp"
    2525#include "Helpers/Verbose.hpp"
    26 #include "LinearAlgebra/gslmatrix.hpp"
    2726#include "leastsquaremin.hpp"
    2827#include "Helpers/Info.hpp"
Note: See TracChangeset for help on using the changeset viewer.