// // 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 } }