#include <list>

#include "FixVec.h"
#include "Boundary.h"
#include "DOFAdmin.h"
#include "ElInfo.h"
#include "Error.h"
#include "FiniteElemSpace.h"
#include "Global.h"
#include "Mesh.h"
#include "Quadrature.h"
#include "AbstractFunction.h"
#include "BoundaryManager.h"
#include "ElementVector.h"
#include "Assembler.h"

namespace AMDiS {

  // template<typename T>
  // void DOFVector<T>::coarseRestrict(RCNeighbourList& list, int n) {}

  template<typename T>
  void DOFVectorBase<T>::addElementVector(T factor, 
					  const ElementVector &elVec, 
					  const BoundaryType *bound,
					  bool add)
  {
    FUNCNAME("DOFVector::addElementVector");
    DegreeOfFreedom i, irow;

    int n_row = elVec.getSize();

    for (i = 0; i < n_row; i++) {
      BoundaryCondition *condition = 
	bound ? this->getBoundaryManager()->getBoundaryCondition(bound[i]) : NULL;

      if(!(condition && condition->isDirichlet())) {
	irow = elVec.dofIndices[i];
	(*this)[irow] = (add ? (*this)[irow] : 0.0);
	(*this)[irow] += factor * elVec[i];
      }
    }
  }

  template<typename T>
  DOFVector<T>::DOFVector(const FiniteElemSpace* f,::std::string n)
    : DOFVectorBase<T>(f, n),
      //elementVector(NULL), 
      refineInter(false), coarsenOperation(NO_OPERATION)
  {
    init(f, n);
  } 

  template<typename T>
  void DOFVector<T>::init(const FiniteElemSpace* f,::std::string n)
  {
    name = n;
    feSpace = f;
    if(feSpace && feSpace->getAdmin()) {
      (feSpace->getAdmin())->addDOFIndexed(this);
    }  
    this->boundaryManager = NEW BoundaryManager;
  }

  template<typename T>
  DOFVector<T>::~DOFVector()
  {
    if(feSpace && feSpace->getAdmin()) {
      (feSpace->getAdmin())->removeDOFIndexed(this);
    }

    if (this->boundaryManager) {
      DELETE this->boundaryManager;
    }
  }

  template<typename T>
  DOFVector<T> * DOFVector<T>::traverseVector = NULL;

  template<typename T>
  FastQuadrature *DOFVector<T>::quad_fast = NULL;

  template<typename T>
  double DOFVector<T>::norm = 0.0;

  template<typename T>
  int DOFVector<T>::dim = 0;

