ixaxpy.hpp

Go to the documentation of this file.

#ifndef IXAXPY_HPP_
#define IXAXPY_HPP_

#include "jmath/jblas.hpp"
#include "jmath/indirectArray.hpp"

namespace jafar {
  namespace jmath {
    namespace ublasExtra {

      using namespace std;
      using namespace jblas;
      using namespace boost::numeric::ublas;


      // SPARSE SYMMETRIC MATRIX PRODUCT
      template<class bubMatrix>
      void ixaxpy_offdiag_subprod(sym_mat & S, const ind_array & ia_invariant, const bubMatrix & J_oi,
          const ind_array & ia_in, const ind_array & ia_out)
      {
        project(S, ia_out, ia_invariant) = prod(J_oi, project(S, ia_in, ia_invariant));
      }

      template<class bubMatrix, class bubMatrixY>
      void ixaxpy_diag_subprod(sym_mat & S, const bubMatrix & J, const ind_array & ia_in, const ind_array & ia_out,
          const bubMatrixY Y)
      {
        matrix<double> JSii = prod(J, project(S, ia_in, ia_in));
        ublas::noalias(project(S, ia_out, ia_out)) = prod<sym_mat> (JSii, trans(J));
        project(S, ia_out, ia_out) += Y;
      }

      template<class bubMatrix>
      void ixaxpy_diag_subprod(sym_mat & S, const bubMatrix & J, const ind_array & ia_in, const ind_array & ia_out)
      {
        matrix<double> JSii = prod(J, project(S, ia_in, ia_in));
        ublas::noalias(project(S, ia_out, ia_out)) = prod<sym_mat> (JSii, trans(J));
      }


      template<class bubMatrixJ, class bubMatrixY>
      void ixaxpy_prod(sym_mat & X, const ind_array & ia_invariant, const bubMatrixJ & A_oi, const ind_array & ia_in,
          const ind_array & ia_out, const bubMatrixY Y_oo)
      {
        ixaxpy_offdiag_subprod(X, ia_invariant, A_oi, ia_in, ia_out);
        ixaxpy_diag_subprod(X, A_oi, ia_in, ia_out, Y_oo);
      }

      template<class bubMatrixJ>
      void ixaxpy_prod(sym_mat & X, const ind_array & ia_invariant, const bubMatrixJ & A_oi, const ind_array & ia_in,
          const ind_array & ia_out)
      {
        ixaxpy_offdiag_subprod(X, ia_invariant, A_oi, ia_in, ia_out);
        ixaxpy_diag_subprod(X, A_oi, ia_in, ia_out);
      }

    }
  }
}

#endif /* IXAXPY_HPP_ */