MrhsReshaper.h

Go to the documentation of this file.

// ========================================================================= //
// This file is part of MyraMath, copyright (c) 2014-2019 by Ryan A Chilton  //
// and distributed by MyraCore, LLC. See LICENSE.txt for license terms.      //
// ========================================================================= //

#ifndef MYRAMATH_ITERATIVE_MRHSRESHAPER_H
#define MYRAMATH_ITERATIVE_MRHSRESHAPER_H

#include <myramath/dense/VectorRange.h>
#include <myramath/dense/MatrixRange.h>
#include <myramath/dense/Matrix.h>

// The multiple-rhs (Matrix B/X) iterative solvers just delegate to the single-rhs (Vector b/x) ones,
// by reshaping/concatenating the columns of the Matrix's B/X into (taller) Vector's b/x. In the lucky
// event that B and X are already column-contiguious (LD = I), all this reshaping can be done in-place.
// Otherwise, must make deep copies of B and X that are contiguous, delegate, then copy the solution back.
//
// These two classes encapsulate those steps/decisions for B->b and X->x, to reduce code duplication/errors.

namespace myra {
namespace detail {

// Captures a Matrix-valued forcing data (B), presents it as a Vector (b), either by reshaping or deep copying.
template<class Number> class MrhsReshaperB
  {
  public:

    // Constructor, requires forcing data B.
    MrhsReshaperB(const CMatrixRange<Number>& in_B)
      : B(in_B)
      {
      // Is B column-contiguous? No copy is needed..
      if (B.I == B.LD)
        b = CVectorRange<Number>(B.begin,B.I*B.J);
      // .. otherwise, need to make a copy
      else
        {
        B_copy = B;
        b = CVectorRange<Number>(B_copy.range().begin,B.I*B.J);
        }
      }

    // Presents the columns of B as a (taller) Vector b.
    CVectorRange<Number> range() const
      { return b; }

  private:

    // Original forcing data.
    CMatrixRange<Number> B;

    // A column-contiguous deep copy of B, if necessary.
    Matrix<Number> B_copy;

    // A Vector view of B, aliases either B or B_copy.
    CVectorRange<Number> b;
    
  };

// Captures a Matrix-valued solution data (X), presents it as a Vector (b), either by reshaping or deep copying.
template<class Number> class MrhsReshaperX
  {
  public:

    // Constructor, requires solution data X.
    MrhsReshaperX(const MatrixRange<Number>& in_X)
      : X(in_X)
      {
      // Is X column-contiguous? No copy is needed..
      if (X.I == X.LD)
        x = VectorRange<Number>(X.begin,X.I*X.J);
      // .. otherwise, need to make a copy
      else
        {
        X_copy = X;
        x = VectorRange<Number>(X_copy.range().begin,X.I*X.J);
        }
      }

    // Presents the columns of X as a (taller) Vector x.
    VectorRange<Number> range() const
      { return x; }

    // After linear solution is done, we might need to copy from X_copy back into the original X.
   ~MrhsReshaperX()
      {
      if (X.I != X.LD)
        X.assign(X_copy);
      }
  
  private:
  
    // Original solution data.
    MatrixRange<Number> X;

    // A column-contiguous deep copy of X, if necessary.
    Matrix<Number> X_copy;

    // A Vector view of X, aliases either X or X_copy.
    VectorRange<Number> x;
    
  };

} // namespace detail
} // namespace

#endif