Changeset 3bc926 for src/LinearAlgebra/Matrix.cpp

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.

File:

• src/LinearAlgebra/Matrix.cpp (modified) (12 diffs)

src/LinearAlgebra/Matrix.cpp

@@ -26 +26 @@
 #include "Helpers/Assert.hpp"
 #include "LinearAlgebra/Vector.hpp"
-#include "VectorContent.hpp"
-#include "MatrixContent.hpp"
+#include "LinearAlgebra/VectorContent.hpp"
+#include "LinearAlgebra/MatrixContent.hpp"
 
 #include <gsl/gsl_blas.h>
@@ -41 +41 @@
 Matrix::Matrix()
 {
-  content = new MatrixContent();
-  content->content = gsl_matrix_calloc(NDIM, NDIM);
-  createViews();
-}
-
-Matrix::Matrix(const double* src){
-  content = new MatrixContent();
-  content->content = gsl_matrix_alloc(NDIM, NDIM);
-  set(0,0, src[0]);
-  set(1,0, src[1]);
-  set(2,0, src[2]);
-
-  set(0,1, src[3]);
-  set(1,1, src[4]);
-  set(2,1, src[5]);
-
-  set(0,2, src[6]);
-  set(1,2, src[7]);
-  set(2,2, src[8]);
-  createViews();
-}
-
-Matrix::Matrix(const Matrix &src){
-  content = new MatrixContent();
-  content->content = gsl_matrix_alloc(NDIM, NDIM);
-  gsl_matrix_memcpy(content->content,src.content->content);
-  createViews();
-}
-
-Matrix::Matrix(MatrixContent* _content) :
-  content(_content)
-{
+  content = new MatrixContent(NDIM, NDIM);
+  createViews();
+}
+
+Matrix::Matrix(const double* src)
+{
+  content = new MatrixContent(NDIM, NDIM, src);
+  createViews();
+}
+
+Matrix::Matrix(const Matrix &src)
+{
+  content = new MatrixContent(src.content);
+  createViews();
+}
+
+Matrix::Matrix(const MatrixContent &src)
+{
+  content = new MatrixContent(src);
+  createViews();
+}
+
+Matrix::Matrix(MatrixContent* _content)
+{
+  content = new MatrixContent(_content);
   createViews();
 }
@@ -86 +79 @@
   }
   delete diagonal_ptr;
-  gsl_matrix_free(content->content);
   delete content;
 }
@@ -107 +99 @@
 
 void Matrix::setIdentity(){
-  for(int i=NDIM;i--;){
-    for(int j=NDIM;j--;){
-      set(i,j,i==j);
-    }
-  }
+  content->setIdentity();
 }
 
 void Matrix::setZero(){
-  for(int i=NDIM;i--;){
-    for(int j=NDIM;j--;){
-      set(i,j,0.);
-    }
-  }
+  content->setZero();
 }
 
@@ -135 +119 @@
 }
 
