dealii_solver_rol_vector.hpp

#ifndef __DEALII_SOLVER_ROL_VECTOR_H__
#define __DEALII_SOLVER_ROL_VECTOR_H__

#include "ROL_Types.hpp"
#include "ROL_Vector_SimOpt.hpp"


#include "optimization/rol_to_dealii_vector.hpp"


template<typename Real = double>
class dealiiSolverVectorWrappingROL
{
private:
    ROL::Ptr<ROL::Vector<Real>> rol_vector_ptr;


    void print(const ROL::Vector<Real> &rol_vector) const
    {
        const ROL::Vector_SimOpt<Real> *vec_split12 = dynamic_cast<const ROL::Vector_SimOpt<Real> *>(&rol_vector);
        if (vec_split12 == NULL) {
            PHiLiP::ROL_vector_to_dealii_vector_reference(rol_vector).print(std::cout,14);
        } else {
            const auto vec_1 = vec_split12->get_1();
            print(*vec_1);
            const auto vec_2 = vec_split12->get_2();
            print(*vec_2);
        }
    }
public:
    using value_type = Real;


    dealiiSolverVectorWrappingROL() = default;

    explicit dealiiSolverVectorWrappingROL(ROL::Ptr<ROL::Vector<Real>> input_vector)
    : rol_vector_ptr(input_vector)
    {};

    ROL::Ptr<ROL::Vector<Real>> getVector()
    {
        return rol_vector_ptr;
    }

    ROL::Ptr<const ROL::Vector<Real>> getVector() const
    {
        return rol_vector_ptr;
    }
    
    void reinit (const dealiiSolverVectorWrappingROL &model_vector,
                 const bool leave_elements_uninitialized = false)
    {
        rol_vector_ptr = model_vector.getVector()->clone();
        if (!leave_elements_uninitialized) {
            (*this) *= 0.0;
        }
    }
    Real operator * (const dealiiSolverVectorWrappingROL &v) const
    {
        return rol_vector_ptr->dot( *(v.getVector()) );
    }

    dealiiSolverVectorWrappingROL & operator=(const Real a)
    {
        rol_vector_ptr->setScalar(a);
        return *this;
    }

    dealiiSolverVectorWrappingROL & operator=(const dealiiSolverVectorWrappingROL &x)
    {
        rol_vector_ptr = (x.getVector())->clone();
        return *this;
    }
    
    dealiiSolverVectorWrappingROL & operator*=(const Real a)
    {
        rol_vector_ptr->scale(a);
        return *this;
    }

    void add (const dealiiSolverVectorWrappingROL &x)
    {
        this->add(1.0,x);
    }

    void add (const Real  a,
              const dealiiSolverVectorWrappingROL &x)
    {
        rol_vector_ptr->axpy(a, *(x.getVector()) );
    }

    void sadd (const Real  a,
               const Real  b,
               const dealiiSolverVectorWrappingROL &x)
    {
        (*this) *= a;
        this->add(b, x);
    }

    void equ (const Real  a,
              const dealiiSolverVectorWrappingROL &x)
    {
        rol_vector_ptr->set( *(x.getVector()) );
        (*this) *= a;
    }

    Real add_and_dot (const Real  a,
                        const dealiiSolverVectorWrappingROL &x,
                        const dealiiSolverVectorWrappingROL &v)
    {
        this->add(a, x);
        return (*this) * v;
    }

    Real l2_norm () const
    {
        return std::sqrt( (*this) * (*this) );
    }

    Teuchos::RCP<dealiiSolverVectorWrappingROL> basis(int i) const
    {
        ROL::Ptr<ROL::Vector<Real>> rol_basis = rol_vector_ptr->basis(i);

        Teuchos::RCP<dealiiSolverVectorWrappingROL> e = Teuchos::rcp(new dealiiSolverVectorWrappingROL(rol_basis));

        return e;
    }

    int size() const
    {
        return rol_vector_ptr->dimension();
    }

    void print() const
    {
        print(*rol_vector_ptr);
        //     const auto &vec_split12 = dynamic_cast<const ROL::Vector_SimOpt<Real>&>(*rol_vector_ptr);
        //     const auto des = vec_split12.get_1();
        //     const auto &des_split = dynamic_cast<const ROL::Vector_SimOpt<Real>&>(*des);

        //     const auto sim_des = des_split.get_1();
        //     const auto ctl_des = des_split.get_2();
        //     const auto con = vec_split12.get_2();

        //     PHiLiP::ROL_vector_to_dealii_vector_reference(*sim_des).print(std::cout,14);
        //     PHiLiP::ROL_vector_to_dealii_vector_reference(*ctl_des).print(std::cout,14);
        //     PHiLiP::ROL_vector_to_dealii_vector_reference(*con).print(std::cout,14);
    }


    Real operator [](int i) const
    {
        const auto &vec_split12 = dynamic_cast<const ROL::Vector_SimOpt<Real>&>(*rol_vector_ptr);
        const auto des = vec_split12.get_1();
        const auto &des_split = dynamic_cast<const ROL::Vector_SimOpt<Real>&>(*des);

        const auto sim_des = des_split.get_1();
        const auto ctl_des = des_split.get_2();
        const auto con = vec_split12.get_2();

        const auto n1 = sim_des->dimension();
        const auto n2 = ctl_des->dimension() + n1;
        //const auto n3 = con->dimension() + n2;

        if (i < n1) {
            return PHiLiP::ROL_vector_to_dealii_vector_reference(*sim_des)[i];
        } else if (i < n2) {
            return PHiLiP::ROL_vector_to_dealii_vector_reference(*ctl_des)[i-n1];
        } else {
            return PHiLiP::ROL_vector_to_dealii_vector_reference(*con)[i-n2];
        }
    }

};

#endif