tmatrix.hh

Go to the documentation of this file.

 
#ifndef tmatrix_hh
#define tmatrix_hh
 
#include <iostream>
#include <map>
#include "toolbox/array.hh"
#include "toolbox/set.hh"
#include "basics/outputOperator.hh"
#include "basics/typedefs.hh"
 
#if __cplusplus >= 201103L
#include <vector>
#endif
 
namespace concepts {
 
  // forward declaration
  template<class F>
  class ElementMatrix;
 
  template<class F>
  class Vector;
 
//   template<class F>
//   class SparseMatrix;  
  
  template<class F, int dim>
  class TColumnBlock;
 
  // *************************************************************** TColumn **
 
  template<class F>
    class TColumn : public virtual OutputOperator {
      public:
        // FIXME: of what use is a TColumn which has no entries, and can never get entries?
        TColumn() : val_(0), lnk_(0), n_(0), idx_(0)  {}
        TColumn(uint n, uint idx, TColumn<F>* lnk = 0, F def = 0);
        TColumn(const TColumn<F>& t);
        virtual ~TColumn();
 
        F operator[](uint i) const { return val_[i]; }
        F& operator[](uint i) { return val_[i]; }
 
        const F* values() const { return (const F*)val_; }
        F* values() { return (F*)val_; }
 
        uint index() const { return idx_; }
 
        TColumn<F>* link() const { return lnk_; }
 
        void append(TColumn<F>* T);
 
        void clear();
 
        uint n() const { return n_; }
 
        virtual TColumn<F>* clone() const;
      protected:
        virtual ::std::ostream& info(::std::ostream& os) const;
        uint successors_() const;
        Array<F> val_;
      private:
        TColumn<F>* lnk_;
 
        uint n_;
 
        uint idx_;
    };
 
  template<class F>
    ::std::ostream& operator<<(::std::ostream& os, TColumn<F>* T) {
      if (T != 0) {
        os << *T;
        if (T->link()) os << ", " << T->link();
      } else os << "TColumn(0)";
      return os;
    }
 
  // ********************************************************* TColumnTensor **
 
  template<class F, int dim>
    class TColumnTensor : public TColumn<F> {
      public:
        TColumnTensor(const uint* n, const uint idx, TColumn<F>* lnk = 0);
        TColumnTensor(const TColumnTensor<F, dim>& t);
        TColumnTensor(const TColumnBlock<F, dim>& t, uint k);
        TColumnTensor(const TColumnTensor<F, dim>& t, const uint* n);
 
        F operator[](uint* i) const {
#ifdef DEBUG
          for (uint j = 0; j < dim; ++j)
            conceptsAssert3(i[j] < n_[j], Assertion(), "j = " << j << ", i[j] = "
                << i[j] << ", n_[j] = " << n_[j]);
#endif
          return TColumn<F>::operator[](
              dim == 2 ? i[0]+i[1]*n_[0] 
              : i[0]+i[1]*n_[0]+i[2]*n_[0]*n_[1]); 
        }
        F& operator[](uint* i) {
#ifdef DEBUG
          for (uint j = 0; j < dim; ++j)
            conceptsAssert3(i[j] < n_[j], Assertion(), "j = " << j << ", i[j] = "
                << i[j] << ", n_[j] = " << n_[j]);
#endif
          return TColumn<F>::operator[](dim == 2 ? i[0]+i[1]*n_[0] :
              i[0]+i[1]*n_[0]+i[2]*n_[0]*n_[1]); }
 
          const uint* n() const { return n_; }
 
          virtual TColumnTensor<F, dim>* clone() const;
      protected:
          virtual ::std::ostream& info(::std::ostream& os) const;
      private: 
          uint n_[dim];
          void recursiveOut_(const int d, uint* n, ::std::ostream& os) const;
    };
 
  // ********************************************************** TColumnBlock **
 
  template<class F, int dim>
    class TColumnBlock : public TColumn<F> {
      public:
        TColumnBlock(const uint* n, const uint idx,
            TColumn<F>* lnk = 0);
        TColumnBlock(const TColumnBlock<F, dim>& t);
        TColumnBlock(const TColumnBlock<F, dim>& t, const uint* n);
 
        F operator[](uint* i) const {
          return TColumn<F>::operator[](ldof(i));
        }
 
