Method_EMA: WIP implementation for sparse matrices

core/include/engine/Eigenmodes.hpp

@@ -33,7 +33,14 @@ bool Hessian_Full_Spectrum(
 bool Hessian_Partial_Spectrum(
     const std::shared_ptr<Data::Parameters_Method> parameters, const vectorfield & spins, const vectorfield & gradient,
     const MatrixX & hessian, int n_modes, MatrixX & tangent_basis, MatrixX & hessian_constrained, VectorX & eigenvalues,
-    MatrixX & eigenvectors );
+    MatrixX & eigenvcetors );
+
+bool Sparse_Hessian_Partial_Spectrum(
+    const std::shared_ptr<Data::Parameters_Method> parameters, const vectorfield & spins, const vectorfield & gradient,
+    const SpMatrixX & hessian, int n_modes, SpMatrixX & tangent_basis, SpMatrixX & hessian_constrained, VectorX & eigenvalues,
+    MatrixX & eigenvectors
+);
+
 }; // end namespace Eigenmodes
 } // end namespace Engine
 #endif

core/include/engine/Manifoldmath.hpp

@@ -87,6 +87,9 @@ void spherical_to_cartesian_coordinate_basis( const vectorfield & vf, MatrixX &
 //
 void spherical_coordinate_christoffel_symbols( const vectorfield & vf, MatrixX & gamma_theta, MatrixX & gamma_phi );
 
+void sparse_hessian_bordered_3N(
+    const vectorfield & image, const vectorfield & gradient, const SpMatrixX & hessian, SpMatrixX & hessian_out );
+
 // Calculate Hessian for a vectorfield constrained to unit length, at any extremum (i.e. where vectors || gradient)
 void hessian_bordered(
     const vectorfield & image, const vectorfield & gradient, const MatrixX & hessian, MatrixX & tangent_basis,

core/src/engine/Eigenmodes.cpp

@@ -4,6 +4,8 @@
 // #include <engine/Backend_par.hpp>
 
 #include <SymEigsSolver.h> // Also includes <MatOp/DenseSymMatProd.h>
+#include <MatOp/SparseSymMatProd.h> // Also includes <MatOp/DenseSymMatProd.h>
+
 #include <Eigen/Dense>
 #include <Eigen/Eigenvalues>
 
@@ -65,7 +67,6 @@ void Calculate_Eigenmodes( std::shared_ptr<Data::Spin_System> system, int idx_im
 
     // Calculate the Eigenmodes
     vectorfield gradient( nos );
-    MatrixX hessian( 3 * nos, 3 * nos );
 
     // The gradient (unprojected)
     system->hamiltonian->Gradient( spins_initial, gradient );
@@ -75,18 +76,41 @@ void Calculate_Eigenmodes( std::shared_ptr<Data::Spin_System> system, int idx_im
     //     g[idx] = mask[idx]*g[idx];
     // });
     Vectormath::set_c_a( 1, gradient, gradient, system->geometry->mask_unpinned );
+
 
-    // The Hessian (unprojected)
-    system->hamiltonian->Hessian( spins_initial, hessian );
-
-    // Get the eigenspectrum
-    MatrixX hessian_constrained = MatrixX::Zero( 2 * nos, 2 * nos );
-    MatrixX tangent_basis       = MatrixX::Zero( 3 * nos, 2 * nos );
     VectorX eigenvalues;
     MatrixX eigenvectors;
-    bool successful = Eigenmodes::Hessian_Partial_Spectrum(
-        system->ema_parameters, spins_initial, gradient, hessian, n_modes, tangent_basis, hessian_constrained,
-        eigenvalues, eigenvectors );
+    SpMatrixX tangent_basis = SpMatrixX(3*nos, 2*nos);
+
+    // TODO: do this properly
+    bool sparse = true;
+    bool successful;
+    if(sparse)
+    {
+        // The Hessian (unprojected)
+        SpMatrixX hessian( 3 * nos, 3 * nos );
+        system->hamiltonian->Sparse_Hessian( spins_initial, hessian );
+        // Get the eigenspectrum
+        SpMatrixX hessian_constrained = SpMatrixX( 2 * nos, 2 * nos );
+
+        successful = Eigenmodes::Sparse_Hessian_Partial_Spectrum(
+            system->ema_parameters, spins_initial, gradient, hessian, n_modes, tangent_basis, hessian_constrained,
+            eigenvalues, eigenvectors );
+
+    } else {
+        // The Hessian (unprojected)
+        MatrixX hessian( 3 * nos, 3 * nos );
+        system->hamiltonian->Hessian( spins_initial, hessian );
+        // Get the eigenspectrum
+        MatrixX hessian_constrained = MatrixX::Zero( 2 * nos, 2 * nos );
+        MatrixX _tangent_basis      = MatrixX(tangent_basis);
+
+        successful = Eigenmodes::Hessian_Partial_Spectrum(
+            system->ema_parameters, spins_initial, gradient, hessian, n_modes, _tangent_basis, hessian_constrained,
+            eigenvalues, eigenvectors );
+
+        tangent_basis = _tangent_basis.sparseView();
+    }
 
     if( successful )
     {
@@ -216,5 +240,46 @@ bool Hessian_Partial_Spectrum(
     return ( hessian_spectrum.info() == Spectra::SUCCESSFUL ) && ( nconv > 0 );
 }
 
+bool Sparse_Hessian_Partial_Spectrum(
+    const std::shared_ptr<Data::Parameters_Method> parameters, const vectorfield & spins, const vectorfield & gradient,
+    const SpMatrixX & hessian, int n_modes, SpMatrixX & tangent_basis, SpMatrixX & hessian_constrained, VectorX & eigenvalues,
+    MatrixX & eigenvectors )
+{
+    int nos = spins.size();
+
+    // Restrict number of calculated modes to [1,2N)
+    n_modes = std::max( 1, std::min( 2 * nos - 2, n_modes ) );
+
+    // Calculate the final Hessian to use for the minimum mode
+    Manifoldmath::sparse_tangent_basis_spherical(spins, tangent_basis);
+
+    SpMatrixX hessian_constrained_3N = SpMatrixX(3*nos, 3*nos);
+    Manifoldmath::sparse_hessian_bordered_3N( spins, gradient, hessian, hessian_constrained_3N );
+
+    hessian_constrained = tangent_basis.transpose() * hessian_constrained_3N * tangent_basis;
+
+    // TODO: Pinning (see non-sparse function for)
+
+    // Create the Spectra Matrix product operation
+    Spectra::SparseSymMatProd<scalar> op( hessian_constrained );
+    // Create and initialize a Spectra solver
+    Spectra::SymEigsSolver<scalar, Spectra::SMALLEST_ALGE, Spectra::SparseSymMatProd<scalar>> hessian_spectrum(
+        &op, n_modes, 2 * nos );
+    hessian_spectrum.init();
+
+    // Compute the specified spectrum, sorted by smallest real eigenvalue
+    int nconv = hessian_spectrum.compute( 1000, 1e-10, int( Spectra::SMALLEST_ALGE ) );
+
+    // Extract real eigenvalues
+    eigenvalues = hessian_spectrum.eigenvalues().real();
+
+    // Retrieve the real eigenvectors
+    eigenvectors = hessian_spectrum.eigenvectors().real();
+
+    // Return whether the calculation was successful
+    return ( hessian_spectrum.info() == Spectra::SUCCESSFUL ) && ( nconv > 0 );
+}
+
+
 } // namespace Eigenmodes
 } // namespace Engine

core/src/engine/Manifoldmath.cpp

@@ -469,6 +469,36 @@ void spherical_to_cartesian_christoffel_symbols( const vectorfield & vf, MatrixX
     }
 }
 
+void sparse_hessian_bordered_3N(
+    const vectorfield & image, const vectorfield & gradient, const SpMatrixX & hessian, SpMatrixX & hessian_out )
+{
+    // Calculates a 3Nx3N matrix in the bordered Hessian approach and transforms it into the tangent basis,
+    // making the result a 2Nx2N matrix. The bordered Hessian's Lagrange multipliers assume a local extremum.
+
+    int nos = image.size();
+    VectorX lambda( nos );
+    for( int i = 0; i < nos; ++i )
+        lambda[i] = image[i].normalized().dot( gradient[i] );
+
+    // Construct hessian_out
+    typedef Eigen::Triplet<scalar> T;
+    std::vector<T> tripletList;
+    tripletList.reserve( hessian.nonZeros() + 3 * nos );
+
+    // Iterate over non zero entries of hesiian
+    for( int k = 0; k < hessian.outerSize(); ++k )
+    {
+        for( SpMatrixX::InnerIterator it( hessian, k ); it; ++it )
+        {
+            tripletList.push_back( T( it.row(), it.col(), it.value() ) );
+        }
+        int j = k % 3;
+        int i = ( k - j ) / 3;
+        tripletList.push_back( T( k, k, -lambda[i] ) ); // Correction to the diagonal
+    }
+    hessian_out.setFromTriplets( tripletList.begin(), tripletList.end() );
+}
+
 void hessian_bordered(
     const vectorfield & image, const vectorfield & gradient, const MatrixX & hessian, MatrixX & tangent_basis,
     MatrixX & hessian_out )
@@ -499,6 +529,7 @@ void hessian_bordered(
     hessian_out = tangent_basis.transpose() * tmp_3N * tangent_basis;
 }
 
+
 void hessian_projected(
     const vectorfield & image, const vectorfield & gradient, const MatrixX & hessian, MatrixX & tangent_basis,
     MatrixX & hessian_out )