ConjugateGradient.h

#pragma once
#include "LinearOperator.h"
#include "AsynchronousTask.h"
 
namespace CoRiBoDynamics{
    class ConjugateGradientSolver{
    public:
        explicit ConjugateGradientSolver(LinearOperator* A);
        void SetThreadPool(CoreBoundThreadPool* ThreadPool);
 
        void solve(const Vector& B, Vector& X, double& tolerance, int& iterations);
 
        void ParallellizeVectorOperations(bool b){m_parallelize_vector_operations = b;}
 
    protected:
        /* |a|^2 */
        double CG_squared_norm(const Vector& a);
 
        /* b = a-b */
        void CG_flip_and_add(const Vector& a, Vector& b);
 
        /*  b = a + beta*b */
        void CG_scale_and_add(const Vector& a, double beta, Vector& b);
        
        /*
            x = x + alpha*p
            r = r - alpha tmp
            return |r|^2
        */
        double MultiOP(Vector& p, Vector& tmp, double alpha, Vector& x, Vector& r);
 
        LinearOperator* m_A;
        Vector m_residual;
        Vector m_p;
        Vector m_z;
        Vector m_tmp;
 
        CoreBoundThreadPool* m_ThreadPool;
        bool m_parallelize_vector_operations;
 
 
        class VectorTask : public CoreBoundThreadPool::Task {
        public:
            VectorTask();
            void SetIndex(int ix0, int ix1);
 
            virtual void Execute();
 
            void CG_squared_norm(const Vector& a);
            void CG_flip_and_add(const Vector& a, Vector& b);
            void CG_scale_and_add(const Vector& a, double beta, Vector& b);
 
            void MultiOP(Vector& p, Vector& tmp, double alpha, Vector& x, Vector& r);
            
            double GetAccumulator();
 
        protected:
            enum OPERATION {NOP, SQUARE_NORM, FLIP_AND_ADD, SCALE_AND_ADD, MULTI_OP} m_operation;
 
            int m_ix0;
            int m_N;
 
            const Vector* m_ca;
            const Vector* m_cb;
            Vector* m_a;
            Vector* m_b;
            double m_alpha;
            double m_beta;
            double m_accumulator;
        };
        std::vector<VectorTask> m_vector_task;
    };
}