operator/diagonal.hh

Go to the documentation of this file.

// diagonal matrix and solver for it

#ifndef diagonal_hh
#define diagonal_hh

#include "matrix.hh"

namespace concepts {

  // ******************************************************** DiagonalMatrix **

  template<typename F>
  class DiagonalMatrix : public Matrix<F> {
  public:
    typedef typename Realtype<F>::type r_type;
    typedef typename Cmplxtype<F>::type c_type;

    template<class G>
    DiagonalMatrix(const Space<G>& spc);
    DiagonalMatrix(uint dim);
    template<class G>
    DiagonalMatrix(const Space<G>& spc, const Array<F> entries);
    DiagonalMatrix(const Matrix<F>& matrix);
    virtual ~DiagonalMatrix();

    virtual void operator()(const Function<r_type>& fncY,
                            Function<F>& fncX);
    virtual void operator()(const Function<c_type>& fncY,
                            Function<c_type>& fncX);

    virtual F operator()(const uint i, const uint j) const;
    virtual F& operator()(const uint i, const uint j);

    virtual void transpMult(const Vector<r_type>& fncY, Vector<F>& fncX);
    virtual void transpMult(const Vector<c_type>& fncY,
                            Vector<c_type>& fncX);

    DiagonalMatrix<F>& operator=(F c);

    DiagonalMatrix<F>& operator+=(const Matrix<F>& d);

    template<class H, class I>
    void addInto(Matrix<H>& dest, const I fact);
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    Array<F> entries_;
    F dummy_;

    template<typename H, typename I>
    void apply_(const Function<H>& fncY, Function<I>& fncX);

    template<typename H, typename I>
    void applyT_(const Vector<H>& fncY, Vector<I>& fncX);
  };

  template<typename F>
  template<typename H, typename I>
  void DiagonalMatrix<F>::apply_(const Function<H>& fncY, Function<I>& fncX) {
    conceptsAssert(fncY.dim() == this->dimY(), Assertion());
    conceptsAssert(fncX.dim() == this->dimX(), Assertion());
    const Vector<H>* vecY = dynamic_cast<const Vector<H>*>(&fncY);
    Vector<I>* vecX = dynamic_cast<Vector<I>*>(&fncX);
    if ((vecY != 0) && (vecY != 0)) {
      applyT_(*vecY, *vecX);
      return;
    }

    // dynamic_casts did not work
    for (uint i = 0; i < this->dimX(); ++i)
      fncX(i) = fncY(i) * entries_[i];
  }

  template<typename F>
  template<typename H, typename I>
  void DiagonalMatrix<F>::applyT_(const Vector<H>& fncY, Vector<I>& fncX)
  {
    const H* vecYdata(fncY);
    I* vecXdata(fncX);
    const F* entryData(entries_);
    for (uint i = 0; i < this->dimX(); ++i)
      *vecXdata++ = *vecYdata++ * *entryData++;
  }

  template<typename F>
  template<class H, class I>
  void DiagonalMatrix<F>::addInto(Matrix<H>& dest, const I fact) {
    conceptsAssert(dest.dimX()  == dest.dimY(), Assertion());
    conceptsAssert(this->dimX() == dest.dimX(), Assertion());
    for (uint j = 0; j < this->dimX(); ++j)
      dest(j, j) += entries_[j] * fact;
  }

  // ******************************************************** DiagonalSolver **

  template <typename F>
  class DiagonalSolver : public Operator<F> {
  public:
    typedef typename Realtype<F>::type r_type;
    typedef typename Cmplxtype<F>::type c_type;

    DiagonalSolver(const DiagonalMatrix<F>& m);
    virtual ~DiagonalSolver();

    virtual void operator()(const Function<r_type>& rhs, Function<F>& sol);
    virtual void operator()(const Function<c_type>& rhs,
                            Function<c_type>& sol);
    void operator()(const Matrix<r_type>& mX, Matrix<F>& mY);
    void operator()(const Matrix<c_type>& mX, Matrix<c_type>& mY);

  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    DiagonalMatrix<F> inverse_;

    template <typename H, typename I>
    void apply_(const Matrix<H>& mX, Matrix<I>& mY);
  };

  // This construction allow the solving of real right hand side for
  // both, real and complex diagonal matrices.
  template <typename F>
  void DiagonalSolver<F>::operator()(const Function<r_type>& rhs,
                                     Function<F>& sol) {
    inverse_(rhs, sol);
  }

  // This construction allow the solving of complex right hand side for
  // both, real and complex diagonal matrices.
  template <typename F>
  void DiagonalSolver<F>::operator()(const Function<c_type>& rhs,
                                     Function<c_type>& sol) {
    inverse_(rhs, sol);
  }

  template <typename F>
  template <typename H, typename I>
  void DiagonalSolver<F>::apply_(const Matrix<H>& mX, Matrix<I>& mY) {
    // number of rows
    const uint n = this->dimX();

    // number of rows should be the same as for the solver
    conceptsAssert(mX.dimX() == n, Assertion());
    conceptsAssert(mY.dimX() == n, Assertion());
    // the resulting matrix should have the same number of columns
    // then the applicated matrix
    conceptsAssert(mX.dimY() == mY.dimY(), Assertion());

    for(uint i = 0; i < n; i++) {
      F& d = inverse_(i,i);
      for(uint j = 0; j < n; j++) {
        const H val = mX(i,j);
        if (val != (H)0.0)
          mY(i,j) += d * val;
      }
    }
  }

} // namespace concepts

#endif // diagonal_hh