SparseMatrixBuilder.h

#pragma once
 
#include "eigen_matrix_defs.h"
#include "ObjectFactoryStack.h"
#include <vector>
#include <map>
#include "LinearOperator.h"
#include "AsynchronousTask.h"
#include "MatrixConnectome.h"
#include <Eigen/StdVector>
 
namespace CoRiBoDynamics
{
class SparseMatrixBuilder: public LinearOperator
{
public:
    explicit SparseMatrixBuilder(CoreBoundThreadPool* ThreadPool);
    ~SparseMatrixBuilder();
 
    void StartBuilding(int size); 
    void InitializeElementMatrix(int i, const mat6& M);
    void AddToMatrixThreadSafe( int i, int j, const mat6& M ); 
 
    virtual double ApplyLinearOperator( const Vector& x, Vector& y ); 
    virtual double ApplyApproximateInverse( const Vector& b, Vector& x ); 
    virtual void ComputeApproximateInverse(); 
 
    virtual int SystemSize();
 
    mat6 GetElementMatrix(int index);
 
    void SetDropSize(int size){m_connectome.set_drop_size(size);}
 
protected:
    void AddToMatrix(int i, int j, const mat6& M); 
 
    /*
        Build Eigen-format subset-matrix by picking all 6x6 elements from the given iterator
    */
    void FillSubMatrix(SparseMat& result, std::vector<int>::iterator elements_begin, std::vector<int>::iterator elements_end);
 
 
    /*
        Extra-diagonal elements are split into 'm_connection_indices' and 'm_connection_matrices'
        since m_connection_indices is kept explicitly sorted during creation, we avoid copy-move of 6x6 matrices by splitting
    */
    std::vector<MatrixConnectome::CRix> m_connection_indices; // list of all extra-diagonal connections in the system, with and index into m_connection_matrices to find the actual numerical values of the connection. sorted in terms of collumn index first and row index second
    std::vector<mat6,Eigen::aligned_allocator<mat6>> m_connection_matrices; // list of all extra diagonal connection matrices
    std::mutex m_ConnectionElementLock; // all potential concurrent access to m_connection_indices and m_connection_matrices must be protected
 
    int m_size; 
    std::mutex* m_MatrixColumnLocks;
 
    MatrixConnectome m_connectome;
 
    mat6* m_A_0_diagonal; 
    mat6* m_A_1_diagonal; 
 
    mat6* m_U_0_diagonal; 
    mat6* m_U_1_diagonal; 
 
    bool* m_set; // which elements in m_A_1_diagonal are non zero?
 
    class ConnectionSolverTask : public CoreBoundThreadPool::Task {
    public:
        void Execute();
        std::vector<int> elements;
        Eigen::SimplicialLLT<SparseMat, Eigen::Upper, Eigen::NaturalOrdering<int>>* m_sparse_solver;
        //Eigen::LLT<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>, Eigen::Upper> m_dense_solver;
        SparseMat connection_matrix;
        const Vector* global_B;
        Vector b;
        Vector x;
        ConnectionSolverTask() { m_sparse_solver = nullptr; }
        ~ConnectionSolverTask() { delete m_sparse_solver; }
 
        double accumulator;
        enum State {COMPUTE,SOLVE} state;
        SparseMatrixBuilder* m_builder;
    };
 
    std::vector<ConnectionSolverTask> m_connection_solver;
 
    class SubMultiplicationTask: public CoreBoundThreadPool::Task{
    public:
        explicit SubMultiplicationTask();
        void setSubMatrix(int start_index, int end_index, mat6* A0, mat6* A1);
        void MultiplySubDomain(const Vector* x, Vector* y);
 
        double GetAccumulator(){return m_accumulator;}
    protected:
        void Execute();
 
        int m_ix0; 
        int m_ix1; 
        mat6* m_A0; 
        mat6* m_A1; 
        const Vector* m_x;
        Vector* m_y;
        double m_accumulator; // sub-matrix contribution to the conjugate dot product x^t * A * x
    };
 
 
    class TridiagonalSubSolverTask: public CoreBoundThreadPool::Task {
    public:
        explicit TridiagonalSubSolverTask();
        void SetMemoryRoot(mat6* U0_root, mat6* U1_root);
        ~TridiagonalSubSolverTask();
        void InvertSubDomain(int start_index, int end_index, mat6* A0, mat6* A1);
        void SolveSubDomain(const Vector* b, Vector* x);
 
        double GetAccumulator(){return m_accumulator;}
    protected:
        enum State {COMPUTE,SOLVE};
        void Compute(); // 6x6 block tri-diagonal matrix decomposition
        void Solve(); // back-substitution solution of the linear system for the current b vector
 
        void Execute();
 
        int m_ix0; 
        int m_ix1; 
        int m_num;
        mat6* m_A0; 
        mat6* m_A1; 
        const Vector* m_b;
        Vector* m_x;
        mat6* m_U1;
        mat6* m_U0;
        State m_state;
 
        mat6* m_U1_root;
        mat6* m_U0_root;
 
        double m_accumulator; // sub-matrix contribution to the conjugate dot product x^t * A^-1 * x
    };
 
    std::vector<TridiagonalSubSolverTask> m_TridiagonalSubSolverTask; 
    std::vector<SubMultiplicationTask> m_SubMultiplicationTask; 
 
    CoreBoundThreadPool* m_ThreadPool; 
 
};
}