Changeset 3bc926


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:
cca9ef
Parents:
9eb7580
git-author:
Frederik Heber <heber@…> (11/15/10 12:36:18)
git-committer:
Frederik Heber <heber@…> (12/04/10 23:54:32)
Message:

Extended MatrixContent class.

  • MatrixContent formerly has been just a structure to contain the gsl_matrix due to forward declaration reasons.
  • now MatrixContent is a true wrapper to gsl_matrix, i.e. all functionality that is now specific to 3 dimensions has been shifted from Matrix over to MatrixContent.
Location:
src
Files:
1 added
8 edited

Legend:

Unmodified
Added
Removed

  • src/LinearAlgebra/Makefile.am

    r9eb7580 r3bc926  
    1919  linearsystemofequations.cpp \
    2020  Matrix.cpp \
     21  MatrixContent.cpp \
    2122  Plane.cpp \
    2223  Space.cpp \
    … …  
    3435  linearsystemofequations.hpp \
    3536  Matrix.hpp \
     37  MatrixContent.hpp \
    3638  Plane.hpp \
    3739  Space.hpp \
    3840  vector_ops.hpp \
    3941  Vector.hpp \
     42  VectorContent.hpp \
    4043  VectorInterface.hpp
    4144

  • src/LinearAlgebra/Matrix.cpp

    r9eb7580 r3bc926  
    2626#include "Helpers/Assert.hpp"
    2727#include "LinearAlgebra/Vector.hpp"
    28 #include "VectorContent.hpp"
    29 #include "MatrixContent.hpp"
     28#include "LinearAlgebra/VectorContent.hpp"
     29#include "LinearAlgebra/MatrixContent.hpp"
    3030
    3131#include <gsl/gsl_blas.h>
    … …  
    4141Matrix::Matrix()
    4242{
    43   content = new MatrixContent();
    44   content->content = gsl_matrix_calloc(NDIM, NDIM);
    45   createViews();
    46 }
    47 
    48 Matrix::Matrix(const double* src){
    49   content = new MatrixContent();
    50   content->content = gsl_matrix_alloc(NDIM, NDIM);
    51   set(0,0, src[0]);
    52   set(1,0, src[1]);
    53   set(2,0, src[2]);
    54 
    55   set(0,1, src[3]);
    56   set(1,1, src[4]);
    57   set(2,1, src[5]);
    58 
    59   set(0,2, src[6]);
    60   set(1,2, src[7]);
    61   set(2,2, src[8]);
    62   createViews();
    63 }
    64 
    65 Matrix::Matrix(const Matrix &src){
    66   content = new MatrixContent();
    67   content->content = gsl_matrix_alloc(NDIM, NDIM);
    68   gsl_matrix_memcpy(content->content,src.content->content);
    69   createViews();
    70 }
    71 
    72 Matrix::Matrix(MatrixContent* _content) :
    73   content(_content)
    74 {
     43  content = new MatrixContent(NDIM, NDIM);
     44  createViews();
     45}
     46
     47Matrix::Matrix(const double* src)
     48{
     49  content = new MatrixContent(NDIM, NDIM, src);
     50  createViews();
     51}
     52
     53Matrix::Matrix(const Matrix &src)
     54{
     55  content = new MatrixContent(src.content);
     56  createViews();
     57}
     58
     59Matrix::Matrix(const MatrixContent &src)
     60{
     61  content = new MatrixContent(src);
     62  createViews();
     63}
     64
     65Matrix::Matrix(MatrixContent* _content)
     66{
     67  content = new MatrixContent(_content);
    7568  createViews();
    7669}
    … …  
    8679  }
    8780  delete diagonal_ptr;
    88   gsl_matrix_free(content->content);
    8981  delete content;
    9082}
    … …  
    10799
    108100void Matrix::setIdentity(){
    109   for(int i=NDIM;i--;){
    110     for(int j=NDIM;j--;){
    111       set(i,j,i==j);
    112     }
    113   }
     101  content->setIdentity();
    114102}
    115103
    116104void Matrix::setZero(){
    117   for(int i=NDIM;i--;){
    118     for(int j=NDIM;j--;){
    119       set(i,j,0.);
    120     }
    121   }
     105  content->setZero();
    122106}
    123107
    … …  
    135119}
    136120
    137 Matrix &Matrix::operator=(const Matrix &src){
    138   if(&src!=this){
    139     gsl_matrix_memcpy(content->content,src.content->content);
    140   }
    141   return *this;
    142 }
    143 
    144 const Matrix &Matrix::operator+=(const Matrix &rhs){
    145   gsl_matrix_add(content->content, rhs.content->content);
    146   return *this;
    147 }
    148 
    149 const Matrix &Matrix::operator-=(const Matrix &rhs){
    150   gsl_matrix_sub(content->content, rhs.content->content);
    151   return *this;
    152 }
    153 
    154 const Matrix &Matrix::operator*=(const Matrix &rhs){
    155   (*this) = (*this)*rhs;
    156   return *this;
    157 }
    158 
    159 const Matrix Matrix::operator+(const Matrix &rhs) const{
     121Matrix &Matrix::operator=(const Matrix &src)
     122{
     123  // MatrixContent checks for self-assignment
     124  *content = *(src.content);
     125  return *this;
     126}
     127
     128const Matrix &Matrix::operator+=(const Matrix &rhs)
     129{
     130  *content += *(rhs.content);
     131  return *this;
     132}
     133
     134const Matrix &Matrix::operator-=(const Matrix &rhs)
     135{
     136  *content -= *(rhs.content);
     137  return *this;
     138}
     139
     140const Matrix &Matrix::operator*=(const Matrix &rhs)
     141{
     142  (*content) *= (*rhs.content);
     143  return *this;
     144}
     145
     146const Matrix Matrix::operator+(const Matrix &rhs) const
     147{
    160148  Matrix tmp = *this;
    161149  tmp+=rhs;
    … …  
    163151}
    164152
    165 const Matrix Matrix::operator-(const Matrix &rhs) const{
     153const Matrix Matrix::operator-(const Matrix &rhs) const
     154{
    166155  Matrix tmp = *this;
    167156  tmp-=rhs;
    … …  
    169158}
    170159
    171 const Matrix Matrix::operator*(const Matrix &rhs) const{
    172   gsl_matrix *res = gsl_matrix_alloc(NDIM, NDIM);
    173   gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, content->content, rhs.content->content, 0.0, res);
    174   MatrixContent *content= new MatrixContent();
    175   content->content = res;
    176   return Matrix(content);
    177 }
    178 
    179 double &Matrix::at(size_t i, size_t j){
    180   ASSERT(i>=0&&i<NDIM,"Index i for Matrix access out of range");
    181   ASSERT(j>=0&&j<NDIM,"Index j for Matrix access out of range");
    182   return *gsl_matrix_ptr (content->content, i, j);
    183 }
    184 
    185 const double Matrix::at(size_t i, size_t j) const{
    186   ASSERT(i>=0&&i<NDIM,"Index i for Matrix access out of range");
    187   ASSERT(j>=0&&j<NDIM,"Index j for Matrix access out of range");
    188   return gsl_matrix_get(content->content, i, j);
    189 }
    190 
    191 Vector &Matrix::row(size_t i){
     160const Matrix Matrix::operator*(const Matrix &rhs) const
     161{
     162  Matrix tmp(content);
     163  tmp *= rhs;
     164  return tmp;
     165}
     166
     167double &Matrix::at(size_t i, size_t j)
     168{
     169  return content->at(i,j);
     170}
     171
     172const double Matrix::at(size_t i, size_t j) const
     173{
     174  return content->at(i,j);
     175}
     176
     177Vector &Matrix::row(size_t i)
     178{
    192179  ASSERT(i>=0&&i<NDIM,"Index i for Matrix access out of range");
    193180  return *rows_ptr[i];
    194181}
    195182
    196 const Vector &Matrix::row(size_t i) const{
     183const Vector &Matrix::row(size_t i) const
     184{
    197185  ASSERT(i>=0&&i<NDIM,"Index i for Matrix access out of range");
    198186  return *rows_ptr[i];
    199187}
    200188
    201 Vector &Matrix::column(size_t i){
     189Vector &Matrix::column(size_t i)
     190{
    202191  ASSERT(i>=0&&i<NDIM,"Index i for Matrix access out of range");
    203192  return *columns_ptr[i];
    204193}
    205194
    206 const Vector &Matrix::column(size_t i) const{
     195const Vector &Matrix::column(size_t i) const
     196{
    207197  ASSERT(i>=0&&i<NDIM,"Index i for Matrix access out of range");
    208198  return *columns_ptr[i];
    209199}
    210200
    211 Vector &Matrix::diagonal(){
     201Vector &Matrix::diagonal()
     202{
    212203  return *diagonal_ptr;
    213204}
    214205
    215 const Vector &Matrix::diagonal() const{
     206const Vector &Matrix::diagonal() const
     207{
    216208  return *diagonal_ptr;
    217209}
    218210
    219 void Matrix::set(size_t i, size_t j, const double value){
    220   ASSERT(i>=0&&i<NDIM,"Index i for Matrix access out of range");
    221   ASSERT(j>=0&&j<NDIM,"Index j for Matrix access out of range");
    222   gsl_matrix_set(content->content,i,j,value);
     211void Matrix::set(size_t i, size_t j, const double value)
     212{
     213  content->set(i,j, value);
    223214}
    224215
    … …  
    253244}
    254245
    255 Matrix Matrix::transpose() const{
    256   MatrixContent *newContent = new MatrixContent();
    257   newContent->content = gsl_matrix_alloc(NDIM, NDIM);
    258   gsl_matrix_transpose_memcpy(newContent->content, content->content);
    259   Matrix res = Matrix(newContent);
     246Matrix Matrix::transpose() const
     247{
     248  std::cout << "const Matrix::transpose()." << std::endl;
     249  Matrix res = Matrix(const_cast<const MatrixContent *>(content)->transpose());
    260250  return res;
    261251}
    … …  
    263253Matrix &Matrix::transpose()
    264254{
    265   double tmp;
    266   for (int i=0;i<NDIM;i++)
    267     for (int j=i+1;j<NDIM;j++) {
    268       tmp = at(j,i);
    269       at(j,i) = at(i,j);
    270       at(i,j) = tmp;
    271     }
    272   return *this;
    273 }
    274 
     255  std::cout << "Matrix::transpose()." << std::endl;
     256  content->transpose();
     257  return *this;
     258}
    275259
    276260Vector Matrix::transformToEigenbasis()
    277261{
    278   gsl_eigen_symmv_workspace *T = gsl_eigen_symmv_alloc(NDIM);
    279   gsl_vector *eval = gsl_vector_alloc(NDIM);
    280   gsl_matrix *evec = gsl_matrix_alloc(NDIM, NDIM);
    281   gsl_eigen_symmv(content->content, eval, evec, T);
    282   gsl_eigen_symmv_free(T);
    283   gsl_matrix_memcpy(content->content, evec);
    284   gsl_matrix_free(evec);
     262  gsl_vector *eval = content->transformToEigenbasis();
    285263  Vector evalues(gsl_vector_get(eval,0), gsl_vector_get(eval,1), gsl_vector_get(eval,2));
    286264  gsl_vector_free(eval);
    … …  
    288266}
    289267
    290 const Matrix &Matrix::operator*=(const double factor){
    291   gsl_matrix_scale(content->content, factor);
    292   return *this;
    293 }
    294 
    295 const Matrix operator*(const double factor,const Matrix& mat){
     268const Matrix &Matrix::operator*=(const double factor)
     269    {
     270  *content *= factor;
     271  return *this;
     272}
     273
     274const Matrix operator*(const double factor,const Matrix& mat)
     275{
    296276  Matrix tmp = mat;
    297   tmp*=factor;
    298   return tmp;
    299 }
    300 
    301 const Matrix operator*(const Matrix &mat,const double factor){
     277  return tmp *= factor;
     278}
     279
     280const Matrix operator*(const Matrix &mat,const double factor)
     281{
    302282  return factor*mat;
    303283}
    304284
    305 bool Matrix::operator==(const Matrix &rhs) const{
    306   for(int i=NDIM;i--;){
    307     for(int j=NDIM;j--;){
    308       if(fabs(at(i,j)-rhs.at(i,j))>MYEPSILON){
    309         return false;
    310       }
    311     }
    312   }
    313   return true;
     285bool Matrix::operator==(const Matrix &rhs) const
     286{
     287  return (*content == *(rhs.content));
    314288}
    315289
    … …  
    333307};
    334308
    335 ostream &operator<<(ostream &ost,const Matrix &mat){
     309ostream &operator<<(ostream &ost,const Matrix &mat)
     310{
    336311  for(int i = 0;i<NDIM;++i){
    337312    ost << "\n";
    … …  
    345320}
    346321
    347 Vector operator*(const Matrix &mat,const Vector &vec){
    348   Vector res;
    349   gsl_blas_dgemv( CblasNoTrans, 1.0, mat.content->content, vec.content->content, 0.0, res.content->content);
    350   return res;
    351 }
    352 
    353 Vector &operator*=(Vector& lhs,const Matrix &mat){
     322Vector operator*(const Matrix &mat,const Vector &vec)
     323{
     324  return (*mat.content) * vec;
     325}
     326
     327Vector &operator*=(Vector& lhs,const Matrix &mat)
     328{
    354329  lhs = mat*lhs;
    355330  return lhs;

  • src/LinearAlgebra/Matrix.hpp

    r9eb7580 r3bc926  
    4343  Matrix(const double*);
    4444  Matrix(const Matrix&);
     45  Matrix(const MatrixContent&);
    4546  virtual ~Matrix();
    4647

  • src/LinearAlgebra/MatrixContent.hpp

    r9eb7580 r3bc926  
    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 Matrix;
     43public:
     44  MatrixContent(size_t rows, size_t columns);
     45  MatrixContent(size_t rows, size_t columns, const double *src);
     46  MatrixContent(gsl_matrix *&src);
     47  MatrixContent(const MatrixContent &src);
     48  MatrixContent(const MatrixContent *src);
     49  ~MatrixContent();
     50
     51  // Accessing
     52  double &at(size_t i, size_t j);
     53  const double at(size_t i, size_t j) const;
     54  void set(size_t i, size_t j, const double value);
     55
     56  // Transformations
     57  MatrixContent transpose() const;
     58  MatrixContent &transpose();
     59  gsl_vector* transformToEigenbasis();
     60
     61  // Initializing
     62  void setIdentity();
     63  void setZero();
     64
     65  // Operators
     66  MatrixContent &operator=(const MatrixContent &src);
     67  const MatrixContent &operator+=(const MatrixContent &rhs);
     68  const MatrixContent &operator-=(const MatrixContent &rhs);
     69  const MatrixContent &operator*=(const MatrixContent &rhs);
     70  const MatrixContent operator*(const MatrixContent &rhs) const;
     71  const MatrixContent &operator*=(const double factor);
     72  bool operator==(const MatrixContent &rhs) const;
     73
     74
     75private:
    2276  gsl_matrix * content;
     77  const size_t rows;      // store for internal purposes
     78  const size_t columns;  // store for internal purposes
    2379};
    2480
     81const MatrixContent operator*(const double factor,const MatrixContent& mat);
     82const MatrixContent operator*(const MatrixContent &mat,const double factor);
     83Vector operator*(const MatrixContent &mat,const Vector &vec);
     84
    2585#endif /* MATRIXCONTENT_HPP_ */

  • src/LinearAlgebra/Vector.hpp

    r9eb7580 r3bc926  
    2323class Vector;
    2424class Matrix;
     25class MatrixContent;
    2526struct VectorContent;
    2627
    … …  
    3132 */
    3233class Vector : public Space{
    33   friend Vector operator*(const Matrix&,const Vector&);
     34  friend Vector operator*(const MatrixContent&,const Vector&);
    3435  friend class Matrix;
    3536public:

  • src/LinearAlgebra/gslmatrix.cpp

    r9eb7580 r3bc926  
    2727
    2828#include <cassert>
     29#include <iostream>
    2930#include <gsl/gsl_linalg.h>
     31#include <gsl/gsl_blas.h>
     32#include <gsl/gsl_eigen.h>
     33#include <gsl/gsl_matrix.h>
     34#include <gsl/gsl_vector.h>
    3035
    3136/** Constructor of class GSLMatrix.
    … …  
    4651  matrix = gsl_matrix_alloc(rows, columns);
    4752  gsl_matrix_memcpy (matrix, src->matrix);
     53};
     54
     55/** Copy constructor of class GSLMatrix.
     56 * We take over pointer to gsl_matrix and set the parameter pointer to NULL
     57 * afterwards.
     58 * \param m row count
     59 * \param n column count
     60 * \param *&src source gsl_matrix
     61 */
     62GSLMatrix::GSLMatrix(size_t m, size_t n, gsl_matrix *&src) :
     63    rows(m),
     64    columns(n)
     65{
     66  matrix = src;
     67  src = NULL;
    4868};
    4969
    … …  
    90110 * \return (m,n)-th element of matrix
    91111 */
    92 double GSLMatrix::Get(size_t m, size_t n)
     112double GSLMatrix::Get(size_t m, size_t n) const
    93113{
    94114  return gsl_matrix_get (matrix, m, n);
    … …  
    121141 * \return const pointer to \a m-th element
    122142 */
    123 const double *GSLMatrix::const_Pointer(size_t m, size_t n)
     143const double *GSLMatrix::const_Pointer(size_t m, size_t n) const
    124144{
    125145  return gsl_matrix_const_ptr (matrix, m, n);
    126146};
     147
     148/** Returns the number of rows of the matrix.
     149 * \return number of rows
     150 */
     151size_t GSLMatrix::getRowCount() const
     152{
     153  return rows;
     154}
     155
     156/** Returns the number of columns of the matrix.
     157 * \return number of columns
     158 */
     159size_t GSLMatrix::getColumnCount() const
     160{
     161  return columns;
     162}
    127163
    128164/* ========================== Initializing =============================== */
    … …  
    204240 * \return true - is null, false - else
    205241 */
    206 bool GSLMatrix::IsNull()
     242bool GSLMatrix::IsNull() const
    207243{
    208244  return gsl_matrix_isnull (matrix);
    … …  
    212248 * \return true - is positive, false - else
    213249 */
    214 bool GSLMatrix::IsPositive()
     250bool GSLMatrix::IsPositive() const
    215251{
    216252  return gsl_matrix_ispos (matrix);
    … …  
    220256 * \return true - is negative, false - else
    221257 */
    222 bool GSLMatrix::IsNegative()
     258bool GSLMatrix::IsNegative() const
    223259{
    224260  return gsl_matrix_isneg (matrix);
    … …  
    228264 * \return true - is non-negative, false - else
    229265 */
    230 bool GSLMatrix::IsNonNegative()
     266bool GSLMatrix::IsNonNegative() const
    231267{
    232268  return gsl_matrix_isnonneg (matrix);
    … …  
    237273 * \return true - matrix is positive-definite, false - else
    238274 */
    239 bool GSLMatrix::IsPositiveDefinite()
     275bool GSLMatrix::IsPositiveDefinite() const
    240276{
    241277  if (rows != columns)  // only possible for square matrices.
    … …  
    249285 * if matrix is square, uses LU decomposition.
    250286 */
    251 double GSLMatrix::Determinant()
     287double GSLMatrix::Determinant() const
    252288{
    253289  int signum = 0;
    … …  
    259295};
    260296
     297
     298/* ============================ Properties ============================== */
     299
     300const GSLMatrix &GSLMatrix::operator+=(const GSLMatrix &rhs)
     301{
     302  gsl_matrix_add(matrix, rhs.matrix);
     303  return *this;
     304}
     305
     306const GSLMatrix &GSLMatrix::operator-=(const GSLMatrix &rhs)
     307{
     308  gsl_matrix_sub(matrix, rhs.matrix);
     309  return *this;
     310}
     311
     312const GSLMatrix &GSLMatrix::operator*=(const GSLMatrix &rhs)
     313{
     314  (*this) = (*this)*rhs;
     315  return *this;
     316}
     317
     318const GSLMatrix GSLMatrix::operator+(const GSLMatrix &rhs) const
     319{
     320  GSLMatrix tmp = *this;
     321  tmp+=rhs;
     322  return tmp;
     323}
     324
     325const GSLMatrix GSLMatrix::operator-(const GSLMatrix &rhs) const
     326{
     327  GSLMatrix tmp = *this;
     328  tmp-=rhs;
     329  return tmp;
     330}
     331
     332const GSLMatrix GSLMatrix::operator*(const GSLMatrix &rhs) const
     333{
     334  gsl_matrix *res = gsl_matrix_alloc(rhs.rows, rhs.columns);
     335  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, matrix, rhs.matrix, 0.0, res);
     336  return GSLMatrix(rhs.rows, rhs.columns, res);
     337}
     338
     339/** Print matrix in a matlab style manner.
     340 * \param &ost reference to output stream
     341 * \param &mat reference to matrix to print
     342 * \return reference to obtained output stream for concatenation
     343 */
     344ostream &operator<<(ostream &ost,const GSLMatrix &mat)
     345{
     346  ost << "[";
     347  for(size_t i = 0;i<mat.rows;++i){
     348    for(size_t j = 0; j<mat.columns;++j){
     349      ost << mat.Get(i,j);
     350      if(j!=mat.columns-1)
     351        ost << " ";
     352    }
     353    if(i!=mat.rows-1)
     354      ost << "; ";
     355  }
     356  ost << "]";
     357  return ost;
     358}
     359

  • src/LinearAlgebra/gslmatrix.hpp

    r9eb7580 r3bc926  
    99#define GSLMATRIX_HPP_
    1010
    11 using namespace std;
    12 
    1311/*********************************************** includes ***********************************/
    1412
    … …  
    1816#endif
    1917
     18#include <iosfwd>
     19
    2020#include <gsl/gsl_matrix.h>
     21
    2122
    2223/****************************************** forward declarations *****************************/
    … …  
    2829class GSLMatrix {
    2930  friend class LinearSystemOfEquations;
     31  friend std::ostream &operator<<(std::ostream &ost,const GSLMatrix &mat);
    3032
    3133public:
    3234  GSLMatrix(size_t m, size_t n);
    3335  GSLMatrix(const GSLMatrix * const src);
     36  GSLMatrix(size_t m, size_t n, gsl_matrix *&src);
    3437  ~GSLMatrix();
    3538
    3639  // Accessing
    3740  void SetFromDoubleArray(double * x);
    38   double Get(size_t m, size_t n);
     41  double Get(size_t m, size_t n) const;
    3942  void Set(size_t m, size_t n, double x);
    4043  double *Pointer(size_t m, size_t n);
    41   const double *const_Pointer(size_t m, size_t n);
     44  const double *const_Pointer(size_t m, size_t n) const;
     45  size_t getRowCount() const;
     46  size_t getColumnCount() const;
    4247
    4348  // Initializing
    … …  
    5055  bool SwapColumns(size_t i, size_t j);
    5156  bool SwapRowColumn(size_t i, size_t j);
     57
     58  // Properties
    5259  bool Transpose();
    53   bool IsNull();
    54   bool IsPositive();
    55   bool IsNegative();
    56   bool IsNonNegative();
    57   bool IsPositiveDefinite();
    58   double Determinant();
     60  bool IsNull() const;
     61  bool IsPositive() const;
     62  bool IsNegative() const;
     63  bool IsNonNegative() const;
     64  bool IsPositiveDefinite() const;
     65  double Determinant() const;
    5966
    60   GSLMatrix& operator=(const GSLMatrix& rhs);
     67  // Operators
     68  GSLMatrix &operator=(const GSLMatrix &src);
     69  const GSLMatrix &operator+=(const GSLMatrix &rhs);
     70  const GSLMatrix &operator-=(const GSLMatrix &rhs);
     71  const GSLMatrix &operator*=(const GSLMatrix &rhs);
     72  const GSLMatrix operator+(const GSLMatrix &rhs) const;
     73  const GSLMatrix operator-(const GSLMatrix &rhs) const;
     74  const GSLMatrix operator*(const GSLMatrix &rhs) const;
    6175
    6276private:
    … …  
    6781};
    6882
     83std::ostream &operator<<(std::ostream &ost,const GSLMatrix &mat);
    6984
    7085#endif /* GSLMATRIX_HPP_ */

  • src/unittests/MatrixUnittest.cpp

    r9eb7580 r3bc926  
    3838void MatrixUnittest::setUp(){
    3939  zero = new Matrix();
     40  for(int i =NDIM;i--;) {
     41    for(int j =NDIM;j--;) {
     42      zero->at(i,j)=0.;
     43    }
     44  }
    4045  one = new Matrix();
    4146  for(int i =NDIM;i--;){
    … …  
    172177  // transpose of unit is unit
    173178  res.setIdentity();
    174   (const Matrix)res.transpose();
     179  res.transpose();
    175180  CPPUNIT_ASSERT_EQUAL(res,*one);
    176181
    … …  
    185190
    186191  res =(*zero) *(*zero);
     192  std::cout << *zero << " times " << *zero << " is " << res << std::endl;
    187193  CPPUNIT_ASSERT_EQUAL(res,*zero);
    188194  res =(*zero) *(*one);
Note: See TracChangeset for help on using the changeset viewer.