// $Id: dematrix.h,v 1.1.1.1 2000/11/02 08:47:16 fritzi Exp $

// This file is part of the DEAL Library
// DEAL is Copyright(1995) by
// Roland Becker, Guido Kanschat, Franz-Theo Suttmeier


#ifndef __dematrix_h
#define __dematrix_h

/*
CLASS
   dEMatrix
   */
template <class NUM, class VECTOR>
class dEMatrix
{
protected:
  Triangulation& tr;
  int dim_im, dim_rg;
public:
  //////////
  dEMatrix(int m, int n, Triangulation& t) : tr(t), dim_im(m), dim_rg(n) {}
  //////////
  dEMatrix(int n, Triangulation& t) : tr(t), dim_im(n), dim_rg(n) {}
  
  //////////
  void reinit(int n) { dim_im = dim_rg = n; }
  //////////
  void reinit(int m, int n) { dim_im = m; dim_rg = n; }

  //////////
  int m() const { return dim_im; }
  //////////
  int n() const { return dim_rg; }

  //////////
  void vmult(VECTOR& dst, const VECTOR& src)
  {
    dst = 0.;
    CellPtr<NUM> c;
    for (c=tr.first_cell(); c.cell() ;c=tr.next_cell(c))
    {
      if (!c.cell()->sleep()) c->vmult(dst,src);
    }
  }

  //////////
  void Tvmult(VECTOR& dst, const VECTOR& src)
  {
    dst = 0.;
    CellPtr<NUM> c;
    for (c=tr.first_cell(); c.cell() ;c=tr.next_cell(c))
    {
      if (!c.cell()->sleep()) c->Tvmult(dst,src);
    }
  }
  //////////
  double residual(VECTOR& dst, const VECTOR& src, const VECTOR& b)
  {
    dst=0.;
    vmult(dst,src);
    dst.sadd(-1.,1.,b);
    return sqrt(dst*dst);
  }
};

#define EMATRIX_DECL(VECTOR) \
void vmult(VECTOR& dst, const VECTOR& src);\
void Tvmult(VECTOR& dst, const VECTOR& src);


#endif