-Matrix &Matrix::operator=(const Matrix &src){
-  if(&src!=this){
-    gsl_matrix_memcpy(content->content,src.content->content);
-  }
-  return *this;
-}
-
-const Matrix &Matrix::operator+=(const Matrix &rhs){
-  gsl_matrix_add(content->content, rhs.content->content);
-  return *this;
-}
-
-const Matrix &Matrix::operator-=(const Matrix &rhs){
-  gsl_matrix_sub(content->content, rhs.content->content);
-  return *this;
-}
-
-const Matrix &Matrix::operator*=(const Matrix &rhs){
-  (*this) = (*this)*rhs;
-  return *this;
-}
-
-const Matrix Matrix::operator+(const Matrix &rhs) const{
+Matrix &Matrix::operator=(const Matrix &src)
+{
+  // MatrixContent checks for self-assignment
+  *content = *(src.content);
+  return *this;
+}
+
+const Matrix &Matrix::operator+=(const Matrix &rhs)
+{
+  *content += *(rhs.content);
+  return *this;
+}
+
+const Matrix &Matrix::operator-=(const Matrix &rhs)
+{
+  *content -= *(rhs.content);
+  return *this;
+}
+
+const Matrix &Matrix::operator*=(const Matrix &rhs)
+{
+  (*content) *= (*rhs.content);
+  return *this;
+}
+
+const Matrix Matrix::operator+(const Matrix &rhs) const
+{
   Matrix tmp = *this;
   tmp+=rhs;
@@ -163 +151 @@
 }
 
-const Matrix Matrix::operator-(const Matrix &rhs) const{
+const Matrix Matrix::operator-(const Matrix &rhs) const
+{
   Matrix tmp = *this;
   tmp-=rhs;
@@ -169 +158 @@
 }
 
-const Matrix Matrix::operator*(const Matrix &rhs) const{
-  gsl_matrix *res = gsl_matrix_alloc(NDIM, NDIM);
-  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, content->content, rhs.content->content, 0.0, res);
-  MatrixContent *content= new MatrixContent();
-  content->content = res;
-  return Matrix(content);
-}
-
-double &Matrix::at(size_t i, size_t j){
-  ASSERT(i>=0&&i<NDIM,"Index i for Matrix access out of range");
-  ASSERT(j>=0&&j<NDIM,"Index j for Matrix access out of range");
-  return *gsl_matrix_ptr (content->content, i, j);
-}
-
-const double Matrix::at(size_t i, size_t j) const{
-  ASSERT(i>=0&&i<NDIM,"Index i for Matrix access out of range");
-  ASSERT(j>=0&&j<NDIM,"Index j for Matrix access out of range");
-  return gsl_matrix_get(content->content, i, j);
-}
-
-Vector &Matrix::row(size_t i){
+const Matrix Matrix::operator*(const Matrix &rhs) const
+{
+  Matrix tmp(content);
+  tmp *= rhs;
+  return tmp;
+}
+
+double &Matrix::at(size_t i, size_t j)
+{
+  return content->at(i,j);
+}
+
+const double Matrix::at(size_t i, size_t j) const
+{
+  return content->at(i,j);
+}
+
+Vector &Matrix::row(size_t i)
+{
   ASSERT(i>=0&&i<NDIM,"Index i for Matrix access out of range");
   return *rows_ptr[i];
 }
 
-const Vector &Matrix::row(size_t i) const{
+const Vector &Matrix::row(size_t i) const
+{
   ASSERT(i>=0&&i<NDIM,"Index i for Matrix access out of range");
   return *rows_ptr[i];
 }
 
-Vector &Matrix::column(size_t i){
+Vector &Matrix::column(size_t i)
+{
   ASSERT(i>=0&&i<NDIM,"Index i for Matrix access out of range");
   return *columns_ptr[i];
 }
 
-const Vector &Matrix::column(size_t i) const{
+const Vector &Matrix::column(size_t i) const
+{
   ASSERT(i>=0&&i<NDIM,"Index i for Matrix access out of range");
   return *columns_ptr[i];
 }
 
-Vector &Matrix::diagonal(){
+Vector &Matrix::diagonal()
+{
   return *diagonal_ptr;
 }
 
-const Vector &Matrix::diagonal() const{
+const Vector &Matrix::diagonal() const
+{
   return *diagonal_ptr;
 }
 
-void Matrix::set(size_t i, size_t j, const double value){
-  ASSERT(i>=0&&i<NDIM,"Index i for Matrix access out of range");
-  ASSERT(j>=0&&j<NDIM,"Index j for Matrix access out of range");
-  gsl_matrix_set(content->content,i,j,value);
+void Matrix::set(size_t i, size_t j, const double value)
+{
+  content->set(i,j, value);
 }
 
@@ -253 +244 @@
 }
 
-Matrix Matrix::transpose() const{
-  MatrixContent *newContent = new MatrixContent();
-  newContent->content = gsl_matrix_alloc(NDIM, NDIM);
-  gsl_matrix_transpose_memcpy(newContent->content, content->content);
-  Matrix res = Matrix(newContent);
+Matrix Matrix::transpose() const
+{
+  std::cout << "const Matrix::transpose()." << std::endl;
+  Matrix res = Matrix(const_cast<const MatrixContent *>(content)->transpose());
   return res;
 }
@@ -263 +253 @@
 Matrix &Matrix::transpose()
 {
-  double tmp;
-  for (int i=0;i<NDIM;i++)
-    for (int j=i+1;j<NDIM;j++) {
-      tmp = at(j,i);
-      at(j,i) = at(i,j);
-      at(i,j) = tmp;
-    }
-  return *this;
-}
-
+  std::cout << "Matrix::transpose()." << std::endl;
+  content->transpose();
+  return *this;
+}
 
 Vector Matrix::transformToEigenbasis()
 {
-  gsl_eigen_symmv_workspace *T = gsl_eigen_symmv_alloc(NDIM);
-  gsl_vector *eval = gsl_vector_alloc(NDIM);
-  gsl_matrix *evec = gsl_matrix_alloc(NDIM, NDIM);
-  gsl_eigen_symmv(content->content, eval, evec, T);
-  gsl_eigen_symmv_free(T);
-  gsl_matrix_memcpy(content->content, evec);
-  gsl_matrix_free(evec);
+  gsl_vector *eval = content->transformToEigenbasis();
   Vector evalues(gsl_vector_get(eval,0), gsl_vector_get(eval,1), gsl_vector_get(eval,2));
   gsl_vector_free(eval);
@@ -288 +266 @@
 }
 
-const Matrix &Matrix::operator*=(const double factor){
-  gsl_matrix_scale(content->content, factor);
-  return *this;
-}
-
-const Matrix operator*(const double factor,const Matrix& mat){
+const Matrix &Matrix::operator*=(const double factor)
+    {
+  *content *= factor;
+  return *this;
+}
+
+const Matrix operator*(const double factor,const Matrix& mat)
+{
   Matrix tmp = mat;
-  tmp*=factor;
-  return tmp;
-}
-
-const Matrix operator*(const Matrix &mat,const double factor){
+  return tmp *= factor;
+}
+
+const Matrix operator*(const Matrix &mat,const double factor)
+{
   return factor*mat;
 }
 
-bool Matrix::operator==(const Matrix &rhs) const{
-  for(int i=NDIM;i--;){
-    for(int j=NDIM;j--;){
-      if(fabs(at(i,j)-rhs.at(i,j))>MYEPSILON){
-        return false;
-      }
-    }
-  }
-  return true;
+bool Matrix::operator==(const Matrix &rhs) const
+{
+  return (*content == *(rhs.content));
 }
 
@@ -333 +307 @@
 };
 
-ostream &operator<<(ostream &ost,const Matrix &mat){
+ostream &operator<<(ostream &ost,const Matrix &mat)
+{
   for(int i = 0;i<NDIM;++i){
     ost << "\n";
@@ -345 +320 @@
 }
 
-Vector operator*(const Matrix &mat,const Vector &vec){
-  Vector res;
-  gsl_blas_dgemv( CblasNoTrans, 1.0, mat.content->content, vec.content->content, 0.0, res.content->content);
-  return res;
-}
-
-Vector &operator*=(Vector& lhs,const Matrix &mat){
+Vector operator*(const Matrix &mat,const Vector &vec)
+{
+  return (*mat.content) * vec;
+}
+
+Vector &operator*=(Vector& lhs,const Matrix &mat)
+{
   lhs = mat*lhs;
   return lhs;