Changeset c5038e

Timestamp:

Oct 12, 2011, 2:10:56 PM (13 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:

29cbe9

Parents:

c8b17b

git-author:

Frederik Heber <heber@…> (09/09/11 21:24:02)

git-committer:

Frederik Heber <heber@…> (10/12/11 14:10:56)

Message:

Made MatrixContent and VectorContent dimensions const by placing into extra class.

• new classes MatrixDimension and VectorDimension.
• they only displays getters to the outside. Thus, we may savely serialize them without having to disable the constness.

Location:

LinearAlgebra/src/LinearAlgebra

Files:

• LinearSystemOfEquations.cpp (modified) (3 diffs)
• MatrixContent.cpp (modified) (43 diffs)
• MatrixContent.hpp (modified) (5 diffs)
• MatrixVector_ops.cpp (modified) (2 diffs)
• VectorContent.cpp (modified) (29 diffs)
• VectorContent.hpp (modified) (6 diffs)

LinearAlgebra/src/LinearAlgebra/LinearSystemOfEquations.cpp

@@ -97 +97 @@
 
   // copy solution
-  for (size_t i=0;i<x->dimension;i++) {
+  for (size_t i=0;i<x->getDimension();i++) {
     array[i] = x->at(i);
   }
@@ -113 +113 @@
 
   // copy solution
-  for (size_t i=0;i<x->dimension;i++)
+  for (size_t i=0;i<x->getDimension();i++)
     v[i] = x->at(i);
   return status;
@@ -128 +128 @@
   int s;
   if (IsSymmetric) { // use LU
-    gsl_permutation * p = gsl_permutation_alloc (x->dimension);
+    gsl_permutation * p = gsl_permutation_alloc (x->getDimension());
     gsl_linalg_LU_decomp (A->content, p, &s);
     gsl_linalg_LU_solve (A->content, p, b->content, x->content);

LinearAlgebra/src/LinearAlgebra/MatrixContent.cpp

@@ -33 +33 @@
 #include <set>
 
-using namespace std;
-
 
 /** Constructor for class MatrixContent.
@@ -41 +39 @@
  */
 MatrixContent::MatrixContent(size_t _rows, size_t _columns) :
-    rows(_rows),
-    columns(_columns),
+    MatrixDimension(_rows,_columns),
     free_content_on_exit(true)
 {
-  content = gsl_matrix_calloc(rows, columns);
+  content = gsl_matrix_calloc(getRows(), getColumns());
 }
 
@@ -56 +53 @@
  */
 MatrixContent::MatrixContent(size_t _rows, size_t _columns, MatrixBaseCase) :
-    rows(_rows),
-    columns(_columns),
+    MatrixDimension(_rows,_columns),
     free_content_on_exit(true)
 {}
@@ -67 +63 @@
  */
 MatrixContent::MatrixContent(size_t _rows, size_t _columns, const double *src) :
-    rows(_rows),
-    columns(_columns),
+    MatrixDimension(_rows,_columns),
     free_content_on_exit(true)
 {
-  content = gsl_matrix_calloc(rows, columns);
+  content = gsl_matrix_calloc(getRows(), getColumns());
   set(0,0, src[0]);
   set(1,0, src[1]);
@@ -95 +90 @@
  */
 MatrixContent::MatrixContent(size_t _row, size_t _column, std::istream &inputstream) :
+  MatrixDimension(_row,_column),
   content(NULL),
-  rows(_row),
-  columns(_column),
   free_content_on_exit(true)
 {
    // allocate matrix and set contents
-  content = gsl_matrix_alloc(rows, columns);
+  content = gsl_matrix_alloc(getRows(), getColumns());
 
   size_t row = 0;
@@ -125 +119 @@
 
     // check number of columns parsed
-    ASSERT(line_of_values.size() == columns,
+    ASSERT(line_of_values.size() == getColumns(),
         "MatrixContent::MatrixContent() - row "
         +toString(row)+" has a different number of columns "
         +toString(line_of_values.size())+" than others before "
-        +toString(columns)+".");
-
-    for (size_t column = 0; column < columns; ++column)
+        +toString(getColumns())+".");
+
+    for (size_t column = 0; column < getColumns(); ++column)
       set(row, column, line_of_values[column]);
     ++row;
   } while (inputstream.good());
   // check number of rows parsed
-  ASSERT(row == rows,
+  ASSERT(row == getRows(),
       "MatrixContent::MatrixContent() - insufficent number of rows "
-      +toString(row)+" < "+toString(rows)+" read.");
+      +toString(row)+" < "+toString(getRows())+" read.");
 }
 
@@ -147 +141 @@
  */
 MatrixContent::MatrixContent(gsl_matrix *&src) :
-    rows(src->size1),
-    columns(src->size2),
+    MatrixDimension(src->size1,src->size2),
     free_content_on_exit(true)
 {
@@ -161 +154 @@
  */
 MatrixContent::MatrixContent(const MatrixContent &src) :
-    rows(src.rows),
-    columns(src.columns),
+    MatrixDimension(src.getRows(),src.getColumns()),
     free_content_on_exit(true)
 {
-  content = gsl_matrix_alloc(src.rows, src.columns);
+  content = gsl_matrix_alloc(src.getRows(), src.getColumns());
   gsl_matrix_memcpy(content,src.content);
 }
@@ -173 +165 @@
  */
 MatrixContent::MatrixContent(const MatrixContent *src) :
-    rows(src->rows),
-    columns(src->columns),
+    MatrixDimension(src->getRows(),src->getColumns()),
     free_content_on_exit(true)
 {
   ASSERT(src != NULL, "MatrixContent::MatrixContent - pointer to source matrix is NULL!");
-  content = gsl_matrix_alloc(src->rows, src->columns);
+  content = gsl_matrix_alloc(src->getRows(), src->getColumns());
   gsl_matrix_memcpy(content,src->content);
 }
@@ -190 +181 @@
 }
 
-/** Getter for MatrixContent::rows.
- * \return MatrixContent::rows
- */
-const size_t MatrixContent::getRows() const
-{
-  return rows;
-}
-
-/** Getter for MatrixContent::columns.
- * \return MatrixContent::columns
- */
-const size_t MatrixContent::getColumns() const
-{
-  return columns;
-}
-
 /** Return a VectorViewContent of the \a column -th column vector.
  *
@@ -213 +188 @@
 VectorContent *MatrixContent::getColumnVector(size_t column) const
 {
-  ASSERT(column < columns,
+  ASSERT(column < getColumns(),
       "MatrixContent::getColumnVector() - requested column "+toString(column)
-      +" greater than dimension "+toString(columns));
+      +" greater than dimension "+toString(getColumns()));
   return (new VectorViewContent(gsl_matrix_column(content,column)));
 }
@@ -225 +200 @@
 VectorContent *MatrixContent::getRowVector(size_t row) const
 {
-  ASSERT(row < rows,
+  ASSERT(row < getRows(),
       "MatrixContent::getColumnVector() - requested row "+toString(row)
-      +" greater than dimension "+toString(rows));
+      +" greater than dimension "+toString(getRows()));
   return (new VectorViewContent(gsl_matrix_row(content,row)));
 }
@@ -243 +218 @@
 void MatrixContent::setIdentity()
 {
-  for(int i=rows;i--;){
-    for(int j=columns;j--;){
+  for(int i=getRows();i--;){
+    for(int j=getColumns();j--;){
       set(i,j,(double)(i==j));
     }
@@ -254 +229 @@
 void MatrixContent::setZero()
 {
-  for(int i=rows;i--;){
-    for(int j=columns;j--;){
+  for(int i=getRows();i--;){
+    for(int j=getColumns();j--;){
       set(i,j,0.);
     }
@@ -266 +241 @@
 void MatrixContent::setValue(double _value)
 {
-  for(int i=rows;i--;){
-    for(int j=columns;j--;){
+  for(int i=getRows();i--;){
+    for(int j=getColumns();j--;){
       set(i,j,_value);
     }
@@ -312 +287 @@
 const MatrixContent &MatrixContent::operator*=(const MatrixContent &rhs)
 {
-  ASSERT(rhs.columns == rhs.rows,
-      "MatrixContent::operator*=() - rhs matrix is not square: "+toString(rhs.columns)+" != "+toString(rhs.rows)+".");
-  ASSERT(columns == rhs.rows,
-      "MatrixContent::operator*=() - columns dimension differ: "+toString(columns)+" != "+toString(rhs.rows)+".");
+  ASSERT(rhs.getColumns() == rhs.getRows(),
+      "MatrixContent::operator*=() - rhs matrix is not square: "+toString(rhs.getColumns())+" != "+toString(rhs.getRows())+".");
+  ASSERT(getColumns() == rhs.getRows(),
+      "MatrixContent::operator*=() - columns dimension differ: "+toString(getColumns())+" != "+toString(rhs.getRows())+".");
   (*this) = (*this)*rhs;
   return *this;
@@ -326 +301 @@
 const MatrixContent MatrixContent::operator*(const MatrixContent &rhs) const
 {
-  ASSERT (columns == rhs.rows,
+  ASSERT (getColumns() == rhs.getRows(),
       "MatrixContent::operator*() - dimensions not match for matrix product (a,b)*(b,c) = (a,c):"
-      "("+toString(rows)+","+toString(columns)+")*("+toString(rhs.rows)+","+toString(rhs.columns)+")");
-  gsl_matrix *res = gsl_matrix_alloc(rows, rhs.columns);
+      "("+toString(getRows())+","+toString(getColumns())+")*("+toString(rhs.getRows())+","+toString(rhs.getColumns())+")");
+  gsl_matrix *res = gsl_matrix_alloc(getRows(), rhs.getColumns());
   gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, content, rhs.content, 0.0, res);
   // gsl_matrix is taken over by constructor, hence no free
@@ -344 +319 @@
 const MatrixContent MatrixContent::operator&(const MatrixContent &rhs) const
 {
-  ASSERT ((rows == rhs.rows) && (columns == rhs.columns),
+  ASSERT ((getRows() == rhs.getRows()) && (getColumns() == rhs.getColumns()),
       "MatrixContent::operator&() - dimensions not match for matrix product (a,b) != (b,c):"
-      "("+toString(rows)+","+toString(columns)+") != ("+toString(rhs.rows)+","+toString(rhs.columns)+")");
-  gsl_matrix *res = gsl_matrix_alloc(rows, rhs.columns);
-  for (size_t i=0;i<rows;++i)
-    for (size_t j=0;j<columns;++j)
+      "("+toString(getRows())+","+toString(getColumns())+") != ("+toString(rhs.getRows())+","+toString(rhs.getColumns())+")");
+  gsl_matrix *res = gsl_matrix_alloc(getRows(), rhs.getColumns());
+  for (size_t i=0;i<getRows();++i)
+    for (size_t j=0;j<getColumns();++j)
       gsl_matrix_set(res, i,j, gsl_matrix_get(content, i,j)*gsl_matrix_get(rhs.content, i,j));
   // gsl_matrix is taken over by constructor, hence no free
@@ -367 +342 @@
 const MatrixContent &MatrixContent::operator&=(const MatrixContent &rhs)
 {
-  ASSERT ((rows == rhs.rows) && (columns == rhs.columns),
+  ASSERT ((getRows() == rhs.getRows()) && (getColumns() == rhs.getColumns()),
       "MatrixContent::operator&() - dimensions not match for matrix product (a,b) != (b,c):"
-      "("+toString(rows)+","+toString(columns)+") != ("+toString(rhs.rows)+","+toString(rhs.columns)+")");
-  for (size_t i=0;i<rows;++i)
-    for (size_t j=0;j<columns;++j)
+      "("+toString(getRows())+","+toString(getColumns())+") != ("+toString(rhs.getRows())+","+toString(rhs.getColumns())+")");
+  for (size_t i=0;i<getRows();++i)
+    for (size_t j=0;j<getColumns();++j)
       gsl_matrix_set(content, i,j, gsl_matrix_get(content, i,j)*gsl_matrix_get(rhs.content, i,j));
   return *this;
@@ -385 +360 @@
 double &MatrixContent::at(size_t i, size_t j)
 {
-  ASSERT((i>=0) && (i<rows),
-      "MatrixContent::at() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(rows)+"]");
-  ASSERT((j>=0) && (j<columns),
-      "MatrixContent::at() - Index j="+toString(j)+" for Matrix access out of range [0,"+toString(columns)+"]");
+  ASSERT((i>=0) && (i<getRows()),
+      "MatrixContent::at() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(getRows())+"]");
+  ASSERT((j>=0) && (j<getColumns()),
+      "MatrixContent::at() - Index j="+toString(j)+" for Matrix access out of range [0,"+toString(getColumns())+"]");
   return *gsl_matrix_ptr (content, i, j);
 }
@@ -399 +374 @@
 const double MatrixContent::at(size_t i, size_t j) const
 {
-  ASSERT((i>=0) && (i<rows),
-      "MatrixContent::at() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(rows)+"]");
-  ASSERT((j>=0) && (j<columns),
-      "MatrixContent::at() - Index j="+toString(j)+" for Matrix access out of range [0,"+toString(columns)+"]");
+  ASSERT((i>=0) && (i<getRows()),
+      "MatrixContent::at() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(getRows())+"]");
+  ASSERT((j>=0) && (j<getColumns()),
+      "MatrixContent::at() - Index j="+toString(j)+" for Matrix access out of range [0,"+toString(getColumns())+"]");
   return gsl_matrix_get(content, i, j);
 }
@@ -435 +410 @@
 void MatrixContent::set(size_t i, size_t j, const double value)
 {
-  ASSERT((i>=0) && (i<rows),
-      "MatrixContent::set() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(rows)+"]");
-  ASSERT((j>=0) && (j<columns),
-      "MatrixContent::set() - Index j="+toString(j)+" for Matrix access out of range [0,"+toString(columns)+"]");
+  ASSERT((i>=0) && (i<getRows()),
+      "MatrixContent::set() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(getRows())+"]");
+  ASSERT((j>=0) && (j<getColumns()),
+      "MatrixContent::set() - Index j="+toString(j)+" for Matrix access out of range [0,"+toString(getColumns())+"]");
   gsl_matrix_set(content,i,j,value);
 }
@@ -448 +423 @@
 void MatrixContent::setFromDoubleArray(double * x)
 {
-  gsl_matrix_view m = gsl_matrix_view_array (x, rows, columns);
+  gsl_matrix_view m = gsl_matrix_view_array (x, getRows(), getColumns());
   gsl_matrix_memcpy (content, &m.matrix);
 };
@@ -478 +453 @@
 bool MatrixContent::SwapRowColumn(size_t i, size_t j)
 {
-  ASSERT (rows == columns,
+  ASSERT (getRows() == getColumns(),
       "MatrixContent::SwapRowColumn() - The matrix must be square for swapping row against column to be possible.");
   return (gsl_matrix_swap_rowcol (content, i, j) == GSL_SUCCESS);
@@ -488 +463 @@
 MatrixContent MatrixContent::transposed() const
 {
-  gsl_matrix *res = gsl_matrix_alloc(columns, rows);  // column and row dimensions exchanged!
+  gsl_matrix *res = gsl_matrix_alloc(getColumns(), getRows());  // column and row dimensions exchanged!
   gsl_matrix_transpose_memcpy(res, content);
   MatrixContent newContent(res);
@@ -500 +475 @@
 MatrixContent &MatrixContent::transpose()
 {
-  ASSERT( rows == columns,
-      "MatrixContent::transpose() - cannot transpose onto itself as matrix not square: "+toString(rows)+"!="+toString(columns)+"!");
+  ASSERT( getRows() == getColumns(),
+      "MatrixContent::transpose() - cannot transpose onto itself as matrix not square: "+toString(getRows())+"!="+toString(getColumns())+"!");
   double tmp;
-  for (size_t i=0;i<rows;i++)
-    for (size_t j=i+1;j<rows;j++) {
+  for (size_t i=0;i<getRows();i++)
+    for (size_t j=i+1;j<getRows();j++) {
       tmp = at(j,i);
       at(j,i) = at(i,j);
@@ -519 +494 @@
  *  -# V is a NxN orthogonal, square matrix.
  *
- * The SVD is only done if MatrixContent::rows is larger equal to MatrixContent::columns, i.e.
+ * The SVD is only done if MatrixDimension::rows is larger equal to MatrixDimension::columns, i.e.
  * N >= M.
  *
@@ -529 +504 @@
 void MatrixContent::performSingularValueDecomposition(MatrixContent &V, VectorContent &S)
 {
-  ASSERT(rows >= columns,
+  ASSERT(getRows() >= getColumns(),
       "MatrixContent::performSingularValueDecomposition() - row count "
-      +toString(rows)+" must be larger than or equal to column count "+toString(columns)+".");
+      +toString(getRows())+" must be larger than or equal to column count "+toString(getColumns())+".");
   ASSERT(V.getRows() == V.getColumns(),
       "MatrixContent::performSingularValueDecomposition() - Matrix V must be square!");
-  ASSERT(V.getRows() == columns,
+  ASSERT(V.getRows() == getColumns(),
       "MatrixContent::performSingularValueDecomposition() - Given matrix V does not have the right dimension "
-      +toString(V.getRows())+"!="+toString(rows)+".");
-  ASSERT(S.getDimension() == columns,
+      +toString(V.getRows())+"!="+toString(getRows())+".");
+  ASSERT(S.getDimension() == getColumns(),
       "MatrixContent::performSingularValueDecomposition() - Vector S does not have the right dimension "
-      +toString(S.getDimension())+"!="+toString(columns)+".");
+      +toString(S.getDimension())+"!="+toString(getColumns())+".");
   gsl_matrix *A = content;
-  gsl_vector *work = gsl_vector_alloc(columns);
+  gsl_vector *work = gsl_vector_alloc(getColumns());
   gsl_linalg_SV_decomp(A,V.content,S.content,work);
   gsl_vector_free(work);
@@ -548 +523 @@
 /** Performs a SVD on the contained matrix and solves the linear system \f$Ax=b\f$ with given \a b.
  *
- * The SVD is only done if MatrixContent::rows is larger equal to MatrixContent::columns.
+ * The SVD is only done if MatrixDimension::rows is larger equal to MatrixDimension::columns.
  *
  * \note The matrix content is overwritten in the process.
@@ -559 +534 @@
 VectorContent MatrixContent::solveOverdeterminedLinearEquation(MatrixContent &V, VectorContent &S, const VectorContent &b)
 {
-  ASSERT(b.getDimension() == rows,
+  ASSERT(b.getDimension() == getRows(),
       "MatrixContent::performSingularValueDecomposition() - Given vector b does not have the right dimension "
-      +toString(b.getDimension())+"!="+toString(rows)+".");
+      +toString(b.getDimension())+"!="+toString(getRows())+".");
   ASSERT(V.getRows() == V.getColumns(),
       "MatrixContent::performSingularValueDecomposition() - Matrix V must be square!");
-  ASSERT(V.getRows() == columns,
+  ASSERT(V.getRows() == getColumns(),
       "MatrixContent::performSingularValueDecomposition() - Given matrix V does not have the right dimension "
-      +toString(V.getRows())+"!="+toString(rows)+".");
-  ASSERT(S.getDimension() == columns,
+      +toString(V.getRows())+"!="+toString(getRows())+".");
+  ASSERT(S.getDimension() == getColumns(),
       "MatrixContent::performSingularValueDecomposition() - Vector S does not haver the right dimension "
-      +toString(S.getDimension())+"!="+toString(columns)+".");
-
-  gsl_vector *x = gsl_vector_alloc(columns);
+      +toString(S.getDimension())+"!="+toString(getColumns())+".");
+
+  gsl_vector *x = gsl_vector_alloc(getColumns());
   gsl_linalg_SV_solve(content,V.content,S.content,b.content,x);
   return VectorContent(x, true);
@@ -578 +553 @@
 /** Performs a SVD on the contained matrix and solves the linear system \f$Ax=b\f$ with given \a b.
  *
- * The SVD is only done if MatrixContent::rows is larger equal to MatrixContent::columns.
+ * The SVD is only done if MatrixDimension::rows is larger equal to MatrixDimension::columns.
  *
  * \note The matrix content is overwritten in the process.
@@ -587 +562 @@
 VectorContent MatrixContent::solveOverdeterminedLinearEquation(const VectorContent &b)
 {
-  ASSERT(b.getDimension() == rows,
+  ASSERT(b.getDimension() == getRows(),
       "MatrixContent::performSingularValueDecomposition() - Given vector b does not have the right dimension "
-      +toString(b.getDimension())+"!="+toString(rows)+".");
-  MatrixContent V(columns, columns);
-  VectorContent S(columns);
+      +toString(b.getDimension())+"!="+toString(getRows())+".");
+  MatrixContent V(getColumns(), getColumns());
+  VectorContent S(getColumns());
   performSingularValueDecomposition(V,S);
   return solveOverdeterminedLinearEquation(V,S,b);
@@ -604 +579 @@
 gsl_vector* MatrixContent::transformToEigenbasis()
 {
-  if (rows == columns) { // symmetric
-    gsl_eigen_symmv_workspace *T = gsl_eigen_symmv_alloc(rows);
-    gsl_vector *eval = gsl_vector_alloc(rows);
-    gsl_matrix *evec = gsl_matrix_alloc(rows, rows);
+  if (getRows() == getColumns()) { // symmetric
+    gsl_eigen_symmv_workspace *T = gsl_eigen_symmv_alloc(getRows());
+    gsl_vector *eval = gsl_vector_alloc(getRows());
+    gsl_matrix *evec = gsl_matrix_alloc(getRows(), getRows());
     gsl_eigen_symmv(content, eval, evec, T);
     gsl_eigen_symmv_free(T);
@@ -615 +590 @@
   } else { // non-symmetric
     // blow up gsl_matrix in content to square matrix, fill other components with zero
-    const size_t greaterDimension = rows > columns ? rows : columns;
+    const size_t greaterDimension = getRows() > getColumns() ? getRows() : getColumns();
     gsl_matrix *content_square = gsl_matrix_alloc(greaterDimension, greaterDimension);
     for (size_t i=0; i<greaterDimension; i++) {
       for (size_t j=0; j<greaterDimension; j++) {
-        const double value = ((i < rows) && (j < columns)) ? gsl_matrix_get(content,i,j) : 0.;
+        const double value = ((i < getRows()) && (j < getColumns())) ? gsl_matrix_get(content,i,j) : 0.;
         gsl_matrix_set(content_square, i,j, value);
       }
@@ -629 +604 @@
 
     // solve eigenvalue problem
-    gsl_eigen_nonsymmv_workspace *T = gsl_eigen_nonsymmv_alloc(rows);
+    gsl_eigen_nonsymmv_workspace *T = gsl_eigen_nonsymmv_alloc(getRows());
     gsl_vector_complex *eval = gsl_vector_complex_alloc(greaterDimension);
     gsl_matrix_complex *evec = gsl_matrix_complex_alloc(greaterDimension, greaterDimension);
@@ -656 +631 @@
 
     // copy real-parts of complex eigenvalues and eigenvectors (column-wise orientation)
-    gsl_vector *eval_real = gsl_vector_alloc(columns);
+    gsl_vector *eval_real = gsl_vector_alloc(getColumns());
     size_t I=0;
     for (size_t i=0; i<greaterDimension; i++) { // only copy real space part
@@ -731 +706 @@
 bool MatrixContent::IsPositiveDefinite() const
 {
-  if (rows != columns)  // only possible for square matrices.
+  if (getRows() != getColumns())  // only possible for square matrices.
     return false;
   else
@@ -744 +719 @@
 {
   int signum = 0;
-  ASSERT(rows == columns,
+  ASSERT(getRows() == getColumns(),
       "MatrixContent::Determinant() - determinant can only be calculated for square matrices.");
-  gsl_permutation *p = gsl_permutation_alloc(rows);
+  gsl_permutation *p = gsl_permutation_alloc(getRows());
   gsl_linalg_LU_decomp(content, p, &signum);
   gsl_permutation_free(p);
@@ -815 +790 @@
 void MatrixContent::write(std::ostream &ost) const
 {
-  for (size_t row = 0; row < rows; ++row) {
-    for (size_t column = 0; column < columns; ++column)
+  for (size_t row = 0; row < getRows(); ++row) {
+    for (size_t column = 0; column < getColumns(); ++column)
       ost << at(row, column) << "\t";
     ost << std::endl;
@@ -863 +838 @@
 MatrixContent const operator+(const MatrixContent& a, const MatrixContent& b)
 {
-  ASSERT(a.rows == b.rows, "MatrixContent::operator+() - row counts have to match: "
-      +toString(a.rows)+" != "+toString(b.rows)+"!");
-  ASSERT(a.columns == b.columns, "MatrixContent::operator+() - column counts have to match: "
-      +toString(a.columns)+" != "+toString(b.columns)+"!");
+  ASSERT(a.getRows() == b.getRows(), "MatrixContent::operator+() - row counts have to match: "
+      +toString(a.getRows())+" != "+toString(b.getRows())+"!");
+  ASSERT(a.getColumns() == b.getColumns(), "MatrixContent::operator+() - column counts have to match: "
+      +toString(a.getColumns())+" != "+toString(b.getColumns())+"!");
   MatrixContent x(a);
   x += b;
@@ -879 +854 @@
 MatrixContent const operator-(const MatrixContent& a, const MatrixContent& b)
 {
-  ASSERT(a.rows == b.rows, "MatrixContent::operator+() - row counts have to match: "
-      +toString(a.rows)+" != "+toString(b.rows)+"!");
-  ASSERT(a.columns == b.columns, "MatrixContent::operator+() - column counts have to match: "
-      +toString(a.columns)+" != "+toString(b.columns)+"!");
+  ASSERT(a.getRows() == b.getRows(), "MatrixContent::operator+() - row counts have to match: "
+      +toString(a.getRows())+" != "+toString(b.getRows())+"!");
+  ASSERT(a.getColumns() == b.getColumns(), "MatrixContent::operator+() - column counts have to match: "
+      +toString(a.getColumns())+" != "+toString(b.getColumns())+"!");
   MatrixContent x(a);
   x -= b;
@@ -894 +869 @@
 bool MatrixContent::operator==(const MatrixContent &rhs) const
     {
-  if ((rows == rhs.rows) && (columns == rhs.columns)) {
-    for(int i=rows;i--;){
-      for(int j=columns;j--;){
+  if ((getRows() == rhs.getRows()) && (getColumns() == rhs.getColumns())) {
+    for(int i=getRows();i--;){
+      for(int j=getColumns();j--;){
         if(fabs(at(i,j)-rhs.at(i,j))>LINALG_MYEPSILON()){
           return false;
@@ -911 +886 @@
 {
   ost << "\\begin{pmatrix}";
-  for (size_t i=0;i<mat.rows;i++) {
-    for (size_t j=0;j<mat.columns;j++) {
+  for (size_t i=0;i<mat.getRows();i++) {
+    for (size_t j=0;j<mat.getColumns();j++) {
       ost << mat.at(i,j) << " ";
-      if (j != mat.columns-1)
+      if (j != mat.getColumns()-1)
         ost << "& ";
     }
-    if (i != mat.rows-1)
+    if (i != mat.getRows()-1)
       ost << "\\\\ ";
   }

LinearAlgebra/src/LinearAlgebra/MatrixContent.hpp

@@ -14 +14 @@
 #endif
 
+#include "CodePatterns/Assert.hpp"
 
 /** MatrixContent is a wrapper for gsl_matrix.
@@ -46 +47 @@
 
 class VectorContent;
+namespace boost {
+  namespace serialization {
+    class access;
+  }
+}
 
 /** Dummy structure to create a unique signature.
@@ -52 +58 @@
 struct MatrixBaseCase{};
 
-class MatrixContent
+/** This class safely stores matrix dimensions away.
+ *
+ * This way we simulate const-ness of the dimensions while still allowing
+ * serialization to modify them.
+ *
+ */
+struct MatrixDimension {
+public:
+  MatrixDimension(size_t _rows, size_t _columns) : rows(_rows), columns(_columns) {}
+
+  size_t getRows() const {return rows;}
+  size_t getColumns() const {return columns;}
+
+private:
+  void setRows(size_t _rows) {rows = _rows;}
+  void setColumns(size_t _columns) {columns = _columns;}
+
+  friend class boost::serialization::access;
+  // serialization (due to gsl_vector separate versions needed)
+  template<class Archive>
+  void serialize(Archive & ar, const unsigned int version)
+  {
+    ar & rows;
+    ar & columns;
+  }
+
+private:
+  size_t rows;      // store for internal purposes
+  size_t columns;  // store for internal purposes
+};
+
+class MatrixContent : public MatrixDimension
 {
   friend std::ostream & operator<<(std::ostream &ost, const MatrixContent &mat);
@@ -121 +158 @@
   bool operator==(const MatrixContent &rhs) const;
 
-  const size_t getRows() const;
-  const size_t getColumns() const;
-
   VectorContent *getColumnVector(size_t column) const;
   VectorContent *getRowVector(size_t row) const;
@@ -136 +170 @@
 
   gsl_matrix * content;
-  const size_t rows;      // store for internal purposes
-  const size_t columns;  // store for internal purposes
 
 private:

LinearAlgebra/src/LinearAlgebra/MatrixVector_ops.cpp

@@ -36 +36 @@
 VectorContent const operator*(const VectorContent& vec, const MatrixContent& mat)
 {
-  ASSERT( mat.columns == vec.dimension,
-      "operator*() - column dimension of matrix ("+toString(mat.columns)
-      +") and dimension of this vector ("+toString(vec.dimension)+") not equal!");
-  VectorContent result(mat.rows);
+  ASSERT( mat.getColumns() == vec.getDimension(),
+      "operator*() - column dimension of matrix ("+toString(mat.getColumns())
+      +") and dimension of this vector ("+toString(vec.getDimension())+") not equal!");
+  VectorContent result(mat.getRows());
   gsl_blas_dgemv( CblasNoTrans, 1.0, mat.content, vec.content, 0.0, result.content);
   return result;
@@ -51 +51 @@
 VectorContent const operator*(const MatrixContent& mat, const VectorContent& vec)
 {
-  ASSERT( mat.columns == vec.dimension,
-      "operator*() - column dimension of this matrix ("+toString(mat.columns)
-      +") and dimension of vector ("+toString(vec.dimension)+") not equal!");
-  VectorContent result(mat.rows);
+  ASSERT( mat.getColumns() == vec.getDimension(),
+      "operator*() - column dimension of this matrix ("+toString(mat.getColumns())
+      +") and dimension of vector ("+toString(vec.getDimension())+") not equal!");
+  VectorContent result(mat.getRows());
   gsl_blas_dgemv( CblasNoTrans, 1.0, mat.content, vec.content, 0.0, result.content);
   return result;

LinearAlgebra/src/LinearAlgebra/VectorContent.cpp

@@ -37 +37 @@
  */
 VectorContent::VectorContent(size_t _dim) :
-    dimension(_dim),
+    VectorDimension(_dim),
     free_content_on_exit(true)
 {
-  content = gsl_vector_calloc(dimension);
+  content = gsl_vector_calloc(getDimension());
 }
 
@@ -51 +51 @@
  */
 VectorContent::VectorContent(VectorBaseCase) :
+    VectorDimension(0),
     free_content_on_exit(false)
 {}
@@ -59 +60 @@
  */
 VectorContent::VectorContent(const VectorContent * const src) :
-    dimension(src->dimension),
+    VectorDimension(src->getDimension()),
     free_content_on_exit(true)
 {
-  content = gsl_vector_alloc(dimension);
+  content = gsl_vector_alloc(getDimension());
   gsl_vector_memcpy (content, src->content);
 };
@@ -71 +72 @@
  */
 VectorContent::VectorContent(const VectorContent & src) :
-    dimension(src.dimension),
+    VectorDimension(src.getDimension()),
     free_content_on_exit(true)
 {
-  content = gsl_vector_alloc(dimension);
+  content = gsl_vector_alloc(getDimension());
   gsl_vector_memcpy (content, src.content);
 };
@@ -84 +85 @@
  */
 VectorContent::VectorContent(const size_t _dimension, std::istream &inputstream) :
-  dimension(_dimension),
+  VectorDimension(_dimension),
   content(NULL),
   free_content_on_exit(true)
 {
    // allocate vector and set contents
-  content = gsl_vector_alloc(dimension);
+  content = gsl_vector_alloc(getDimension());
 
   size_t row = 0;
@@ -114 +115 @@
 
     // check number of columns parsed
-    ASSERT(line_of_values.size()+index <= dimension,
+    ASSERT(line_of_values.size()+index <= getDimension(),
         "VectorContent::VectorContent() - exceedingly many entries in row "
         +toString(row)+": current number is "
         +toString(index)+" and expecting now "
-        +toString(line_of_values.size())+" which exceeds "+toString(dimension)+".");
+        +toString(line_of_values.size())+" which exceeds "+toString(getDimension())+".");
 
     for (std::vector<double>::const_iterator iter = line_of_values.begin();
@@ -126 +127 @@
   } while (inputstream.good());
   // check number of rows parsed
-  ASSERT(index == dimension,
+  ASSERT(index == getDimension(),
       "VectorContent::VectorContent() - insufficent number of components "
-      +toString(index)+" < "+toString(dimension)+" read.");
+      +toString(index)+" < "+toString(getDimension())+" read.");
 }
 
@@ -138 +139 @@
  */
 VectorContent::VectorContent(gsl_vector * _src, bool _free_content_on_exit) :
-  dimension(_src->size),
+  VectorDimension(_src->size),
   free_content_on_exit(_free_content_on_exit)
 {
@@ -163 +164 @@
 double &VectorContent::at(size_t i)
 {
-  ASSERT((i>=0) && (i<dimension),
-      "VectorContent::at() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(dimension)+"]");
+  ASSERT((i>=0) && (i<getDimension()),
+      "VectorContent::at() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(getDimension())+"]");
   return *gsl_vector_ptr (content, i);
 }
@@ -174 +175 @@
 const double VectorContent::at(size_t i) const
 {
-  ASSERT((i>=0) && (i<dimension),
-      "VectorContent::at() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(dimension)+"]");
+  ASSERT((i>=0) && (i<getDimension()),
+      "VectorContent::at() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(getDimension())+"]");
   return gsl_vector_get(content, i);
 }
 
 /** These functions return a pointer to the \a m-th element of a vector.
- *  If \a m lies outside the allowed range of 0 to VectorContent::dimension-1 then the error handler is invoked and a null pointer is returned.
+ *  If \a m lies outside the allowed range of 0 to VectorDimension::dimension-1 then the error handler is invoked and a null pointer is returned.
  * \param m m-th element
  * \return pointer to \a m-th element
@@ -190 +191 @@
 
 /** These functions return a constant pointer to the \a m-th element of a vector.
- *  If \a m lies outside the allowed range of 0 to VectorContent::dimension-1 then the error handler is invoked and a null pointer is returned.
+ *  If \a m lies outside the allowed range of 0 to VectorDimension::dimension-1 then the error handler is invoked and a null pointer is returned.
  * \param m m-th element
  * \return const pointer to \a m-th element
@@ -205 +206 @@
 VectorContent& VectorContent::operator=(const VectorContent& src)
 {
-  ASSERT(dimension == src.dimension,
+  ASSERT(getDimension() == src.getDimension(),
       "Dimensions have to be the same to assign VectorContent onto another:"
-      +toString(dimension)+" != "+toString(src.dimension)+"!");
+      +toString(getDimension())+" != "+toString(src.getDimension())+"!");
   // check for self assignment
   if(&src!=this){
@@ -230 +231 @@
 void VectorContent::setFromDoubleArray(double * x)
 {
-  gsl_vector_view m = gsl_vector_view_array (x, dimension);
+  gsl_vector_view m = gsl_vector_view_array (x, getDimension());
   gsl_vector_memcpy (content, &m.vector);
 };
@@ -285 +286 @@
 double &VectorContent::operator[](size_t i)
 {
-  ASSERT((i>=0) && (i<dimension),
-      "VectorContent::operator[]() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(dimension)+"]");
+  ASSERT((i>=0) && (i<getDimension()),
+      "VectorContent::operator[]() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(getDimension())+"]");
   return *gsl_vector_ptr (content, i);
 }
@@ -296 +297 @@
 const double VectorContent::operator[](size_t i) const
 {
-  ASSERT((i>=0) && (i<dimension),
-      "VectorContent::operator[]() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(dimension)+"]");
+  ASSERT((i>=0) && (i<getDimension()),
+      "VectorContent::operator[]() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(getDimension())+"]");
   return gsl_vector_get(content, i);
 }
@@ -310 +311 @@
 {
   bool status = true;
-  ASSERT(dimension == b.dimension, "Dimenions of VectorContents to compare differ");
-  for (size_t i=0;i<dimension;i++)
+  ASSERT(getDimension() == b.getDimension(),
+      "VectorContent::operator==() - Dimensions of VectorContents to compare differ: "
+      +toString(getDimension())+" != "+toString(b.getDimension())+".");
+  for (size_t i=0;i<getDimension();i++)
     status = status && (fabs(at(i) - b.at(i)) <= LINALG_MYEPSILON());
   return status;
@@ -323 +326 @@
 const VectorContent& VectorContent::operator+=(const VectorContent& b)
 {
-  ASSERT(dimension == b.dimension, "Dimenions of VectorContents to compare differ");
-  for (size_t i=0;i<dimension;i++)
+  ASSERT(getDimension() == b.getDimension(),
+      "VectorContent::operator=+() - Dimensions of VectorContents to compare differ: "
+      +toString(getDimension())+" != "+toString(b.getDimension())+".");
+  for (size_t i=0;i<getDimension();i++)
     at(i) = at(i)+b.at(i);
   return *this;
@@ -336 +341 @@
 const VectorContent& VectorContent::operator-=(const VectorContent& b)
 {
-  ASSERT(dimension == b.dimension, "Dimenions of VectorContents to compare differ");
-  for (size_t i=0;i<dimension;i++)
+  ASSERT(getDimension() == b.getDimension(),
+      "VectorContent::operator-=() - Dimensions of VectorContents to compare differ: "
+      +toString(getDimension())+" != "+toString(b.getDimension())+".");
+  for (size_t i=0;i<getDimension();i++)
     at(i) = at(i)-b.at(i);
   return *this;
@@ -349 +356 @@
 const VectorContent& VectorContent::operator*=(const double m)
 {
-  for (size_t i=0;i<dimension;i++)
+  for (size_t i=0;i<getDimension();i++)
     at(i) = at(i)*m;
   return *this;
@@ -361 +368 @@
 VectorContent const operator+(const VectorContent& a, const VectorContent& b)
 {
-  ASSERT(a.dimension == b.dimension, "VectorContent::operator+() - dimensions have to match: "+toString(a.dimension)+" != "+toString(b.dimension)+"!");
+  ASSERT(a.getDimension() == b.getDimension(),
+      "VectorContent::operator+() - dimensions have to match: "
+      +toString(a.getDimension())+" != "+toString(b.getDimension())+"!");
   VectorContent x(a);
-  for (size_t i=0;i<a.dimension;i++)
+  for (size_t i=0;i<a.getDimension();i++)
     x.at(i) = a.at(i)+b.at(i);
   return x;
@@ -375 +384 @@
 VectorContent const operator-(const VectorContent& a, const VectorContent& b)
 {
-  ASSERT(a.dimension == b.dimension, "VectorContent::operator-() - dimensions have to match: "+toString(a.dimension)+" != "+toString(b.dimension)+"!");
+  ASSERT(a.getDimension() == b.getDimension(),
+      "VectorContent::operator-() - dimensions have to match: "
+      +toString(a.getDimension())+" != "+toString(b.getDimension())+"!");
   VectorContent x(a);
-  for (size_t i=0;i<a.dimension;i++)
+  for (size_t i=0;i<a.getDimension();i++)
     x.at(i) = a.at(i)-b.at(i);
   return x;
@@ -390 +401 @@
 {
   VectorContent x(a);
-  for (size_t i=0;i<a.dimension;i++)
+  for (size_t i=0;i<a.getDimension();i++)
     x.at(i) = a.at(i)*m;
   return x;
@@ -403 +414 @@
 {
   VectorContent x(a);
-  for (size_t i=0;i<a.dimension;i++)
+  for (size_t i=0;i<a.getDimension();i++)
     x.at(i) = a.at(i)*m;
   return x;
@@ -411 +422 @@
 {
   ost << "[";
-  for (size_t i=0;i<m.dimension;i++) {
+  for (size_t i=0;i<m.getDimension();i++) {
     ost << m.at(i);
-    if (i != m.dimension-1)
+    if (i != m.getDimension()-1)
       ost << " ";
   }
@@ -429 +440 @@
 {
   double result = 0.;
-  for (size_t i = dimension; i--; )
+  for (size_t i = getDimension(); i--; )
     result += fabs(at(i));
   return (result <= LINALG_MYEPSILON());
@@ -440 +451 @@
 {
   double NormValue = 0.;
-  for (size_t i=dimension;--i;)
+  for (size_t i=getDimension();--i;)
     NormValue += at(i)*at(i);
   return (fabs(NormValue - 1.) <= LINALG_MYEPSILON());
@@ -484 +495 @@
 {
   double res = 0.;
-  for (int i=dimension;i--;)
+  for (int i=getDimension();i--;)
     res += (at(i)-y[i])*(at(i)-y[i]);
   return (res);
@@ -570 +581 @@
 void VectorContent::write(std::ostream &ost) const
 {
-  for (size_t index = 0; index < dimension; ++index)
+  for (size_t index = 0; index < getDimension(); ++index)
     ost << at(index) << "\t";
   ost << std::endl;
@@ -587 +598 @@
   return res;
 };
-
-/** Getter for VectorContent::dimension.
- * @return VectorContent::dimension
- */
-const size_t VectorContent::getDimension() const
-{
-  return dimension;
-}

LinearAlgebra/src/LinearAlgebra/VectorContent.hpp

@@ -14 +14 @@
 #endif
 
+#include "CodePatterns/Assert.hpp"
 
 /**
@@ -32 +33 @@
 class Vector;
 class MatrixContent;
+struct VectorViewContent;
+namespace boost {
+  namespace serialization {
+    class access;
+  }
+}
+
+/** This class safely stores vector dimension away.
+ *
+ * This way we simulate const-ness of the dimension while still allowing
+ * serialization to modify them.
+ *
+ */
+struct VectorDimension {
+  friend struct VectorViewContent;
+public:
+  VectorDimension(size_t _dimension) : dimension(_dimension) {}
+
+  size_t getDimension() const {return dimension;}
+
+private:
+  void setDim(size_t _dimension) {dimension = _dimension;}
+
+  friend class boost::serialization::access;
+  // serialization (due to gsl_vector separate versions needed)
+  template<class Archive>
+  void serialize(Archive & ar, const unsigned int version)
+  {
+    ar & dimension;
+  }
+
+private:
+  size_t dimension;      // store for internal purposes
+};
 
 /** Dummy structure to create a unique signature.
@@ -38 +73 @@
 struct VectorBaseCase{};
 
-class VectorContent{
+class VectorContent : public VectorDimension
+{
   friend std::ostream & operator<< (std::ostream& ost, const VectorContent &m);
   friend VectorContent const operator*(const VectorContent& a, const double m);
@@ -58 +94 @@
   virtual ~VectorContent();
 
+  static size_t preparseVectorDimensions(std::istream &inputstream);
+
   // Accessing
   double &at(size_t m);
@@ -65 +103 @@
   double *Pointer(size_t m) const;
   const double *const_Pointer(size_t m) const;
+  void write(std::ostream &ost) const;
 
   // Assignment operator
@@ -106 +145 @@
   const double operator*(const VectorContent& b) const;
 
-  const size_t getDimension() const;
-
-  static size_t preparseVectorDimensions(std::istream &inputstream);
-  void write(std::ostream &ost) const;
-
-  size_t dimension;
   gsl_vector *content;
 private: