NNLS_solver.h

/* Note: This solver is based on a modified version of the algorithm descrived in "SOLVING 
 * LEAST SQUARES PROBLEMS". This is because in the context of Energy Conserving Sampling
 * and Weighting a modified exit condition was introduced. The details regarding the exit
 * condition (EQUATION 13) are described in the following paper:
 * https://onlinelibrary.wiley.com/doi/full/10.1002/nme.5332
 * More details on the ECSW hyper-reduction technique can be found in:
 * https://onlinelibrary.wiley.com/doi/full/10.1002/nme.4820
 * 
 * Additionally Functionality added to the Eigen NNLS:
 * -Option to use Gradient Exit Condition (used in Eigen/original textbook) or Residual Exit Condition (default - Introduced in ECSW work) (discussed further in NNLS_solver.cpp)
 * -Option to use an iterative linear solver or direct solve (default) for LS problem in the algorithm
 * -Option to input a transposed matrix
*/

#ifndef NNLS_H
#define NNLS_H

#include <iostream>
#include <string>
#include <AztecOO_config.h>
#include <Epetra_MpiComm.h>
#include <Epetra_ConfigDefs.h>
#include <Epetra_Map.h>
#include <Epetra_CrsMatrix.h>
#include <Epetra_Vector.h>
#include <Epetra_MultiVector.h>
#include <Epetra_Import.h>
#include <Epetra_LinearProblem.h>
#include <EpetraExt_MatrixMatrix.h>
#include <AztecOO.h>
#include <Amesos.h>
#include <Amesos_BaseSolver.h>
#include <eigen/Eigen/Dense>
#include "parameters/all_parameters.h"
#include "reduced_order/multi_core_helper_functions.h"

using Eigen::Matrix;
namespace PHiLiP {
class NNLSSolver
{
public:
    NNLSSolver(
        const Parameters::AllParameters *const parameters_input,
        const dealii::ParameterHandler &parameter_handler_input,
        const Epetra_CrsMatrix &A, 
        Epetra_MpiComm &Comm, 
        Epetra_Vector &b);

    NNLSSolver(
        const Parameters::AllParameters *const parameters_input,
        const dealii::ParameterHandler &parameter_handler_input,
        const Epetra_CrsMatrix &A,
        const bool is_input_A_matrix_transposed, 
        Epetra_MpiComm &Comm, 
        Epetra_Vector &b);

    NNLSSolver(
        const Parameters::AllParameters *const parameters_input,
        const dealii::ParameterHandler &parameter_handler_input,
        const Epetra_CrsMatrix &A, 
        Epetra_MpiComm &Comm, 
        Epetra_Vector &b, 
        bool grad_exit_crit);

    NNLSSolver(
        const Parameters::AllParameters *const parameters_input,
        const dealii::ParameterHandler &parameter_handler_input,
        const Epetra_CrsMatrix &A,
        const bool is_input_A_matrix_transposed, 
        Epetra_MpiComm &Comm, 
        Epetra_Vector &b, 
        bool grad_exit_crit);
    
    NNLSSolver(        
        const Parameters::AllParameters *const parameters_input,
        const dealii::ParameterHandler &parameter_handler_input,
        const Epetra_CrsMatrix &A, 
        Epetra_MpiComm &Comm, 
        Epetra_Vector &b, 
        bool iter_solver, 
        int LS_iter, 
        double LS_tol);

    NNLSSolver(        
        const Parameters::AllParameters *const parameters_input,
        const dealii::ParameterHandler &parameter_handler_input,
        const Epetra_CrsMatrix &A,
        const bool is_input_A_matrix_transposed, 
        Epetra_MpiComm &Comm, 
        Epetra_Vector &b, 
        bool iter_solver, 
        int LS_iter, 
        double LS_tol);

    NNLSSolver(
        const Parameters::AllParameters *const parameters_input,
        const dealii::ParameterHandler &parameter_handler_input,
        const Epetra_CrsMatrix &A, 
        const bool is_input_A_matrix_transposed,
        Epetra_MpiComm &Comm, 
        Epetra_Vector &b, 
        bool grad_exit_crit, 
        bool iter_solver, 
        int LS_iter, 
        double LS_tol);
 
    ~NNLSSolver() {};

    bool solve();

    Epetra_Vector & get_solution() {return multi_x_;}

    void starting_solution(Epetra_Vector &start) {
        Epetra_Vector start_single_core = allocate_vector_to_single_core(start);
        P.flip();
        sub_into_x(start_single_core);
        P.flip();}
  
    const Parameters::AllParameters *const all_parameters; 

    const dealii::ParameterHandler &parameter_handler;

protected:
    Epetra_MpiComm Comm_;
    const Epetra_CrsMatrix A_;
    Epetra_Vector b_;
    Epetra_Vector x_;
    Epetra_Vector multi_x_;
    int LS_iter_;
    double LS_tol_; 
    int numInactive_;
    std::vector<bool> Z;
    std::vector<bool> P;

public:
    bool iter_solver_;
    bool grad_exit_crit_;
    const bool is_input_A_matrix_transposed_;
    int iter_;

protected:
    Eigen::VectorXd index_set;

private:
    void epetra_permutation_matrix(Epetra_CrsMatrix &P_mat);

    void positive_set_matrix(Epetra_CrsMatrix &P_mat);

    void sub_into_x(Epetra_Vector &temp);

    void add_into_x(Epetra_Vector &temp, double alpha);

    void move_to_active_set(int idx);

    void move_to_inactive_set(int idx);

};
} // PHiLiP namespace
#endif  // NNLS_H