        F& operator[](uint* i) {
          return TColumn<F>::operator[](ldof(i));
        }
 
        void write(TColumnTensor<F,dim>& t, uint k);
 
        uint ldof(uint* i) const;
 
        inline const uint* n() const { return n_;} 
        inline const uint n(uint i) const { return n_[i]; }
 
        virtual TColumnBlock<F, dim>* clone() const;
      protected:
        virtual ::std::ostream& info(::std::ostream& os) const;
      private: 
        uint n_[dim*dim];
        void recursiveOut_(const int d, uint* n, ushort& position, ::std::ostream& os) const;
    };
 
  // *********************************************************** TColumnsSet **
 
  template<class F, int dim>
    class TColumnSet : public virtual OutputOperator {
      public:
        // Destructor
        ~TColumnSet();
        bool operator()(uint key);
        TColumn<F>* appendT(uint key, TColumn<F>* T);
        TColumn<F>* T(uint key) const;
        void p(uint key, ushort p[dim]);
        void p(uint key, ushort p);
        Array<ushort> p(uint key) const;
        void clear();
      protected:
        virtual ::std::ostream& info(::std::ostream& os) const;
      private:
        ::std::map<uint, TColumn<F>*> T_;
        ::std::map<uint, Array<ushort> > p_;
    };
 
  // *********************************************************** TMatrixBase **
 
  template<class F> 
    class TMatrixBase : public virtual OutputOperator {
      public:
        TMatrixBase(const uint n);
        TMatrixBase(const uint m, const uint n);
        virtual ~TMatrixBase();
 
        virtual void operator()(const ElementMatrix<F>& A, 
                                ElementMatrix<F>& B) const = 0;
        virtual void operator()(const ElementMatrix<::std::complex<F> >& A, 
                                ElementMatrix<::std::complex<F> >& B) const = 0;
 
        virtual uint index(const uint i) const = 0;
 
        inline uint m() const { return m_; }
        
        inline void setM(uint m) { m_ = m; }
 
        inline uint n() const { return n_; }
 
        virtual void usedIdx(TColumn<bool>& c) const = 0;
 
        virtual void extract(const concepts::Vector<F>& solution,
                             concepts::Array<F>& coeff) const = 0;
        virtual void extract(const concepts::Vector<::std::complex<F> >& solution,
                             concepts::Array<::std::complex<F> >& coeff) const = 0;
      protected:
        virtual ::std::ostream& info(::std::ostream& os) const;
 
        uint n_;
        uint m_;
    };
 
  // **************************************************************** TIndex **
 
  template<class F> 
    class TIndex : public TMatrixBase<F> {
      public:
        TIndex(const uint m, const uint n, uint idx[]);
 
        TIndex(const uint m, const uint n, const Array<uint>& idx);
        virtual ~TIndex();
 
        void operator()(const ElementMatrix<F>& A, ElementMatrix<F>& B) const;
        void operator()(const ElementMatrix<::std::complex<F> >& A,
                        ElementMatrix<::std::complex<F> >& B) const;
 
        inline uint index(const uint i) const { return idx_[i]; }
 
        virtual void usedIdx(TColumn<bool>& c) const;
        virtual void extract(const concepts::Vector<F>& solution,
                             concepts::Array<F>& coeff) const;
        virtual void extract(const concepts::Vector<::std::complex<F> >& solution,
                             concepts::Array<::std::complex<F> >& coeff) const;
      protected:
        virtual ::std::ostream& info(::std::ostream& os) const;
      private:
        Array<uint> idx_;
    };
 
  // *************************************************************** TMatrix **
 
  template<class F>
  class TMatrix : public TMatrixBase<F> { 
  public:
    struct Control : public virtual OutputOperator {
      uint idx;
      uint sz;
    protected:
      virtual ::std::ostream& info(::std::ostream& os) const;
    };
 
    struct Data : public virtual OutputOperator {
      uint idx;
      F data;
    protected:
      virtual ::std::ostream& info(::std::ostream& os) const;
    };
 
    TMatrix(TColumn<F>* T);
    
#if __cplusplus >= 201103L
 
    TMatrix(const uint nShapeFunctions,
            const ::std::vector<::std::pair<uint,uint>>& ids,
            const ::std::vector<F>& coeffs = ::std::vector<F>(0)  );    
#endif // __cplusplus >= 201103L    
    
