DOFIndexed.cc

//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


#include "DOFIndexed.h"
#include "DOFMatrix.h"

#include <boost/numeric/mtl/mtl.hpp>
#include <boost/numeric/mtl/utility/tag.hpp>
#include <boost/numeric/mtl/utility/category.hpp>
#include <boost/numeric/linear_algebra/identity.hpp>


// Defining the interface for MTL4
namespace mtl {

  // Let MTL4 know that DOFIndexed it is a column vector
  namespace traits {
    template <typename T>
    struct category< AMDiS::DOFIndexed<T> > 
    {
      typedef tag::dense_col_vector type;
    };
  }

  namespace ashape {
    template <typename T>
    struct ashape< AMDiS::DOFIndexed<T> > 
    {
      typedef cvec<typename ashape<T>::type> type;
    };
  }

  // Modelling Collection and MutableCollection
  template <typename T>
  struct Collection< AMDiS::DOFIndexed<T> > 
  {
    typedef T               value_type;
    typedef const T&        const_reference;
    typedef std::size_t     size_type;
  };

  template <typename T>
  struct MutableCollection< AMDiS::DOFIndexed<T> > 
    : public Collection< AMDiS::DOFIndexed<T> > 
  {
    typedef T&              reference;
  };


} // namespace mtl



namespace AMDiS {

  // Some free functions used in MTL4

  template <typename T>
  inline std::size_t size(const AMDiS::DOFIndexed<T>& v)
  {
    return v.getSize();
  }

  template <typename T>
  inline std::size_t num_rows(const AMDiS::DOFIndexed<T>& v)
  {
    return v.getSize();
  }

  template <typename T>
  inline std::size_t num_cols(const AMDiS::DOFIndexed<T>& v)
  {
    return 1;
  }
 

  template <typename T>
  inline void set_to_zero(AMDiS::DOFIndexed<T>& v)
  {
    using math::zero;
    T  ref, my_zero(zero(ref));

    std::fill(v.begin(), v.end(), my_zero);
  }


  void mv(MatrixTranspose transpose, 
	  const DOFMatrix &a, 
	  const DOFIndexed<double> &x,
	  DOFIndexed<double> &result,
	  bool add)
  {
    FUNCNAME("DOFIndexed<T>::mv");
    TEST_EXIT(false)("This function is not supported any more\n");
#if 0
    if (transpose == NoTranspose)
      mult(a.getBaseMatrix(), x, result);
    else if (transpose == Transpose)
      mult(trans(const_cast<DOFMatrix::base_matrix_type&>(a.getBaseMatrix())), x, result);
    else 
      ERROR_EXIT("transpose = %d\n", transpose);

    int irow, jcol;
    double sum;

    if (transpose == NoTranspose) {
      DOFMatrix::Iterator rowIterator(const_cast<DOFMatrix*>(&a), USED_DOFS);
      for(rowIterator.reset(); !rowIterator.end(); ++rowIterator) { 
	sum = 0;
	irow = rowIterator.getDOFIndex();
	if(!add) result[irow] = 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;
	    }
	  }
	result[irow] += sum;
      }
    } else {
      ERROR_EXIT("transpose=%d\n", transpose);
    }  
#endif
  }
}