libdune-localfunctions-dev 2.2.1-2 / usr / include / dune / localfunctions / utility / coeffmatrix.hh

#ifndef DUNE_COEFFMATRIX_HH
#define DUNE_COEFFMATRIX_HH
#include <cassert>
#include <iostream>
#include <vector>
#include <dune/common/fvector.hh>
#include <dune/localfunctions/utility/field.hh>
#include <dune/localfunctions/utility/tensor.hh>

namespace Dune 
{
  /*************************************************
   * Default class for storing a coefficient matrix
   * for the PolynomialBasis. Basically a simple
   * CRS structure is used. The additional complexity
   * is due to the storage and efficient evaluation
   * of higher order derivatives. See the remarks
   * in tensor.hh which also hold true for this file.
   *************************************************/
  template <class Field, class Field2>
  struct Mult 
  {
    typedef Field2 BasisEntry;
    static void add(const Field &vec1, const BasisEntry &vec2,
                    BasisEntry &res) 
    {
      res += vec1*vec2;
    }
  };

  template <class Field,class Field2, int dimRange>
  struct Mult< Field,FieldVector<Field2,dimRange> > 
  {
    typedef FieldVector<Field2,dimRange> BasisEntry;
    static void add(const Field &vec1, const BasisEntry &vec2,
                    BasisEntry &res) 
    {
      res.axpy(vec1,vec2);
    }
  };

  template< class F , unsigned int bSize >
  class SparseCoeffMatrix
  {
  public:
    typedef F Field;
    static const unsigned int blockSize = bSize;
    typedef SparseCoeffMatrix<Field,blockSize> This;

    SparseCoeffMatrix() 
    : coeff_(0),
      rows_(0),
      skip_(0),
      numRows_(0),
      numCols_(0)
    {}

    ~SparseCoeffMatrix()
    {
      delete [] coeff_;
      delete [] rows_;
      delete [] skip_;
    }

    const unsigned int size () const
    {
      return numRows_/blockSize;
    }
    const unsigned int baseSize () const
    {
      return numCols_;
    }

    template< class BasisIterator, class FF>
    void mult ( const BasisIterator &x, 
                unsigned int numLsg,
                FF *y ) const
    {
      typedef typename BasisIterator::Derivatives XDerivatives;
      assert( numLsg*blockSize <= (size_t)numRows_ );
      unsigned int row = 0;
      Field *pos = rows_[ 0 ];
      unsigned int *skipIt = skip_;
      XDerivatives val;
      for( size_t i = 0; i < numLsg; ++i)
      {
        for( unsigned int r = 0; r < blockSize; ++r, ++row )
        {
          val = 0;
          BasisIterator itx = x;
          for( ; pos != rows_[ row+1 ]; ++pos, ++skipIt )
          {
            itx += *skipIt;
            val.axpy(*pos,*itx);
          }
          DerivativeAssign<XDerivatives,FF>::apply(r,val,*(y+i*XDerivatives::size*blockSize));
        }
      }
    }
    template< class BasisIterator, class Vector>
    void mult ( const BasisIterator &x, 
                Vector &y ) const
    {
      typedef typename Vector::value_type YDerivatives;
      typedef typename BasisIterator::Derivatives XDerivatives;
      size_t numLsg = y.size();
      assert( numLsg*blockSize <= (size_t)numRows_ );
      unsigned int row = 0;
      Field *pos = rows_[ 0 ];
      unsigned int *skipIt = skip_;
      XDerivatives val;
      for( size_t i = 0; i < numLsg; ++i)
      {
        for( unsigned int r = 0; r < blockSize; ++r, ++row )
        {
          val = 0;
          BasisIterator itx = x;
          for( ; pos != rows_[ row+1 ]; ++pos, ++skipIt )
          {
            itx += *skipIt;
            val.axpy(*pos,*itx);
          }
          DerivativeAssign<XDerivatives,YDerivatives>::apply(r,val,y[i]);
        }
      }
    }
    template <unsigned int deriv, class BasisIterator, class Vector>
    void mult ( const BasisIterator &x, 
                Vector &y ) const
    { 
      typedef typename Vector::value_type YDerivatives;
      typedef typename BasisIterator::Derivatives XDerivatives;
      typedef FieldVector<typename XDerivatives::Field,YDerivatives::dimension> XLFETensor;
      size_t numLsg = y.size();
      assert( numLsg*blockSize <= (size_t)numRows_ );
      unsigned int row = 0;
      Field *pos = rows_[ 0 ];
      unsigned int *skipIt = skip_;
      for( size_t i = 0; i < numLsg; ++i)
      {
        XLFETensor val(typename XDerivatives::Field(0));
        for( unsigned int r = 0; r < blockSize; ++r, ++row )
        {
          BasisIterator itx = x;
          for( ; pos != rows_[ row+1 ]; ++pos, ++skipIt )
          {
            itx += *skipIt;
            LFETensorAxpy<XDerivatives,XLFETensor,deriv>::apply(r,*pos,*itx,val);
          }
        }
        field_cast(val,y[i]);
      }
    }

