fmatrix.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
// $Id$
#ifndef DUNE_FMATRIX_HH
#define DUNE_FMATRIX_HH

#include <math.h>
#include <complex>
#include <iostream>
#include "exceptions.hh"
#include "fvector.hh"
#include "precision.hh"

namespace Dune {

  /**
              @addtogroup DenseMatVec
              @{
   */

  /*! \file

     \brief  This file implements a matrix constructed from a given type
     representing a field and compile-time given number of rows and columns.
   */

  template<class K, int n, int m> class FieldMatrix;

  /** @brief Error thrown if operations of a FieldMatrix fail. */
  class FMatrixError : public Exception {};

  // template meta program for assignment from scalar
  template<int I>
  struct fmmeta_assignscalar {
    template<class T, class K>
    static void assignscalar (T* x, const K& k)
    {
      fmmeta_assignscalar<I-1>::assignscalar(x,k);
      x[I] = k;
    }
  };
  template<>
  struct fmmeta_assignscalar<0> {
    template<class T, class K>
    static void assignscalar (T* x, const K& k)
    {
      x[0] = k;
    }
  };

  // template meta program for operator+=
  template<int I>
  struct fmmeta_plusequal {
    template<class T>
    static void plusequal (T& x, const T& y)
    {
      x[I] += y[I];
      fmmeta_plusequal<I-1>::plusequal(x,y);
    }
  };
  template<>
  struct fmmeta_plusequal<0> {
    template<class T>
    static void plusequal (T& x, const T& y)
    {
      x[0] += y[0];
    }
  };

  // template meta program for operator-=
  template<int I>
  struct fmmeta_minusequal {
    template<class T>
    static void minusequal (T& x, const T& y)
    {
      x[I] -= y[I];
      fmmeta_minusequal<I-1>::minusequal(x,y);
    }
  };
  template<>
  struct fmmeta_minusequal<0> {
    template<class T>
    static void minusequal (T& x, const T& y)
    {
      x[0] -= y[0];
    }
  };

  // template meta program for operator*=
  template<int I>
  struct fmmeta_multequal {
    template<class T, class K>
    static void multequal (T& x, const K& k)
    {
      x[I] *= k;
      fmmeta_multequal<I-1>::multequal(x,k);
    }
  };
  template<>
  struct fmmeta_multequal<0> {
    template<class T, class K>
    static void multequal (T& x, const K& k)
    {
      x[0] *= k;
    }
  };

  // template meta program for operator/=
  template<int I>
  struct fmmeta_divequal {
    template<class T, class K>
    static void divequal (T& x, const K& k)
    {
      x[I] /= k;
      fmmeta_divequal<I-1>::divequal(x,k);
    }
  };
  template<>
  struct fmmeta_divequal<0> {
    template<class T, class K>
    static void divequal (T& x, const K& k)
    {
      x[0] /= k;
    }
  };

  // template meta program for dot
  template<int I>
  struct fmmeta_dot {
    template<class X, class Y, class K>
    static K dot (const X& x, const Y& y)
    {
      return x[I]*y[I] + fmmeta_dot<I-1>::template dot<X,Y,K>(x,y);
    }
  };
  template<>
  struct fmmeta_dot<0> {
    template<class X, class Y, class K>
    static K dot (const X& x, const Y& y)
    {
      return x[0]*y[0];
    }
  };

  // template meta program for umv(x,y)
  template<int I>
  struct fmmeta_umv {
    template<class Mat, class X, class Y, int c>
    static void umv (const Mat& A, const X& x, Y& y)
    {
      typedef typename Mat::row_type R;
      typedef typename Mat::field_type K;
      y[I] += fmmeta_dot<c>::template dot<R,X,K>(A[I],x);
      fmmeta_umv<I-1>::template umv<Mat,X,Y,c>(A,x,y);
    }
  };
  template<>
  struct fmmeta_umv<0> {
    template<class Mat, class X, class Y, int c>
    static void umv (const Mat& A, const X& x, Y& y)
    {
      typedef typename Mat::row_type R;
      typedef typename Mat::field_type K;
      y[0] += fmmeta_dot<c>::template dot<R,X,K>(A[0],x);
    }
  };

  // template meta program for mmv(x,y)
  template<int I>
  struct fmmeta_mmv {
    template<class Mat, class X, class Y, int c>
    static void mmv (const Mat& A, const X& x, Y& y)
    {
      typedef typename Mat::row_type R;
      typedef typename Mat::field_type K;
      y[I] -= fmmeta_dot<c>::template dot<R,X,K>(A[I],x);
      fmmeta_mmv<I-1>::template mmv<Mat,X,Y,c>(A,x,y);
    }
  };
  template<>
  struct fmmeta_mmv<0> {
    template<class Mat, class X, class Y, int c>
    static void mmv (const Mat& A, const X& x, Y& y)
    {
      typedef typename Mat::row_type R;
      typedef typename Mat::field_type K;
      y[0] -= fmmeta_dot<c>::template dot<R,X,K>(A[0],x);
    }
  };

  template<class K, int n, int m, class X, class Y>
  inline void fm_mmv (const FieldMatrix<K,n,m>& A,  const X& x, Y& y)
  {
    for (int i=0; i<n; i++)
      for (int j=0; j<m; j++)
        y[i] -= A[i][j]*x[j];
  }

  template<class K>
  inline void fm_mmv (const FieldMatrix<K,1,1>& A, const FieldVector<K,1>& x, FieldVector<K,1>& y)