  template<typename T>
  double DOFVector<T>::nrm2() const
  {
    FUNCNAME("DOFVector<T>::nrm2()");
    double nrm;
    const DOFAdmin *admin = NULL;

    TEST_EXIT(feSpace && (admin = feSpace->getAdmin()))
      ("pointer is NULL: %8X, %8X\n", feSpace,admin);
    TEST_EXIT(static_cast<int>(vec.size()) >= admin->getUsedSize())
      ("size = %d too small: admin->sizeUsed = %d\n", vec.size(),
       admin->getUsedSize());

    nrm = 0.0;

    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), USED_DOFS);
    for(vecIterator.reset(); !vecIterator.end(); ++vecIterator)
      nrm += (*vecIterator) * (*vecIterator);

    return(sqrt(nrm));
  }

  template<typename T>
  double DOFVector<T>::squareNrm2() const
  {
    FUNCNAME("DOFVector<T>::nrm2()");
    double nrm;
    const DOFAdmin *admin = NULL;

    TEST_EXIT(feSpace && (admin = feSpace->getAdmin()))
      ("pointer is NULL: %8X, %8X\n", feSpace,admin);
    TEST_EXIT(static_cast<int>(vec.size()) >= admin->getUsedSize())
      ("size = %d too small: admin->sizeUsed = %d\n", vec.size(),
       admin->getUsedSize());

    nrm = 0.0;

    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), USED_DOFS);
    for(vecIterator.reset(); !vecIterator.end(); ++vecIterator)
      nrm += (*vecIterator) * (*vecIterator);

    return nrm;
  }

  template<typename T>
  T DOFVector<T>::asum() const
  {
    FUNCNAME("DOFVector<T>::asum");
    double    nrm;
    const DOFAdmin *admin = NULL;
    TEST_EXIT(feSpace && (admin = feSpace->getAdmin()))
      ("pointer is NULL: %8X, %8X\n", feSpace, admin);
    TEST_EXIT(vec.size() >= admin->getUsedSize())
      ("size = %d too small: admin->sizeUsed = %d\n", vec.size(), 
       admin->getUsedSize());

    nrm = 0.0;

    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), USED_DOFS);
    for(vecIterator.reset(); !vecIterator.end(); ++vecIterator)
      nrm += abs(*vecIterator);

    return(nrm);
  }

  // Andreas ergaenzte ..

  template<typename T>
  T DOFVector<T>::sum() const
  {
    FUNCNAME("DOFVector<T>::sum");
    double    nrm;
    const DOFAdmin *admin = NULL;
    TEST_EXIT(feSpace && (admin = feSpace->getAdmin()))
      ("pointer is NULL: %8X, %8X\n", feSpace, admin);
    TEST_EXIT(vec.size() >= admin->getUsedSize())
      ("size = %d too small: admin->sizeUsed = %d\n", vec.size(), 
       admin->getUsedSize());

    nrm = 0.0;

    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), USED_DOFS);
    for(vecIterator.reset(); !vecIterator.end(); ++vecIterator)
      nrm += *vecIterator;

    return(nrm);
  }

  // .. bis hier (Andreas)


  template<typename T>
  void DOFVector<T>::set(T alpha)
  {
    FUNCNAME("DOFVector<T>::set");
    const DOFAdmin *admin = NULL;

    TEST_EXIT(feSpace && (admin = feSpace->getAdmin()))
      ("pointer is NULL: %8X, %8X\n", feSpace, admin);
    TEST_EXIT(static_cast<int>( vec.size()) >= admin->getUsedSize())
      ("size = %d too small: admin->sizeUsed = %d\n", vec.size(),
       admin->getUsedSize());

    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(this), USED_DOFS);
    for(vecIterator.reset(); !vecIterator.end(); ++vecIterator) {
      *vecIterator = alpha ; 
    };
  }

  // template<typename T>
  // void DOFVector<T>::scal(T alpha)
  // {
  //   FUNCNAME("DOFVector<T>::scal");
  //   const DOFAdmin *admin = NULL;

  //   TEST_EXIT( feSpace && (admin = feSpace->getAdmin()))
  //     ("pointer is NULL: %8X, %8X\n", this, admin);
  //   TEST_EXIT(static_cast<int>(vec.size()) >= admin->getUsedSize())
  //     ("size = %d too small: admin->sizeUsed = %d\n", vec.size(), 
  //      admin->getUsedSize());

  //   Iterator vecIterator((DOFIndexed<T>*) this, USED_DOFS);
  //   for(vecIterator.reset(); !vecIterator.end(); ++vecIterator) {
  //     (*vecIterator) *= alpha ; 
  //   };
  // }


  template<typename T>
  void DOFVector<T>::copy(const DOFVector<T>& x)
  {
    FUNCNAME("DOFVector<T>::copy");
    const DOFAdmin *admin = NULL;

    TEST_EXIT(feSpace && x.feSpace)
      ("feSpace is NULL: %8X, %8X\n", feSpace,x.feSpace);
    TEST_EXIT((admin = feSpace->getAdmin()) && (admin == x.feSpace->getAdmin()))
      ("no admin or different admins: %8X, %8X\n",
       feSpace->getAdmin(), x.feSpace->getAdmin());
    TEST_EXIT(static_cast<int>(vec.size()) >= admin->getUsedSize())
      ("size = %d too small: admin->sizeUsed = %d\n", vec.size(), 
       admin->getUsedSize());
    TEST_EXIT(static_cast<int>(x.vec.size()) >= admin->getUsedSize())
      ("x.size = %d too small: admin->sizeUsed = %d\n", x.vec.size(), 
       admin->getUsedSize());

    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(this), USED_DOFS);
    Iterator xIterator(dynamic_cast<DOFVector<T>*>(const_cast<DOFVector<T>*>(&x)), USED_DOFS);
    for(vecIterator.reset(), xIterator.reset();
	!vecIterator.end();
	++vecIterator, ++xIterator) 
      {  
	*vecIterator = *xIterator; 
      };
  }

  // template<typename T>
  // void DOFVector<T>::axpy(T alpha, const DOFVector<T>& x)
  // {
  //   FUNCNAME("DOFVector<T>::axpy");
  //   const DOFAdmin *admin;

  //   TEST_EXIT(feSpace && x.feSpace)
  //     ("feSpace is NULL: %8X, %8X\n", feSpace,x.feSpace);
  //   TEST_EXIT((admin = feSpace->getAdmin()) && (admin == x.feSpace->getAdmin()))
  //     ("no admin or different admins: %8X, %8X\n",
  //      feSpace->getAdmin(), x.feSpace->getAdmin());
  //   TEST_EXIT(static_cast<int>(vec.size()) >= admin->getUsedSize())
  //     ("size = %d too small: admin->size = %d\n", vec.size(), 
  //      admin->getUsedSize());
  //   TEST_EXIT(static_cast<int>(x.vec.size()) >= admin->getUsedSize())
  //     ("x.size = %d too small: admin->size = %d\n", x.vec.size(), 
  //      admin->getUsedSize());

  //   Iterator vecIterator((DOFIndexed<T>*) this, USED_DOFS);
  //   Iterator xIterator((DOFIndexed<T>*) &x, USED_DOFS);
  //   for(vecIterator.reset(), xIterator.reset();
  //       !vecIterator.end();
  //       ++vecIterator, ++xIterator) 
  //   {  
  //     *vecIterator += alpha * (*xIterator); 
  //   };
  // }

  // template<typename T>
  // void DOFVector<T>::axpy(T alpha, 
  // 			const DOFVector<T>& x, 
  // 			const DOFVector<T>& y)
  // {
  //   FUNCNAME("DOFVector<T>::axpy");
  //   const DOFAdmin *admin;

  //   TEST_EXIT(feSpace && x.feSpace && y.feSpace)
  //     ("feSpace is NULL: %8X, %8X, %8X\n", feSpace, x.feSpace, y.feSpace);
  //   TEST_EXIT((admin = feSpace->getAdmin()) && 
  // 	    (admin == x.feSpace->getAdmin()) &&
  // 	    (admin == y.feSpace->getAdmin()))
  //     ("no admin or different admins: %8X, %8X, %8X\n",
  //      feSpace->getAdmin(), x.feSpace->getAdmin(), y.feSpace->getAdmin());
  //   TEST_EXIT(static_cast<int>(vec.size()) >= admin->getUsedSize())
  //     ("size = %d too small: admin->size = %d\n", vec.size(), 
  //      admin->getUsedSize());
  //   TEST_EXIT(static_cast<int>(x.vec.size()) >= admin->getUsedSize())
  //     ("x.size = %d too small: admin->size = %d\n", x.vec.size(), 
  //      admin->getUsedSize());
  //   TEST_EXIT(static_cast<int>(y.vec.size()) >= admin->getUsedSize())
  //     ("y.size = %d too small: admin->size = %d\n", y.vec.size(), 
  //      admin->getUsedSize());

  //   Iterator vecIterator((DOFIndexed<T>*) this, USED_DOFS);
  //   Iterator xIterator((DOFIndexed<T>*) &x, USED_DOFS);
  //   Iterator yIterator((DOFIndexed<T>*) &y, USED_DOFS);
  //   for(vecIterator.reset(), xIterator.reset(), yIterator.reset();
  //       !vecIterator.end();
  //       ++vecIterator, ++xIterator, ++yIterator) 
  //   {  
  //     *vecIterator = alpha * (*xIterator) + (*yIterator); 
  //   };
  // }

  // template<typename T>
  // void DOFVector<T>::xpay(T alpha, const DOFVector<T>& x)
  // {
  //   FUNCNAME("DOFVector<T>::xpay");
  //   const DOFAdmin *admin;

  //   TEST_EXIT(feSpace && x.feSpace)
  //     ("feSpace is NULL: %8X, %8X\n", feSpace,x.feSpace);
  //   TEST_EXIT((admin = feSpace->getAdmin()) && (admin == x.feSpace->getAdmin()))
  //     ("no admin or different admins: %8X, %8X\n",
  //      feSpace->getAdmin(), x.feSpace->getAdmin());
  //   TEST_EXIT(static_cast<int>(vec.size()) >= admin->getUsedSize())
  //     ("size = %d too small: admin->sizeUsed = %d\n",
  //      vec.size(), admin->getUsedSize());
  //   TEST_EXIT(static_cast<int>(x.vec.size()) >= admin->getUsedSize())
  //     ("x.size = %d too small: admin->sizeUsed = %d\n", 
  //      x.vec.size(), admin->getUsedSize());

  //   Iterator vecIterator((DOFIndexed<T>*) this, USED_DOFS);
  //   Iterator xIterator((DOFIndexed<T>*) &x, USED_DOFS);
  //   for(vecIterator.reset(), xIterator.reset();
  //       !vecIterator.end();
  //       ++vecIterator, ++xIterator) 
  //   {
  //     *vecIterator = (*xIterator) + alpha * (*vecIterator); 
  //   };
  // }

  template<typename T>
  T DOFVector<T>::min() const
  {
    FUNCNAME("DOFVector<T>::min");
    T m;
    const DOFAdmin *admin = NULL;

    TEST_EXIT( feSpace && (admin = feSpace->getAdmin()))
      ("pointer is NULL: %8X, %8X\n", this, admin);
    TEST_EXIT((static_cast<int>( vec.size())) >= admin->getUsedSize())
      ("size = %d too small: admin->sizeUsed = %d\n", vec.size(), 
       admin->getUsedSize());
  
    Iterator vecIterator(const_cast<DOFIndexed<T>*>(dynamic_cast<const DOFIndexed<T>*>(this)), USED_DOFS);
    for(vecIterator.reset(), m = *vecIterator; !vecIterator.end(); ++vecIterator)
      {
	m = ::std::min(m, *vecIterator);
      }
    return m;
  }

  template<typename T>
  T DOFVector<T>::max() const 
  {
    FUNCNAME("DOFVector<T>::max");
    T m;
    const DOFAdmin *admin = NULL;

    TEST_EXIT( feSpace && (admin = feSpace->getAdmin()))
      ("pointer is NULL: %8X, %8X\n", this, admin);
    TEST_EXIT((static_cast<int>( vec.size())) >= admin->getUsedSize())
      ("size = %d too small: admin->sizeUsed = %d\n", vec.size(), 
       admin->getUsedSize());

    Iterator vecIterator(const_cast<DOFIndexed<T>*>(dynamic_cast<const DOFIndexed<T>*>(this)), USED_DOFS);
    for(vecIterator.reset(), m = *vecIterator; !vecIterator.end(); ++vecIterator)
      {
	m = ::std::max(m, *vecIterator);
      }
    return m;
  }

  template<typename T>
  void gemv(MatrixTranspose transpose, T alpha,
	    const DOFMatrix& a, const DOFVector<T>& x,
	    T beta, DOFVector<T>& y)
  {
    FUNCNAME("gemv");
    int        j, jcol, ysize;
    T          sum, ax;
    const DOFMatrix::MatrixRow *row;
    const DOFAdmin *cadmin,*radmin;

    TEST_EXIT(a.getRowFESpace() && a.getColFESpace() && x.getFESpace() && y.getFESpace())
      ("feSpace is NULL: %8X, %8X, %8X, %8X\n", a.getRowFESpace(),a.getColFESpace(),x.getFESpace(),y.getFESpace());
    TEST_EXIT((radmin=a.getRowFESpace()->getAdmin()) && 
	      (cadmin = a.getColFESpace()->getAdmin()) && 
	      (((transpose == NoTranspose) && (cadmin == x.getFESpace()->getAdmin()) && 
		(radmin == y.getFESpace()->getAdmin()))||
	       ((transpose == Transpose) && (radmin == x.getFESpace()->getAdmin()) && 
		(cadmin == y.getFESpace()->getAdmin()))))
      ("no admin or different admins: %8X, %8X, %8X, %8X\n",
       a.getRowFESpace()->getAdmin(), a.getColFESpace()->getAdmin(), 
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
  
    if (transpose == NoTranspose) {
      TEST_EXIT(static_cast<int>(x.getSize()) >= cadmin->getUsedSize())
	("x.size = %d too small: admin->sizeUsed = %d\n",
	 x.getSize(), cadmin->getUsedSize());
      TEST_EXIT(static_cast<int>(y.getSize()) >= radmin->getUsedSize())
	("y.size = %d too small: admin->sizeUsed = %d\n",
	 y.getSize(), radmin->getUsedSize());
      TEST_EXIT(static_cast<int>( a.getSize()) >= radmin->getUsedSize())
	("a.size = %d too small: admin->sizeUsed = %d\n",
	 a.getSize(), radmin->getUsedSize());
    }
    else  if (transpose == Transpose) {
      TEST_EXIT(static_cast<int>(x.getSize()) >= radmin->getUsedSize())
	("x.size = %d too small: admin->sizeUsed = %d\n",
	 x.getSize(), radmin->getUsedSize());
      TEST_EXIT(static_cast<int>(y.getSize()) >= cadmin->getUsedSize())
	("y.size = %d too small: admin->sizeUsed = %d\n",
	 y.getSize(), cadmin->getUsedSize());
      TEST_EXIT(static_cast<int>( a.getSize()) >= radmin->getUsedSize())
	("a.size = %d too small: admin->sizeUsed = %d\n",
	 a.getSize(), radmin->getUsedSize());
    }

    ysize = y.getSize();

    typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&y), FREE_DOFS);
    for(vecIterator.reset();
	!vecIterator.end(); ++vecIterator) { 
      *vecIterator = 0;
    };

    if (transpose == NoTranspose) {
      typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&y), USED_DOFS); 
      DOFMatrix::Iterator rowIterator(const_cast<DOFMatrix*>(&a), USED_DOFS);
      for(vecIterator.reset(), rowIterator.reset();
	  !rowIterator.end();
	  ++rowIterator, ++vecIterator) { 
	sum = 0;
	row = &(a[rowIterator.getDOFIndex()]);
	for(::std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
	    colIterator != rowIterator->end();
	    colIterator++) {
	  jcol = colIterator->col;
	  if (jcol >= 0) { // entry used? 
	    sum += (static_cast<T>(colIterator->entry)) * x[jcol];
	  } else {
	    if (jcol == DOFMatrix::NO_MORE_ENTRIES)
	      break;
	  }
	}
	*vecIterator *= beta;
	*vecIterator += alpha * sum;
      };
    } else if (transpose == Transpose) {
      typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&y), USED_DOFS);
      for(vecIterator.reset();
	  !vecIterator.end(); 
	  ++vecIterator) {
	*vecIterator  *= beta ; 
      };
    
      typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&x)), USED_DOFS);
      DOFMatrix::Iterator rowIterator(const_cast<DOFMatrix*>(&a), USED_DOFS);
      for(xIterator.reset(), rowIterator.reset();
	  !rowIterator.end();
	  ++rowIterator, ++xIterator) { 
	ax = alpha * (*xIterator);
	row = &(a[rowIterator.getDOFIndex()]);
	for(::std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
	    colIterator != rowIterator->end();
	    colIterator++) {
	  jcol = colIterator->col;
	  if (jcol >= 0) // entry used?
	    y[jcol] += ax * (static_cast<T>(colIterator->entry));
	  else 
	    if (jcol == DOFMatrix::NO_MORE_ENTRIES) break;
	}
      }
    }
    else {
      ERROR_EXIT("transpose=%d\n", transpose);
    }
  }


  // /****************************************************************************/

  // template<typename T>
  // void DOFVector<T>::mv(MatrixTranspose transpose,
  // 		      const DOFMatrix&a, const DOFVector<T>&x)
  // {
  //   FUNCNAME("DOFVector<T>::mv");
  //   int            jcol;
  //   T              sum, ax;
  //   const DOFAdmin *radmin,*cadmin = NULL;

  //   TEST_EXIT(a.getRowFESpace() && a.getColFESpace() && x.feSpace && feSpace)
  //     ("feSpace is NULL: %8X, %8X, %8X, %8X\n", 
  //      a.getRowFESpace(), a.getColFESpace(), x.feSpace,feSpace);
  //   TEST_EXIT((radmin=a.getRowFESpace()->getAdmin()) && 
  // 	    (cadmin = a.getColFESpace()->getAdmin()) && 
  // 	    (((transpose == NoTranspose) && (cadmin == x.feSpace->getAdmin()) && 
  // 	    (radmin == feSpace->getAdmin()))||
  // 	     ((transpose == Transpose) && (radmin == x.feSpace->getAdmin()) && 
  // 	      (cadmin == feSpace->getAdmin()))))
  //     ("no admin or different admins: %8X, %8X, %8X, %8X\n",
  //      a.getRowFESpace()->getAdmin(), a.getColFESpace()->getAdmin(), 
  //      x.feSpace->getAdmin(), feSpace->getAdmin());

  //   if (transpose == NoTranspose) {
  //     TEST_EXIT(static_cast<int>(x.vec.size()) >= cadmin->getUsedSize())
  //       ("x.size = %d too small: admin->sizeUsed = %d\n",
  //        x.vec.size(), cadmin->getUsedSize());
  //     TEST_EXIT(static_cast<int>(vec.size()) >= radmin->getUsedSize())
  //       ("size = %d too small: admin->sizeUsed = %d\n",
  //        vec.size(), radmin->getUsedSize());
  //     TEST_EXIT(static_cast<int>( a.getSize()) >= radmin->getUsedSize())
  //       ("a.size = %d too small: admin->sizeUsed = %d\n",
  //        a.getSize(), radmin->getUsedSize());

  //     Iterator vecIterator((DOFIndexed<T>*) this, USED_DOFS); 
  //     DOFMatrix::Iterator rowIterator((DOFMatrix*) &a, USED_DOFS);
  //     for(vecIterator.reset(), rowIterator.reset();
  // 	!rowIterator.end();
  // 	++rowIterator, ++vecIterator) { 
  //       sum = 0;
  //       for(::std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
  // 	  colIterator != rowIterator->end();
  // 	  colIterator++) {
  // 	jcol = colIterator->col;
  // 	if (jcol >= 0) { // entry used? 
  // 	  sum += ((T) colIterator->entry) * x[jcol];
  // 	} else {
  // 	  if (jcol == DOFMatrix::NO_MORE_ENTRIES)
  // 	    break;
  // 	}
  //       }
  //       *vecIterator = sum;
  //     }
  //   } else if (transpose == Transpose) {
  //     TEST_EXIT(static_cast<int>(x.vec.size()) >= radmin->getUsedSize())
  //       ("x.size = %d too small: admin->sizeUsed = %d\n",
  //        x.vec.size(), radmin->getUsedSize());
  //     TEST_EXIT(static_cast<int>(vec.size()) >= cadmin->getUsedSize())
  //       ("size = %d too small: admin->sizeUsed = %d\n",
  //        vec.size(), cadmin->getUsedSize());
  //     TEST_EXIT(static_cast<int>( a.getSize()) >= radmin->getUsedSize())
  //       ("a.size = %d too small: admin->sizeUsed = %d\n",
  //        a.getSize(), radmin->getUsedSize());

  //     Iterator vecIterator((DOFIndexed<T>*) this, USED_DOFS);
  //     for(vecIterator.reset();
  // 	!vecIterator.end(); ++vecIterator) {
  //       *vecIterator = 0; 
  //     };

  //     Iterator xIterator((DOFIndexed<T>*) &x, USED_DOFS);
  //     DOFMatrix::Iterator rowIterator((DOFMatrix*) &a, USED_DOFS);
  //     for(xIterator.reset(), rowIterator.reset();
  // 	!rowIterator.end();
  // 	++rowIterator, ++xIterator) { 
  //       ax = (*xIterator);
  //       for(::std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
  // 	  colIterator != rowIterator->end();
  // 	  colIterator++) {
  // 	jcol = colIterator->col;
  // 	if (jcol >= 0) // entry used?
  // 	  vec[jcol] += ax * ((T) colIterator->entry);
  // 	else 
  // 	  if (jcol == DOFMatrix::NO_MORE_ENTRIES) break;
  //       }
  //     }
  //   } else {
  //     ERROR_EXIT("transpose=%d\n", transpose);
  //   }
  // }

  template<typename T>
  void DOFVector<T>::print() const
  {
    FUNCNAME("DOFVector<T>::print");
    int  i, j;
    const DOFAdmin *admin = NULL;
    const char     *format;

    if (feSpace) admin = feSpace->getAdmin();

    MSG("Vec `%s':\n", name.c_str());
    j = 0;
    if (admin) {
      if (admin->getUsedSize() > 100)
	format = "%s(%3d,%10.5le)";
      else if (admin->getUsedSize() > 10)
	format = "%s(%2d,%10.5le)";
      else
	format = "%s(%1d,%10.5le)";

      Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), USED_DOFS);
      for(vecIterator.reset();
	  !vecIterator.end(); ++vecIterator) {
	if ((j % 3) == 0) {
	  if (j) Msg::print("\n");
	  MSG(format, "", vecIterator.getDOFIndex(), *vecIterator);
	}
	else 
	  Msg::print(format, " ", vecIterator.getDOFIndex(), *vecIterator);
	j++;
      };
      Msg::print("\n");
    }
    else
      {
	MSG("no DOFAdmin, print whole vector.\n");
    
	for (i = 0; i < static_cast<int>( vec.size()); i++) {
	  if ((j % 3) == 0)
	    {
	      if (j) Msg::print("\n");
	      MSG("(%d,%10.5e)",i,vec[i]);
	    }
	  else 
	    Msg::print(" (%d,%10.5e)",i,vec[i]);
	  j++;
	}
	Msg::print("\n");
      }
    return;
  }


  template<typename T>
  T DOFVectorBase<T>::evalUh(const DimVec<double>& lambda, 
			     DegreeOfFreedom* dof_indices)
  {
    int i;
    BasisFunction* phi = const_cast<BasisFunction*>(this->getFESpace()->getBasisFcts());
    int numberOfBasFcts = phi->getNumber();
    T val = 0.0;

    for (i = 0; i < numberOfBasFcts; i++)
      val += (*this)[dof_indices[i]]*(*phi->getPhi(i))(lambda);

    return val;
  }

  template<typename T>
  void DOFVector<T>::interpol(AbstractFunction<T, WorldVector<double> > *fct)
  {
    FUNCNAME("interpol");

    TEST_EXIT(interFct = fct)("no function to interpolate\n");

    if (!this->getFESpace())
      {
	MSG("no dof admin in vec %s, skipping interpolation\n",
	    this->getName().c_str());
	return;
      }

    if (!(this->getFESpace()->getAdmin()))
      {
	MSG("no dof admin in feSpace %s, skipping interpolation\n", 
	    this->getFESpace()->getName().c_str());
	return;
      }
  
    if (!(this->getFESpace()->getBasisFcts()))
      {
	MSG("no basis functions in admin of vec %s, skipping interpolation\n",
	    this->getName().c_str());
	return;
      }

    if (!(fct))
      {
	MSG("function that should be interpolated only pointer to NULL, ");
	Msg::print("skipping interpolation\n");
	return;
      }

    traverseVector = this;
    this->getFESpace()->getMesh()->traverse(-1, 
					    Mesh::CALL_LEAF_EL |
					    Mesh::FILL_COORDS, 
					    interpolFct);

    return;
  }

  template<typename T>
  int DOFVector<T>::interpolFct(ElInfo* elinfo)
  {
    int i;
    const BasisFunction *basFct = traverseVector->getFESpace()->getBasisFcts();
    const DOFAdmin* admin = traverseVector->getFESpace()->getAdmin();
    const DegreeOfFreedom   *dof = 
      basFct->getLocalIndices(const_cast<Element *>(elinfo->getElement()), 
			      admin, NULL);
    const T *inter_val = 
      const_cast<BasisFunction*>(basFct)->interpol(elinfo, 
						   0, 
						   NULL,
						   traverseVector->interFct, 
						   NULL);

    int number = basFct->getNumber();
    for (i = 0; i < number; i++)
      (*traverseVector)[dof[i]] = inter_val[i];

    return 0;
  }

  // das hat Andreas eingefuegt: Integral...

  template<typename T>
  double DOFVector<T>::Int(Quadrature* q) const
  {
    FUNCNAME("DOFVector::Int");
  
    Mesh* mesh = feSpace->getMesh();

    int deg;
    dim = mesh->getDim();

    if (!q)
      {
	deg = 2*feSpace->getBasisFcts()->getDegree();
	q = Quadrature::provideQuadrature(dim, deg);
      }

    quad_fast = FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(), *q, INIT_PHI);

    norm = 0.0;
    traverseVector = const_cast<DOFVector<T>*>(this);

    mesh->traverse(-1, 
		   Mesh::CALL_LEAF_EL|
		   Mesh::FILL_COORDS |
		   Mesh::FILL_DET, 
		   Int_fct);

    return norm;  
  }

  template<typename T>
  int DOFVector<T>::Int_fct(ElInfo *elinfo)
  {
    double det, normT;
    const T *uh_vec;
    int iq;

    det = elinfo->getDet();
    //uh_loc = traverseVector->getLocalVector(elinfo->getElement(), NULL);
    //uh_vec = quad_fast->uhAtQp(uh_loc, NULL);

    uh_vec = traverseVector->getVecAtQPs(elinfo,
					 NULL,
					 quad_fast,
					 NULL);

    int numPoints = quad_fast->getNumPoints();
    for (normT = iq = 0; iq < numPoints; iq++)
      {
	normT += quad_fast->getWeight(iq)*(uh_vec[iq]);
      }
    norm += det*normT;

    return 0;
  }

  // ... und die L1-Norm ...

  template<typename T>
  double DOFVector<T>::L1Norm(Quadrature* q) const
  {
    FUNCNAME("DOFVector::L1Norm");
  
    Mesh* mesh = feSpace->getMesh();

    int deg;
    dim = mesh->getDim();

    if (!q)
      {
	deg = 2*feSpace->getBasisFcts()->getDegree();
	q = Quadrature::provideQuadrature(dim, deg);
      }

    quad_fast = FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(),*q,INIT_PHI);

    norm = 0.0;
    traverseVector = const_cast<DOFVector<T>*>(this);

    mesh->traverse(-1, 
		   Mesh::CALL_LEAF_EL|
		   Mesh::FILL_COORDS |
		   Mesh::FILL_DET, 
		   L1Norm_fct);

    return norm;  
  }

  template<typename T>
  int DOFVector<T>::L1Norm_fct(ElInfo *elinfo)
  {
    double det, normT;
    const T *uh_loc, *uh_vec;
    int iq;

    det = elinfo->getDet();
    //   uh_loc = 
    //     traverseVector->getLocalVector(elinfo->getElement(), 
    // 				   NULL);
    //   uh_vec = quad_fast->uhAtQp(uh_loc, NULL);

    uh_vec = traverseVector->getVecAtQPs(elinfo,
					 NULL,
					 quad_fast,
					 NULL);

    int numPoints = quad_fast->getNumPoints();
    for (normT = iq = 0; iq < numPoints; iq++)
      {
	normT += quad_fast->getWeight(iq)*abs(uh_vec[iq]);
      }
    norm += det*normT;

    return 0;
  }


  // bis hierhin gehen Andreas Ergaenzungen...

  template<typename T>
  double DOFVector<T>::L2NormSquare(Quadrature* q) const
  {
    FUNCNAME("DOFVector::L2NormSquare");
  
    Mesh* mesh = feSpace->getMesh();

    int deg;
    dim = mesh->getDim();

    if (!q)
      {
	deg = 2*feSpace->getBasisFcts()->getDegree();
	q = Quadrature::provideQuadrature(dim, deg);
      }

    quad_fast = FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(), 
						      *q, 
						      INIT_PHI);

    norm = 0.0;
    traverseVector = const_cast<DOFVector<T>*>(this);

    mesh->traverse(-1, Mesh::CALL_LEAF_EL|Mesh::FILL_COORDS|Mesh::FILL_DET, L2Norm_fct);

    return norm;  
  }

  template<typename T>
  int DOFVector<T>::L2Norm_fct(ElInfo *elinfo)
  {
    double det, normT;
    const T *uh_vec;
    int iq;

    det = elinfo->getDet();
    //   uh_loc = 
    //     traverseVector->getLocalVector(elinfo->getElement(), 
    // 				   NULL);
    //   uh_vec = quad_fast->uhAtQp(uh_loc, NULL);

    uh_vec = traverseVector->getVecAtQPs(elinfo,
					 NULL,
					 quad_fast,
					 NULL);

    int numPoints = quad_fast->getNumPoints();
    for (normT = iq = 0; iq < numPoints; iq++)
      {
	normT += quad_fast->getWeight(iq)*sqr(uh_vec[iq]);
      }
    norm += det*normT;

    return 0;
  }

  template<typename T>
  int DOFVector<T>::H1Norm_fct(ElInfo *elinfo)
  {
    double norm2, normT;
    //const T *uh_loc;
    const WorldVector<T> *grduh_vec;
    int iq, j;

    double det = elinfo->getDet();
    //const DimVec<WorldVector<double> > &Lambda = elinfo->getGrdLambda();
  
    //det = elinfo->calcGrdLambda(Lambda);
    //uh_loc = 
    //  traverseVector->getLocalVector(elinfo->getElement(), NULL);
    //grduh_vec = 
    //  quad_fast->grdUhAtQp(const_cast<const DimVec<WorldVector<T> >&>(Lambda), uh_loc, NULL);

    grduh_vec = traverseVector->getGrdAtQPs(elinfo, NULL, quad_fast, NULL);

    int dimOfWorld = Global::getGeo(WORLD);

    int numPoints = quad_fast->getNumPoints();
    for (normT = iq = 0; iq < numPoints; iq++)
      {
	for (norm2 = j = 0; j < dimOfWorld; j++)
	  norm2 += sqr(grduh_vec[iq][j]);

	normT += quad_fast->getWeight(iq)*norm2;
      }
    norm += det*normT;

    return 0;
  }


  template<typename T>
  double DOFVector<T>::H1NormSquare(Quadrature *q) const
  {
    FUNCNAME("DOFVector::H1NormSquare");
    int deg;
  
    Mesh *mesh = feSpace->getMesh();
    dim = mesh->getDim();

    if (!q)
      {
	deg = 2*feSpace->getBasisFcts()->getDegree()-2;
	q = Quadrature::provideQuadrature(dim, deg);
      }
    quad_fast = 
      FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(), 
					    *q, 
					    INIT_GRD_PHI);

    norm = 0.0;
    traverseVector = const_cast<DOFVector<T>*>(this);

    mesh->traverse(-1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | 
		   Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA,
		   H1Norm_fct);

    return norm;
  }

  template<typename T>
  void DOFVector<T>::compressDOFIndexed(int first, int last, 
					::std::vector<DegreeOfFreedom> &newDOF)
  {
    int i, j;
    for(i = first; i <= last; i++) {
      if((j = newDOF[i]) >= 0) {
	vec[j] = vec[i];
      }
    }
  }

  template<typename T>
  ElementVector* DOFVectorBase<T>::assemble(T factor, ElInfo *elInfo,
					    const BoundaryType *bound, 
					    Operator *op)
  {
    FUNCNAME("DOFVector::assemble");

    //TEST_EXIT(op || operators.size())("no rhs operator\n");
    if(!(op || this->operators.size())) return NULL;

    Operator *operat = op ? op : this->operators[0];

    //   const DegreeOfFreedom *rowDOF;
    //   const Element         *el = elInfo->getElement();
  
    //   const BasisFunction *psi = operat->getRowFESpace()->getBasisFcts();

    //   DOFAdmin *rowAdmin = operat->getRowFESpace()->getAdmin();

    //   int nRow = psi->getNumber();

    //  if(!elementVector) elementVector = operat->createElementVector();

    //  operat->resetElementVector(elementVector);

    this->elementVector = 
      operat->getAssembler()->initElementVector(this->elementVector, elInfo);

    if(op) {
      op->getElementVector(elInfo, this->elementVector);
    } else {
      ::std::vector<Operator*>::iterator it;
      ::std::vector<double*>::iterator factorIt;
      for(it = this->operators.begin(), factorIt = this->operatorFactor.begin();	
	  it != this->operators.end(); 
	  ++it, ++factorIt) 
	{
	  (*it)->getElementVector(elInfo, 
				  this->elementVector, 
				  *factorIt ? **factorIt : 1.0);
	}
    }

    //   rowDOF = psi->getLocalIndices(el, rowAdmin, NULL);

    addElementVector(factor,
		     *this->elementVector, 
		     bound);

    return this->elementVector;
  }

  template<typename T>
  Flag DOFVectorBase<T>::getAssembleFlag()
  {
    Flag fillFlag(0);
    ::std::vector<Operator*>::iterator op;
    for(op = this->operators.begin(); op != this->operators.end(); ++op) {
      fillFlag |= (*op)->getFillFlag();
    }
    return fillFlag;
  }

  template<typename T>
  DOFVector<T>& DOFVector<T>::operator=(const DOFVector<T>& rhs )
  {
    feSpace=rhs.feSpace;
    vec=rhs.vec;
    this->elementVector=NULL;
    interFct=rhs.interFct;
    refineInter=rhs.refineInter;
    coarsenOperation=rhs.coarsenOperation;
    this->operators=rhs.operators;
    this->operatorFactor=rhs.operatorFactor;
    if (rhs.boundaryManager) {
      if(this->boundaryManager) delete this->boundaryManager; 
      this->boundaryManager = new BoundaryManager(*rhs.boundaryManager);
      //    boundaryManager->setDOFVector(this);
    }
    else 
      this->boundaryManager=NULL;
    return *this;
  }

  template<typename T>
  const DOFVector<T>& operator*=(DOFVector<T>& x, T scal)
  {
    FUNCNAME("operator*=(DOFVector<T>& x, T scal)");
    const DOFAdmin *admin = NULL;

    TEST_EXIT( x.getFESpace() && (admin = x.getFESpace()->getAdmin()))
      ("pointer is NULL: %8X, %8X\n", x.getFESpace(), admin);

    typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&x), USED_DOFS);
    for(vecIterator.reset(); !vecIterator.end(); ++vecIterator) {
      (*vecIterator) *= scal; 
    };
    return x;
  }


  template<typename T>
  const DOFVector<T>& operator+=(DOFVector<T>& x, const DOFVector<T>& y)
  {
    FUNCNAME("operator*=(DOFVector<T>& x, const DOFVector<T>& y)");
    const DOFAdmin *admin;

    TEST_EXIT(x.getFESpace() && y.getFESpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
    TEST_EXIT((admin = x.getFESpace()->getAdmin()) &&
	      (admin == y.getFESpace()->getAdmin()))
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
    TEST_EXIT(x.getSize() == y.getSize())("different sizes\n");

    typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(&x), USED_DOFS);
    typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&y)), USED_DOFS);
    for(xIterator.reset(), yIterator.reset();
	!xIterator.end();
	++xIterator, ++yIterator) 
      {  
	*xIterator += *yIterator; 
      }
    return x;
  }

  // F.H. 17/3/2006
  template<typename T>
  const DOFVector<T>& operator-=(DOFVector<T>& x, const DOFVector<T>& y)
  {
    FUNCNAME("operator*=(DOFVector<T>& x, const DOFVector<T>& y)");
    const DOFAdmin *admin;

    TEST_EXIT(x.getFESpace() && y.getFESpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
    TEST_EXIT((admin = x.getFESpace()->getAdmin()) &&
	      (admin == y.getFESpace()->getAdmin()))
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
    TEST_EXIT(x.getSize() == y.getSize())("different sizes\n");

    typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(&x), USED_DOFS);
    typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&y)), USED_DOFS);
    for(xIterator.reset(), yIterator.reset();
	!xIterator.end();
	++xIterator, ++yIterator) 
      {  
	*xIterator -= *yIterator; 
      }
    return x;
  }




  template<typename T>
  const DOFVector<T>& operator*=(DOFVector<T>& x, const DOFVector<T>& y)
  {
    FUNCNAME("operator*=(DOFVector<T>& x, const DOFVector<T>& y)");
    const DOFAdmin *admin;

    TEST_EXIT(x.getFESpace() && y.getFESpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
    TEST_EXIT((admin = x.getFESpace()->getAdmin()) &&
	      (admin == y.getFESpace()->getAdmin()))
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
    TEST_EXIT(x.getSize() == y.getSize())("different sizes\n");

    typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(&x), USED_DOFS);
    typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&y)), USED_DOFS);
    for(xIterator.reset(), yIterator.reset();
	!xIterator.end();
	++xIterator, ++yIterator) 
      {  
	*xIterator *= *yIterator; 
      }

    return x;
  }
  // F.H. 17/3/2006


  template<typename T>
  T operator*(DOFVector<T>& x, DOFVector<T>& y)
  {
    FUNCNAME("DOFVector<T>::operator*(DOFVector<T>& x, DOFVector<T>& y)");
    T dot;
    const DOFAdmin *admin = NULL;

    TEST_EXIT(x.getFESpace() && y.getFESpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
    TEST_EXIT((admin = x.getFESpace()->getAdmin()) && (admin == y.getFESpace()->getAdmin()))
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
    TEST_EXIT(x.getSize() == y.getSize())("different sizes\n");

    dot = 0;

    typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(&x), USED_DOFS);
    typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(&y), USED_DOFS);
    for (xIterator.reset(), yIterator.reset();
	 !xIterator.end();
	 ++xIterator, ++yIterator) 
      {
	dot += (*xIterator) * (*yIterator);
      };

    return(dot);
  }

  template<typename T>
  void mv(MatrixTranspose transpose,
	  const DOFMatrix &a, 
	  const DOFVector<T>&x,
	  DOFVector<T> &result,
	  bool add)
  {
    FUNCNAME("DOFVector<T>::mv()");

    TEST_EXIT(a.getRowFESpace() && a.getColFESpace() && x.getFESpace() && result.getFESpace())
      ("getFESpace() is NULL: %8X, %8X, %8X, %8X\n", 
       a.getRowFESpace(), a.getColFESpace(), x.getFESpace(), result.getFESpace());

    const DOFAdmin *rowAdmin = a.getRowFESpace()->getAdmin();
    const DOFAdmin *colAdmin = a.getColFESpace()->getAdmin();

    TEST_EXIT((rowAdmin && colAdmin) && 
	      (((transpose == NoTranspose) && (colAdmin == x.getFESpace()->getAdmin()) && 
		(rowAdmin == result.getFESpace()->getAdmin()))||
	       ((transpose == Transpose) && (rowAdmin == x.getFESpace()->getAdmin()) && 
		(colAdmin == result.getFESpace()->getAdmin()))))
      ("no admin or different admins: %8X, %8X, %8X, %8X\n",
       a.getRowFESpace()->getAdmin(), a.getColFESpace()->getAdmin(), 
       x.getFESpace()->getAdmin(), result.getFESpace()->getAdmin());

    if (transpose == NoTranspose) {
      TEST_EXIT(static_cast<int>(x.getSize()) >= colAdmin->getUsedSize())
	("x.size = %d too small: admin->sizeUsed = %d\n",
	 x.getSize(), colAdmin->getUsedSize());
      TEST_EXIT(static_cast<int>( result.getSize()) >= rowAdmin->getUsedSize())
	("size = %d too small: admin->sizeUsed = %d\n",
	 result.getSize(), rowAdmin->getUsedSize());
      TEST_EXIT(static_cast<int>( a.getSize()) >= rowAdmin->getUsedSize())
	("a.size = %d too small: admin->sizeUsed = %d\n",
	 a.getSize(), rowAdmin->getUsedSize());

      // This is the old implementation of the mv-multiplication. It have been changed
      // because of the OpenMP-parallelization:
      //       typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS); 
      //       DOFMatrix::Iterator rowIterator(const_cast<DOFMatrix*>(&a), USED_DOFS);
      //       for(vecIterator.reset(), rowIterator.reset();
      // 	  !rowIterator.end();
      // 	  ++rowIterator, ++vecIterator) { 
      // 	sum = 0;
      // 	if(!add) *vecIterator = 0.0;
      // 	for(::std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
      // 	    colIterator != rowIterator->end();
      // 	    colIterator++) {
      // 	  jcol = colIterator->col;
      // 	  if (jcol >= 0) { // entry used? 
      // 	    sum += (static_cast<double>(colIterator->entry)) * x[jcol];
      // 	  } else {
      // 	    if (jcol == DOFMatrix::NO_MORE_ENTRIES)
      // 	      break;
      // 	  }
      // 	}
      // 	*vecIterator += sum;
      //       }

      int i;
      int maxI = result.getSize();

#ifdef _OPENMP
#pragma omp parallel for schedule(dynamic, 25000) default(shared) private(i)
#endif
      for (i = 0; i < maxI; i++) {
	if (rowAdmin->isDOFFree(i)) {
	  continue;
	}

  	T sum = 0;

 	const ::std::vector<MatEntry> *v = &a[i];

  	for (int j = 0; j < static_cast<int>(v->size()); j++) {
	  const MatEntry *m = &((*v)[j]);
	  int jcol = m->col;

  	  if (jcol >= 0) { // entry used? 
  	    sum += (static_cast<double>(m->entry)) * x[jcol];
  	  } else {
  	    if (jcol == DOFMatrix::NO_MORE_ENTRIES)
  	      break;
  	  }
  	}

  	if (add) {
  	  result[i] += sum;
  	} else {
  	  result[i] = sum;
  	}
      }

    } else if (transpose == Transpose) {
      TEST_EXIT(static_cast<int>(x.getSize()) >= rowAdmin->getUsedSize())
	("x.size = %d too small: admin->sizeUsed = %d\n",
	 x.getSize(), rowAdmin->getUsedSize());
      TEST_EXIT(static_cast<int>( result.getSize()) >= colAdmin->getUsedSize())
	("size = %d too small: admin->sizeUsed = %d\n",
	 result.getSize(), colAdmin->getUsedSize());
      TEST_EXIT(static_cast<int>( a.getSize()) >= rowAdmin->getUsedSize())
	("a.size = %d too small: admin->sizeUsed = %d\n",
	 a.getSize(), rowAdmin->getUsedSize());

      if(!add) {
	typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS);
	for(vecIterator.reset();
	    !vecIterator.end(); ++vecIterator) {
	  *vecIterator = 0; 
	};
      }

      typename DOFVector<T>::Iterator xIterator(const_cast<DOFIndexed<T>*>(dynamic_cast<const DOFIndexed<T>*>(&x)), USED_DOFS);
      DOFMatrix::Iterator rowIterator(const_cast<DOFMatrix*>(&a), USED_DOFS);
      for(xIterator.reset(), rowIterator.reset();
	  !rowIterator.end();
	  ++rowIterator, ++xIterator) { 
	T ax = (*xIterator);
	for(::std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
	    colIterator != rowIterator->end();
	    colIterator++) {
	  int jcol = colIterator->col;
	  if (jcol >= 0) // entry used?
	    result[jcol] += (static_cast<double>(colIterator->entry)) * ax;
	  else 
	    if (jcol == DOFMatrix::NO_MORE_ENTRIES) break;
	}
      }
    } else {
      ERROR_EXIT("transpose = %d\n", transpose);
    }
  }

  template<typename T>
  void axpy(double alpha, 
	    const DOFVector<T>& x, 
	    DOFVector<T>& y)
  {
    FUNCNAME("DOFVector<T>::axpy()");

    TEST_EXIT(x.getFESpace() && y.getFESpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());

    const DOFAdmin *admin = x.getFESpace()->getAdmin();

    TEST_EXIT((admin) && (admin == y.getFESpace()->getAdmin()))
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
    TEST_EXIT(static_cast<int>(x.getSize()) >= admin->getUsedSize())
      ("size = %d too small: admin->size = %d\n", x.getSize(), 
       admin->getUsedSize());
    TEST_EXIT(static_cast<int>(y.getSize()) >= admin->getUsedSize())
      ("y.size = %d too small: admin->size = %d\n", y.getSize(), 
       admin->getUsedSize());

      // This is the old implementation of the mv-multiplication. It have been changed
      // because of the OpenMP-parallelization:
//     typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&x)), USED_DOFS);
//     typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(&y), USED_DOFS);
//     for(xIterator.reset(), yIterator.reset();
// 	!xIterator.end();
// 	++xIterator, ++yIterator) 
//       {  
// 	*yIterator += alpha * (*xIterator); 
//       };

      int i;
      int maxI = y.getSize();

#ifdef _OPENMP
#pragma omp parallel for schedule(dynamic, 25000) default(shared) private(i)
#endif
      for (i = 0; i < maxI; i++) {
	if (!admin->isDOFFree(i)) {
	  y[i] += alpha * x[i];
	}
      }
  }

  template<typename T>
  const DOFVector<T>& operator*(const DOFVector<T>& v, double d)
  {
    static DOFVector<T> result;
    return mult(d, v, result);
  }

  template<typename T>
  const DOFVector<T>& operator*(double d, const DOFVector<T>& v)
  {
    static DOFVector<T> result;
    return mult(d, v, result);
  }

  template<typename T>
  const DOFVector<T>& operator+(const DOFVector<T> &v1 , const DOFVector<T> &v2)
  {
    static DOFVector<T> result;
    return add(v1, v2, result);
  }


  template<typename T>
  void xpay(double alpha, 
	    const DOFVector<T>& x, 
	    DOFVector<T>& y)
  {
    FUNCNAME("DOFVector<T>::xpay()");

    TEST_EXIT(x.getFESpace() && y.getFESpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());

    const DOFAdmin *admin = x.getFESpace()->getAdmin();

    TEST_EXIT(admin && (admin == y.getFESpace()->getAdmin()))
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
    TEST_EXIT(static_cast<int>(x.getSize()) >= admin->getUsedSize())
      ("size = %d too small: admin->size = %d\n", x.getSize(), 
       admin->getUsedSize());
    TEST_EXIT(static_cast<int>(y.getSize()) >= admin->getUsedSize())
      ("y.size = %d too small: admin->size = %d\n", y.getSize(), 
       admin->getUsedSize());

    // This is the old implementation of the mv-multiplication. It have been changed
    // because of the OpenMP-parallelization:
    //     typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&x)), USED_DOFS);
    //     typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(&y), USED_DOFS);
    //     for(xIterator.reset(), yIterator.reset();
    // 	!xIterator.end();
    // 	++xIterator, ++yIterator) 
    //       {  
    // 	*yIterator = alpha *(*yIterator)+ (*xIterator); 
    //       };

      int i;
      int maxI = y.getSize();

#ifdef _OPENMP
#pragma omp parallel for schedule(dynamic, 25000) default(shared) private(i)
#endif
      for (i = 0; i < maxI; i++) {
	if (!admin->isDOFFree(i)) {
	  y[i] = alpha * y[i] + x[i];
	}
      }
  }

  template<typename T>
  inline const DOFVector<T>& mult(double scal,
				  const DOFVector<T>& v,
				  DOFVector<T>& result)
  {
    typename DOFVector<T>::Iterator vIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v)), USED_DOFS);
    typename DOFVector<T>::Iterator rIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS);
    for(vIterator.reset(), rIterator.reset();
	!vIterator.end();
	++vIterator, ++rIterator) 
      {  
	*rIterator = scal * (*vIterator);
      };

    return result;
  }

  template<typename T>
  inline const DOFVector<T>& add(const DOFVector<T>& v,
				 double scal,
				 DOFVector<T>& result)
  {
    typename DOFVector<T>::Iterator vIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v)), USED_DOFS);
    typename DOFVector<T>::Iterator rIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS);
    for(vIterator.reset(), rIterator.reset();
	!vIterator.end();
	++vIterator, ++rIterator) 
      {  
	*rIterator = (*vIterator) + scal;
      };
    return result;
  }

  template<typename T>
  inline const DOFVector<T>& add(const DOFVector<T>& v1,
				 const DOFVector<T>& v2,
				 DOFVector<T>& result)
  {
    typename DOFVector<T>::Iterator v1Iterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v1)), USED_DOFS);
    typename DOFVector<T>::Iterator v2Iterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v2)), USED_DOFS);
    typename DOFVector<T>::Iterator rIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS);
    for(v1Iterator.reset(), v2Iterator.reset(), rIterator.reset();
	!v1Iterator.end();
	++v1Iterator, ++v2Iterator, ++rIterator) 
      {  
	*rIterator = (*v1Iterator) + (*v2Iterator);
      };
    return result;

  }

  template<typename T>
  const T *DOFVectorBase<T>::getLocalVector(const Element *el, T *d) const
  {
    static T* localVec = NULL;
    const DOFAdmin* admin = feSpace->getAdmin();

    int i;
    int nBasFcts = feSpace->getBasisFcts()->getNumber();

    T *result;
  
    if(d) {
      result = d;
    } else {
      if(localVec) delete [] localVec;
      localVec = new T[nBasFcts]; 
      result = localVec;
    }

    const DegreeOfFreedom *localIndices = 
      feSpace->getBasisFcts()->getLocalIndices(el, admin, NULL);

    for(i = 0; i < nBasFcts; i++) {
      result[i] = (*this)[localIndices[i]];
    }

    return result;
  }

  template<typename T>
  const T *DOFVectorBase<T>::getVecAtQPs(const ElInfo         *elInfo, 
					 const Quadrature     *quad,
					 const FastQuadrature *quadFast,
					 T                    *vecAtQPs) const
  {
    FUNCNAME("DOFVector<T>::getVecAtQPs()");
  
    TEST_EXIT(quad || quadFast)("neither quad nor quadFast defined\n");

    if(quad && quadFast) {
      TEST_EXIT(quad == quadFast->getQuadrature())
	("quad != quadFast->quadrature\n");
    }

    TEST_EXIT(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
      ("invalid basis functions");

    Element *el = elInfo->getElement();

    const Quadrature *quadrature = quadFast ? quadFast->getQuadrature() : quad;

    const BasisFunction *basFcts = feSpace->getBasisFcts();

    int numPoints = quadrature->getNumPoints();
    int nBasFcts  = basFcts->getNumber();
    int i, j;

    static T *localvec = NULL;

    T *result;

    if (vecAtQPs) {
      result = vecAtQPs;
    } else {
      if(localvec) delete [] localvec;
      localvec = new T[numPoints];
      for(i = 0; i < numPoints; i++) {
	localvec[i] = 0.0;
      }
      result = localvec;
    }
  
    const T *localVec = getLocalVector(el, NULL);

    for(i = 0; i < numPoints; i++) {
      for(result[i] = j = 0; j < nBasFcts; j++) {
	result[i] +=  localVec[j] * 
	  (quadFast ? 
	   (quadFast->getPhi(i, j)) : 
	   ((*(basFcts->getPhi(j)))(quad->getLambda(i))));
      }
    }

    return const_cast<const T*>(result);
  }

}