    template< class RowMatrix >
    void fill ( const RowMatrix &mat, bool verbose=false )
    {
      numRows_ = mat.rows();
      numCols_ = mat.cols();
      unsigned int size = numRows_*numCols_;

      delete [] coeff_;
      delete [] rows_;
      delete [] skip_;

      Field* coeff = new Field[ size ];
      // we always initialize the next skip entry to zero,
      // including the one following the end, so allocate
      // size+1 entries so we will stay within the bounds.
      unsigned int *skip = new unsigned int[ size+1 ];
      rows_ = new Field*[ numRows_+1 ];
      std::vector<Field> row( numCols_ );

      rows_[ 0 ] = coeff;
      Field *cit = coeff;
      unsigned int *sit = skip;
      for( unsigned int r = 0; r < numRows_; ++r )
      {
        *sit = 0;
        mat.row( r, row );
        for( unsigned int c = 0; c < numCols_; ++c ) 
        {
          const Field &val = row[c];
          if (val < Zero<Field>() || Zero<Field>() < val) 
          {
            *cit = val;
            ++sit;
            ++cit;
            *sit = 1;
          } else
          {
            ++(*sit);
          }
        }
        rows_[ r+1 ] = cit;
      }
      assert( size_t(rows_[numRows_]-rows_[0]) <= size_t(size) );
      size = rows_[numRows_]-rows_[0];
      coeff_ = new Field[ size ];
      skip_ = new unsigned int[ size ];
      for (unsigned int i=0;i<size;++i)
      {
        coeff_[i] = coeff[i];
        skip_[i] = skip[i];
      }
      for (unsigned int i=0;i<=numRows_;++i)
        rows_[ i ] = coeff_ + (rows_[ i ] - coeff);
      
      delete [] coeff;
      delete [] skip;

      if (verbose)
        std::cout << "Entries: " << (rows_[numRows_]-rows_[0]) 
                  << " full: " << numCols_*numRows_
                  << std::endl;
    }
    // b += a*C[k]
    template <class Vector>
    void addRow( unsigned int k, const Field &a, Vector &b) const
    {
      assert(k<numRows_);
      unsigned int j=0;
      unsigned int *skipIt = skip_ + (rows_[ k ]-rows_[ 0 ]);
      for( Field *pos = rows_[ k ]; 
           pos != rows_[ k+1 ]; 
           ++pos, ++skipIt )
      {
        j += *skipIt;
        assert( j < b.size() );
        b[j] += field_cast<typename Vector::value_type>( (*pos)*a );  // field_cast
      }
    }
  private:
    SparseCoeffMatrix ( const This &other )
    : numRows_( other.numRows_ ),
      numCols_( other.numCols_ )
    {
      const unsigned int size = other.rows_[numRows_]-other.rows_[0];
      coeff_ = new Field[ size ];
      rows_ = new Field*[ numRows_+1 ];
      skip_ = new unsigned int[ size ];
      for (unsigned int i=0;i<size;++i)
      {
        coeff_[i] = other.coeff_[i];
        skip_[i] = other.skip_[i];
      }
      for (unsigned int i=0;i<=numRows_;++i)
        rows_[ i ] = coeff_ + (other.rows_[ i ] - other.coeff_);
    }

    This &operator= (const This&); 
    Field *coeff_;
    Field **rows_;
    unsigned int *skip_;
    unsigned int numRows_,numCols_;
  };

}

#endif // DUNE_COEFFMATRIX_HH