    TMatrix(const TMatrix<F>& T, Set<uint>* idx = 0);
        
    virtual ~TMatrix();
 
    template<class G>
    TMatrix(const TMatrix<G>& T);
    void operator()(const ElementMatrix<F>& A, ElementMatrix<F>& B) const;
    void operator()(const ElementMatrix<::std::complex<F> >& A,
                    ElementMatrix<::std::complex<F> >& B) const;
//     void operator()(const SparseMatrix<F>& A,
//                     SparseMatrix<F>& B) const;                   
    void operator()(const Vector<F>& A,
                    Vector<F>& B) const;                    
    void applyOne(const ElementMatrix<F>& A, ElementMatrix<F>& B) const;
 
    virtual uint index(const uint i) const { return ctrl_[i].idx; }
 
    void append(TColumn<F>* T);
 
    const typename TMatrix::Control* control(uint i) const { return &ctrl_[i]; };
 
    const typename TMatrix::Data* data(uint i) const { return data_+i; };
 
    // Returns the T-Matrix as linked list of columns
    TColumn<F>* columns() const;
 
    uint nEntries() const {
      return sz_;
    }
 
    virtual void usedIdx(TColumn<bool>& c) const;
    virtual void usedIdx(Set<uint>& c) const;
    virtual void extract(const concepts::Vector<F>& solution,
                         concepts::Array<F>& coeff) const;
    virtual void extract(const concepts::Vector<::std::complex<F> >& solution,
                         concepts::Array<::std::complex<F> >& coeff) const;
//     Vector<F> extractToVector(const Vector<F>& solution) const;                         
    void extractOne(const concepts::Vector<F>& solution,
                    concepts::Array<F>& coeff) const;
#if __cplusplus >= 201103L
    void concatenate(const TMatrix<F>& other);
#endif // __cplusplus >= 201103L
    
    ::std::ostream& outputLocalMap(::std::ostream& os) const;
  protected:
    virtual ::std::ostream& info(::std::ostream& os) const;
  private:
    Array<Control> ctrl_;
 
    Array<Data> data_;
 
    uint sz_;
  };
 
  template<class F>
  template<class G>
    TMatrix<F>::TMatrix(const TMatrix<G>& T)
    : TMatrixBase<F>(T.n())
      , ctrl_(T.n())
      , data_(T.nEntries())
      , sz_(T.nEntries())
    {
      setM( T.m() );
 
      for(uint i=0; i < this->n(); ++i) {
        ctrl_[i].idx = T.control(i)->idx;
        ctrl_[i].sz = T.control(i)->sz;
      }
      for(uint i=0; i < sz_; ++i) {
        data_[i].idx = T.data(i)->idx;
        data_[i].data = T.data(i)->data;
      }
 
    }
 
 
  // *************************************************************** TMatrixBlock **
 
    template<class F>
class TMatrixBlock: public TMatrix<F> {
public:
  TMatrixBlock(const TMatrix<F>* T0, F weight0 = 1.0, const TMatrix<F>* T1 = 0, F weight1 = 0.0);
 
  virtual ~TMatrixBlock(){};
 
  virtual void usedIdx(TColumn<bool>& c) const{
    TMatrix<F>::usedIdx(c);
  }
 
  virtual void usedIdx(Set<uint>& c) const{
    TMatrix<F>::usedIdx(c);
  }
  virtual void extract(const concepts::Vector<F>& solution,
      concepts::Array<F>& coeff) const{
    TMatrix<F>::extract(solution,coeff);
  }
  virtual void extract(const concepts::Vector<::std::complex<F> >& solution,
      concepts::Array<::std::complex<F> >& coeff) const{
    TMatrix<F>::extract(solution, coeff);
  }
 
  inline uint offsetrow() const{
    return offsetrow_;
  }
 
protected:
  virtual ::std::ostream& info(::std::ostream& os) const;
 
private:
  //routine append TMatrix information at a offset position of row (as longer TColumns)
  void insert_(const TMatrix<F>& T,uint offset, uint size, F weight);
 
  //information about the offset characterising the blocks
  uint offsetrow_;
  uint offsetcol_;
};
 
 
 
 
} // namespace concepts
 
#endif // tmatrix_hh