Changeset 742371

Timestamp:

Dec 4, 2010, 11:56:28 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:

f5bca2

Parents:

a062e1

git-author:

Frederik Heber <heber@…> (11/18/10 18:24:28)

git-committer:

Frederik Heber <heber@…> (12/04/10 23:56:28)

Message:

Subspace Factorizer is almost working.

• for the testing (4x4) matrix, the middle eigenvectors are good, but the smallest and biggest are off by a few percent, yet converged ...

Location:

src

Files:

• LinearAlgebra/MatrixContent.cpp (modified) (1 diff)
• unittests/SubspaceFactorizerUnittest.cpp (modified) (5 diffs)
• unittests/SubspaceFactorizerUnittest.hpp (modified) (1 diff)

src/LinearAlgebra/MatrixContent.cpp

@@ -253 +253 @@
       "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.columns)+".");
+      "MatrixContent::operator*=() - columns dimension differ: "+toString(columns)+" != "+toString(rhs.rows)+".");
   (*this) = (*this)*rhs;
   return *this;

src/unittests/SubspaceFactorizerUnittest.cpp

@@ -25 +25 @@
 
 #include <gsl/gsl_vector.h>
+#include <boost/foreach.hpp>
 #include <boost/shared_ptr.hpp>
 
@@ -126 +127 @@
 /** For given set of row and column indices, we extract the small block matrix.
  * @param bigmatrix big matrix to extract from
- * @param rowindexset row index set
- * @param columnindexset column index set
+ * @param Eigenvectors eigenvectors of the subspaces to obtain matrix in
+ * @param columnindexset index set to pick out of all indices
  * @return small matrix with dimension equal to the number of indices for row and column.
  */
 MatrixContent * getSubspaceMatrix(
     MatrixContent &bigmatrix,
-    const IndexSet &rowindexset,
-    const IndexSet &columnindexset)
+    VectorArray &Eigenvectors,
+    const IndexSet &indexset)
 {
   // check whether subsystem is big enough for both index sets
-  ASSERT(rowindexset.size() <= bigmatrix.getRows(),
+  ASSERT(indexset.size() <= bigmatrix.getRows(),
       "embedSubspaceMatrix() - bigmatrix has less rows "+toString(bigmatrix.getRows())
       +" than needed by index set "
-      +toString(rowindexset.size())+"!");
-  ASSERT(columnindexset.size() <= bigmatrix.getColumns(),
+      +toString(indexset.size())+"!");
+  ASSERT(indexset.size() <= bigmatrix.getColumns(),
       "embedSubspaceMatrix() - bigmatrix has less columns "+toString(bigmatrix.getColumns())
       +" than needed by index set "
-      +toString(columnindexset.size())+"!");
-  MatrixContent *subsystem =  new MatrixContent(rowindexset.size(), columnindexset.size());
+      +toString(indexset.size())+"!");
+
+  // construct small matrix
+  MatrixContent *subsystem =  new MatrixContent(indexset.size(), indexset.size());
   size_t localrow = 0; // local row indices for the subsystem
   size_t localcolumn = 0;
-  for (IndexSet::const_iterator rowindex = rowindexset.begin();
-      rowindex != rowindexset.end();
-      ++rowindex) {
+  BOOST_FOREACH( size_t rowindex, indexset) {
     localcolumn = 0;
-    for (IndexSet::const_iterator columnindex = columnindexset.begin();
-        columnindex != columnindexset.end();
-        ++columnindex) {
-      ASSERT((*rowindex < bigmatrix.getRows()) && (*columnindex < bigmatrix.getColumns()),
+    BOOST_FOREACH( size_t columnindex, indexset) {
+      ASSERT((rowindex < bigmatrix.getRows()) && (columnindex < bigmatrix.getColumns()),
           "current index pair ("
-          +toString(*rowindex)+","+toString(*columnindex)
+          +toString(rowindex)+","+toString(columnindex)
           +") is out of bounds of bigmatrix ("
           +toString(bigmatrix.getRows())+","+toString(bigmatrix.getColumns())
           +")");
-      subsystem->at(localrow,localcolumn) = bigmatrix.at(*rowindex, *columnindex);
+      subsystem->at(localrow,localcolumn) = (*Eigenvectors[rowindex]) * (bigmatrix * (*Eigenvectors[columnindex]));
       localcolumn++;
     }
@@ -170 +169 @@
 /** For a given set of row and columns indices, we embed a small block matrix into a bigger space.
  *
- * @param bigmatrix big matrix to place submatrix into
- * @param rowindexset
- * @param columnindexset
- * @return
- */
-void embedSubspaceMatrix(
-    MatrixContent &bigmatrix,
+ * @param eigenvectors current eigenvectors
+ * @param rowindexset row index set
+ * @param columnindexset column index set
+ * @return bigmatrix with eigenvectors contained
+ */
+MatrixContent * embedSubspaceMatrix(
+    VectorArray &Eigenvectors,
     MatrixContent &subsystem,
-    const IndexSet &rowindexset,
     const IndexSet &columnindexset)
 {
   // check whether bigmatrix is at least as big as subsystem
-  ASSERT(subsystem.getRows() <= bigmatrix.getRows(),
-      "embedSubspaceMatrix() - subsystem has more rows "+toString(subsystem.getRows())
-      +" than bigmatrix "
-      +toString(bigmatrix.getRows())+"!");
-  ASSERT(subsystem.getColumns() <= bigmatrix.getColumns(),
-      "embedSubspaceMatrix() - subsystem has more columns "+toString(subsystem.getColumns())
-      +" than bigmatrix "
-      +toString(bigmatrix.getColumns())+"!");
+  ASSERT(Eigenvectors.size() > 0,
+      "embedSubspaceMatrix() - no Eigenvectors given!");
+  ASSERT(subsystem.getRows() <= Eigenvectors[0]->getDimension(),
+      "embedSubspaceMatrix() - subsystem has more rows "
+      +toString(subsystem.getRows())+" than eigenvectors "
+      +toString(Eigenvectors[0]->getDimension())+"!");
+  ASSERT(subsystem.getColumns() <= Eigenvectors.size(),
+      "embedSubspaceMatrix() - subsystem has more columns "
+      +toString(subsystem.getColumns())+" than eigenvectors "
+      +toString(Eigenvectors.size())+"!");
   // check whether subsystem is big enough for both index sets
-  ASSERT(rowindexset.size() <= subsystem.getRows(),
-      "embedSubspaceMatrix() - subsystem has less rows "+toString(subsystem.getRows())
-      +" than needed by index set "
-      +toString(rowindexset.size())+"!");
-  ASSERT(columnindexset.size() <= subsystem.getColumns(),
-      "embedSubspaceMatrix() - subsystem has less columns "+toString(subsystem.getColumns())
-      +" than needed by index set "
+  ASSERT(subsystem.getColumns() == subsystem.getRows(),
+      "embedSubspaceMatrix() - subsystem is not square "
+      +toString(subsystem.getRows())+" than needed by index set "
+      +toString(subsystem.getColumns())+"!");
+  ASSERT(columnindexset.size() == subsystem.getColumns(),
+      "embedSubspaceMatrix() - subsystem has not the same number of columns "
+      +toString(subsystem.getColumns())+" compared to the index set "
       +toString(columnindexset.size())+"!");
-  size_t localrow = 0; // local row indices for the subsystem
+
+  // construct intermediate matrix
+  MatrixContent *intermediatematrix = new MatrixContent(Eigenvectors[0]->getDimension(), columnindexset.size());
   size_t localcolumn = 0;
-  for (IndexSet::const_iterator rowindex = rowindexset.begin();
-      rowindex != rowindexset.end();
-      ++rowindex) {
-    localcolumn = 0;
-    for (IndexSet::const_iterator columnindex = columnindexset.begin();
-        columnindex != columnindexset.end();
-        ++columnindex) {
-      bigmatrix.at(*rowindex, *columnindex) = subsystem.at(localrow,localcolumn);
-      localcolumn++;
-    }
-    localrow++;
-  }
-}
-
-/** Operator for output to std:: ostream operator of an IndexSet.
+  BOOST_FOREACH(size_t columnindex, columnindexset) {
+    // create two vectors from each row and copy assign them
+    boost::shared_ptr<VectorContent> srceigenvector(Eigenvectors[columnindex]);
+    boost::shared_ptr<VectorContent> desteigenvector(intermediatematrix->getColumnVector(localcolumn));
+    *desteigenvector = *srceigenvector;
+    localcolumn++;
+  }
+  // matrix product with eigenbasis subsystem matrix
+  *intermediatematrix *= subsystem;
+
+  // and place at right columns into bigmatrix
+  MatrixContent *bigmatrix = new MatrixContent(Eigenvectors[0]->getDimension(), Eigenvectors.size());
+  bigmatrix->setZero();
+  localcolumn = 0;
+  BOOST_FOREACH(size_t columnindex, columnindexset) {
+    // create two vectors from each row and copy assign them
+    boost::shared_ptr<VectorContent> srceigenvector(intermediatematrix->getColumnVector(localcolumn));
+    boost::shared_ptr<VectorContent> desteigenvector(bigmatrix->getColumnVector(columnindex));
+    *desteigenvector = *srceigenvector;
+    localcolumn++;
+  }
+
+  return bigmatrix;
+}
+
+/** Prints the scalar product of each possible pair that is not orthonormal.
+ * We use class logger for printing.
+ * @param AllIndices set of all possible indices of the eigenvectors
+ * @param CurrentEigenvectors array of eigenvectors
+ * @return true - all are orthonormal to each other,
+ *        false - some are not orthogonal or not normalized.
+ */
+bool checkOrthogonality(const IndexSet &AllIndices, const VectorArray &CurrentEigenvectors)
+{
+  size_t nonnormalized = 0;
+  size_t nonorthogonal = 0;
+  // check orthogonality
+  BOOST_FOREACH( size_t firstindex, AllIndices) {
+    BOOST_FOREACH( size_t secondindex, AllIndices) {
+      const double scp = (*CurrentEigenvectors[firstindex])*(*CurrentEigenvectors[secondindex]);
+      if (firstindex == secondindex) {
+        if (fabs(scp - 1.) > MYEPSILON) {
+          nonnormalized++;
+          Log() << Verbose(1) << "Vector " << firstindex << " is not normalized, off by "
+              << fabs(1.-(*CurrentEigenvectors[firstindex])*(*CurrentEigenvectors[secondindex])) << std::endl;
+        }
+      } else {
+        if (fabs(scp) > MYEPSILON) {
+          nonorthogonal++;
+          Log() << Verbose(1) << "Scalar product between " << firstindex << " and " << secondindex
+              << " is " << (*CurrentEigenvectors[firstindex])*(*CurrentEigenvectors[secondindex]) << std::endl;
+        }
+      }
+    }
+  }
+
+  if ((nonnormalized == 0) && (nonorthogonal == 0)) {
+    Log() << Verbose(1) << "All vectors are orthonormal to each other." << std::endl;
+    return true;
+  }
+  if ((nonnormalized == 0) && (nonorthogonal != 0))
+    Log() << Verbose(1) << "All vectors are normalized." << std::endl;
+  if ((nonnormalized != 0) && (nonorthogonal == 0))
+    Log() << Verbose(1) << "All vectors are orthogonal to each other." << std::endl;
+  return false;
+}
+
+/** Calculate the sum of the scalar product of each possible pair.
+ * @param AllIndices set of all possible indices of the eigenvectors
+ * @param CurrentEigenvectors array of eigenvectors
+ * @return sum of scalar products between all possible pairs
+ */
+double calculateOrthogonalityThreshold(const IndexSet &AllIndices, const VectorArray &CurrentEigenvectors)
+{
+  double threshold = 0.;
+  // check orthogonality
+  BOOST_FOREACH( size_t firstindex, AllIndices) {
+    BOOST_FOREACH( size_t secondindex, AllIndices) {
+      const double scp = (*CurrentEigenvectors[firstindex])*(*CurrentEigenvectors[secondindex]);
+      if (firstindex == secondindex) {
+        threshold += fabs(scp - 1.);
+      } else {
+        threshold += fabs(scp);
+      }
+    }
+  }
+  return threshold;
+}
+
+/** Operator for output to std::ostream operator of an IndexSet.
  * @param ost output stream
  * @param indexset index set to output
@@ -231 +308 @@
 }
 
+/** Assign eigenvectors of subspace to full eigenvectors.
+ * We use parallelity as relation measure.
+ * @param eigenvalue eigenvalue to assign along with
+ * @param CurrentEigenvector eigenvector to assign, is taken over within
+ *        boost::shared_ptr
+ * @param CurrentEigenvectors full eigenvectors
+ * @param CorrespondenceList list to make sure that each subspace eigenvector
+ *        is assigned to a unique full eigenvector
+ * @param ParallelEigenvectorList list of "similar" subspace eigenvectors per
+ *        full eigenvector, allocated
+ */
+void AssignSubspaceEigenvectors(
+    double eigenvalue,
+    VectorContent *CurrentEigenvector,
+    VectorArray &CurrentEigenvectors,
+    IndexSet &CorrespondenceList,
+    VectorValueList *&ParallelEigenvectorList)
+{
+  Log() << Verbose(1) << "Current Eigenvector is " << *CurrentEigenvector << std::endl;
+
+  // (for now settle with the one we are most parallel to)
+  size_t mostparallel_index = 4;
+  double mostparallel_scalarproduct = 0.;
+  BOOST_FOREACH( size_t indexiter, CorrespondenceList) {
+    Log() << Verbose(2) << "Comparing to old " << indexiter << "th eigenvector " << *(CurrentEigenvectors[indexiter]) << std::endl;
+    const double scalarproduct = (*(CurrentEigenvectors[indexiter])) * (*CurrentEigenvector);
+    Log() << Verbose(2) << "SKP is " << scalarproduct << std::endl;
+    if (fabs(scalarproduct) > mostparallel_scalarproduct) {
+      mostparallel_scalarproduct = fabs(scalarproduct);
+      mostparallel_index = indexiter;
+    }
+  }
+  if (mostparallel_index != 4) {
+    // put into std::list for later use
+    // invert if pointing in negative direction
+    if ((*(CurrentEigenvectors[mostparallel_index])) * (*CurrentEigenvector) < 0) {
+      *CurrentEigenvector *= -1.;
+      Log() << Verbose(1) << "Pushing (inverted) " << *CurrentEigenvector << " into parallel list [" << mostparallel_index << "]" << std::endl;
+    } else {
+      Log() << Verbose(1) << "Pushing " << *CurrentEigenvector << " into parallel list [" << mostparallel_index << "]" << std::endl;
+    }
+    ParallelEigenvectorList[mostparallel_index].push_back(make_pair(boost::shared_ptr<VectorContent>(CurrentEigenvector), eigenvalue));
+    CorrespondenceList.erase(mostparallel_index);
+  }
+}
+
 void SubspaceFactorizerUnittest::EigenvectorTest()
 {
   VectorArray CurrentEigenvectors;
-  VectorList *ParallelEigenvectorList = new VectorList[4];
+  ValueArray CurrentEigenvalues;
+  VectorValueList *ParallelEigenvectorList = new VectorValueList[4];
   IndexSet AllIndices;
 
@@ -258 +382 @@
 
   // set to first guess, i.e. the unit vectors of R^4
-  for (IndexSet::const_iterator index = AllIndices.begin();
-      index != AllIndices.end();
-      ++index) {
+  BOOST_FOREACH( size_t index, AllIndices) {
     boost::shared_ptr<VectorContent> EV(new VectorContent(4));
     EV->setZero();
-    EV->at(*index) = 1.;
+    EV->at(index) = 1.;
     CurrentEigenvectors.push_back(EV);
-  }
-
-  // for every dimension
-  for (size_t dim = 0; dim<3;++dim) {
-    // for every index set of this dimension
-    Log() << Verbose(0) << std::endl << std::endl;
-    Log() << Verbose(0) << "Current dimension is " << dim << std::endl;
-    std::pair<IndexMap::const_iterator,IndexMap::const_iterator> Bounds = Dimension_to_Indexset.equal_range(dim);
-    for (IndexMap::const_iterator IndexsetIter = Bounds.first;
-        IndexsetIter != Bounds.second;
-        ++IndexsetIter) {
-      // show the index set
-      Log() << Verbose(0) << std::endl;
-      Log() << Verbose(1) << "Current index set is " << *(IndexsetIter->second) << std::endl;
-
-      // create transformation matrices from these
-      MatrixContent *subsystem = getSubspaceMatrix(*fourbyfour, *(IndexsetIter->second), *(IndexsetIter->second));
-      Log() << Verbose(2) << "Subsystem matrix is " << *subsystem << std::endl;
-
-      // solve _small_ systems for eigenvalues
-      VectorContent *Eigenvalues = new VectorContent(subsystem->transformToEigenbasis());
-      Log() << Verbose(2) << "Eigenvector matrix is " << *subsystem << std::endl;
-      Log() << Verbose(2) << "Eigenvalues are " << *Eigenvalues << std::endl;
-      delete Eigenvalues;
-
-      // blow up eigenvectors to 4dim column vector again
-      MatrixContent *Eigenvectors = new MatrixContent(4,4);
-      Eigenvectors->setZero();
-      embedSubspaceMatrix(*Eigenvectors, *subsystem, *(IndexsetIter->second), *(IndexsetIter->second));
-      Log() << Verbose(1) << "4x4 Eigenvector matrix is " << *Eigenvectors << std::endl;
-
-      // we don't need the subsystem anymore
-      delete subsystem;
-
-      // copy the index set, we look for one new eigenvector for each old one
-      IndexSet CorrespondenceList((*IndexsetIter->second));
-
-      // go through all eigenvectors in this subspace
-      for (IndexSet::const_iterator iter = (*IndexsetIter->second).begin();
-          iter != (*IndexsetIter->second).end();
-          ++iter) {
-        // we rather copy the column vector, as the containing matrix is destroyed lateron
-        VectorContent *CurrentEigenvector = new VectorContent(Eigenvectors->getColumnVector(*iter));
-        Log() << Verbose(1) << "Current Eigenvector is " << *CurrentEigenvector << std::endl;
-
-        // recognize eigenvectors parallel to existing ones
-        // (for now settle with the one we are most parallel to)
-        size_t mostparallel_index = 4;
-        double mostparallel_scalarproduct = 0.;
-        for (IndexSet::const_iterator indexiter = CorrespondenceList.begin();
-            indexiter != CorrespondenceList.end();
-            ++indexiter) {
-          Log() << Verbose(2) << "Comparing to old " << *indexiter << "th eigenvector " << *(CurrentEigenvectors[*indexiter]) << std::endl;
-          const double scalarproduct = (*(CurrentEigenvectors[*indexiter])) * (*CurrentEigenvector);
-          Log() << Verbose(2) << "SKP is " << scalarproduct << std::endl;
-          if (fabs(scalarproduct) > mostparallel_scalarproduct) {
-            mostparallel_scalarproduct = fabs(scalarproduct);
-            mostparallel_index = *indexiter;
+    CurrentEigenvalues.push_back(0.);
+  }
+
+  size_t run=1; // counting iterations
+  double threshold = 1.;  // containing threshold value
+  while ((threshold > 1e-6) && (run < 200)) {
+    // for every dimension
+    for (size_t dim = 0; dim<3;++dim) {
+      // for every index set of this dimension
+      Log() << Verbose(0) << std::endl << std::endl;
+      Log() << Verbose(0) << "Current dimension is " << dim << std::endl;
+      std::pair<IndexMap::const_iterator,IndexMap::const_iterator> Bounds = Dimension_to_Indexset.equal_range(dim);
+      for (IndexMap::const_iterator IndexsetIter = Bounds.first;
+          IndexsetIter != Bounds.second;
+          ++IndexsetIter) {
+        // show the index set
+        Log() << Verbose(0) << std::endl;
+        Log() << Verbose(1) << "Current index set is " << *(IndexsetIter->second) << std::endl;
+
+        // create transformation matrices from these
+        MatrixContent *subsystem = getSubspaceMatrix(*fourbyfour, CurrentEigenvectors, *(IndexsetIter->second));
+        Log() << Verbose(2) << "Subsystem matrix is " << *subsystem << std::endl;
+
+        // solve _small_ systems for eigenvalues
+        VectorContent *Eigenvalues = new VectorContent(subsystem->transformToEigenbasis());
+        Log() << Verbose(2) << "Eigenvector matrix is " << *subsystem << std::endl;
+        Log() << Verbose(2) << "Eigenvalues are " << *Eigenvalues << std::endl;
+
+        // blow up eigenvectors to 4dim column vector again
+        MatrixContent *Eigenvectors = embedSubspaceMatrix(CurrentEigenvectors, *subsystem, *(IndexsetIter->second));
+        Log() << Verbose(1) << "4x4 Eigenvector matrix is " << *Eigenvectors << std::endl;
+
+        // we don't need the subsystem anymore
+        delete subsystem;
+
+        // go through all eigenvectors in this subspace
+        IndexSet CorrespondenceList((*IndexsetIter->second)); // assure one-to-one and onto assignment
+        size_t localindex = 0;
+        BOOST_FOREACH( size_t iter, (*IndexsetIter->second)) {
+          // recognize eigenvectors parallel to existing ones
+          AssignSubspaceEigenvectors(
+              Eigenvalues->at(localindex),
+              new VectorContent(Eigenvectors->getColumnVector(iter)),
+              CurrentEigenvectors,
+              CorrespondenceList,
+              ParallelEigenvectorList);
+          localindex++;
+        }
+
+        // free eigenvectors
+        delete Eigenvectors;
+        delete Eigenvalues;
+      }
+    }
+
+    // print list of similar eigenvectors
+    BOOST_FOREACH( size_t index, AllIndices) {
+      Log() << Verbose(2) << "Similar to " << index << "th current eigenvector " << *(CurrentEigenvectors[index]) << " are:" << std::endl;
+      BOOST_FOREACH( VectorValueList::value_type &iter, ParallelEigenvectorList[index] ) {
+        Log() << Verbose(2) << *(iter.first) << std::endl;
+      }
+      Log() << Verbose(2) << std::endl;
+    }
+
+    // create new CurrentEigenvectors from averaging parallel ones.
+    BOOST_FOREACH(size_t index, AllIndices) {
+      CurrentEigenvectors[index]->setZero();
+      CurrentEigenvalues[index] = 0.;
+      BOOST_FOREACH( VectorValueList::value_type &iter, ParallelEigenvectorList[index] ) {
+        *CurrentEigenvectors[index] += (*iter.first) * (iter.second);
+        CurrentEigenvalues[index] += (iter.second);
+      }
+      *CurrentEigenvectors[index] *= 1./CurrentEigenvalues[index];
+      CurrentEigenvalues[index] /= (double)ParallelEigenvectorList[index].size();
+      ParallelEigenvectorList[index].clear();
+    }
+
+    // check orthonormality
+    threshold = calculateOrthogonalityThreshold(AllIndices, CurrentEigenvectors);
+    bool dontOrthonormalization = checkOrthogonality(AllIndices, CurrentEigenvectors);
+
+    // orthonormalize
+    if (!dontOrthonormalization) {
+      Log() << Verbose(0) << "Orthonormalizing ... " << std::endl;
+      for (IndexSet::const_iterator firstindex = AllIndices.begin();
+          firstindex != AllIndices.end();
+          ++firstindex) {
+        for (IndexSet::const_iterator secondindex = firstindex;
+            secondindex != AllIndices.end();
+            ++secondindex) {
+          if (*firstindex == *secondindex) {
+            (*CurrentEigenvectors[*secondindex]) *= 1./(*CurrentEigenvectors[*secondindex]).Norm();
+          } else {
+            (*CurrentEigenvectors[*secondindex]) -=
+                ((*CurrentEigenvectors[*firstindex])*(*CurrentEigenvectors[*secondindex]))
+                *(*CurrentEigenvectors[*firstindex]);
           }
         }
-        if (mostparallel_index != 4) {
-          // put into std::list for later use
-          Log() << Verbose(1) << "Pushing " << *CurrentEigenvector << " into parallel list [" << mostparallel_index << "]" << std::endl;
-          ParallelEigenvectorList[mostparallel_index].push_back(boost::shared_ptr<VectorContent>(CurrentEigenvector));
-          CorrespondenceList.erase(mostparallel_index);
-        }
-        // no need to delete CurrentEigenvector as is taken care of by shared_ptr
-      }
-      // free eigenvectors
-      delete Eigenvectors;
-    }
-  }
-
-  // print list of parallel eigenvectors
-  for (IndexSet::const_iterator index = AllIndices.begin();
-      index != AllIndices.end();
-      ++index) {
-    Log() << Verbose(0) << "Similar to " << *index << "th current eigenvector " << *(CurrentEigenvectors[*index]) << " are:" << std::endl;
-    for (VectorList::const_iterator iter = ParallelEigenvectorList[*index].begin();
-        iter != ParallelEigenvectorList[*index].end();
-        ++iter) {
-      Log() << Verbose(0) << **iter << std::endl;
-    }
-    Log() << Verbose(0) << std::endl;
-  }
-
-  // create new CurrentEigenvectors from averaging parallel ones.
-  for (IndexSet::const_iterator index = AllIndices.begin();
-      index != AllIndices.end();
-      ++index) {
-    for (VectorList::const_iterator iter = ParallelEigenvectorList[*index].begin();
-        iter != ParallelEigenvectorList[*index].end();
-        ++iter) {
-      *CurrentEigenvectors[*index] += **iter;
-    }
-    *CurrentEigenvectors[*index] *= 1./(double)(ParallelEigenvectorList[*index].size()+1);
-  }
-
-  // show new ones
-  Log() << Verbose(0) << "Resulting new eigenvectors are:" << std::endl;
-  for (IndexSet::const_iterator index = AllIndices.begin();
-      index != AllIndices.end();
-      ++index) {
-    Log() << Verbose(0) << *CurrentEigenvectors[*index] << std::endl;
+      }
+    }
+
+    // check orthonormality again
+    checkOrthogonality(AllIndices, CurrentEigenvectors);
+
+    // show new ones
+    Log() << Verbose(0) << "Resulting new eigenvectors and -values, run " << run << " are:" << std::endl;
+    BOOST_FOREACH( size_t index, AllIndices) {
+      Log() << Verbose(0) << *CurrentEigenvectors[index] << " with " << CurrentEigenvalues[index] << std::endl;
+    }
+    run++;
   }
 

src/unittests/SubspaceFactorizerUnittest.hpp

@@ -14 +14 @@
 typedef std::multimap< size_t, boost::shared_ptr<IndexSet> > IndexMap;
 typedef std::list< boost::shared_ptr<VectorContent> > VectorList;
+typedef std::list< std::pair<boost::shared_ptr<VectorContent>, double> > VectorValueList;
 typedef std::vector< boost::shared_ptr<VectorContent> > VectorArray;
+typedef std::vector< double > ValueArray;
 
 class MatrixContent;