MatrixVector corrected

AMDiS/CMakeLists.txt

@@ -685,6 +685,11 @@ INSTALL(FILES ${HEADERS}
 	DESTINATION include/amdis/traits/)
 list(APPEND deb_add_dirs "include/amdis/traits")
 
+FILE(GLOB HEADERS "${SOURCE_DIR}/utility/*.h*")
+INSTALL(FILES ${HEADERS} 
+	DESTINATION include/amdis/utility/)
+list(APPEND deb_add_dirs "include/amdis/utility")
+
 FILE(GLOB HEADERS "${SOURCE_DIR}/time/*.h")
 INSTALL(FILES ${HEADERS} 
 	DESTINATION include/amdis/time/)

AMDiS/src/DOFVector.cc

@@ -486,7 +486,7 @@ namespace AMDiS {
 	    elInfo2->worldToCoord(worldVec, &coords2);
 	  
 	    bool isPositive = true;
-	    for (int j = 0; j < coords2.size(); j++) {
+	    for (int j = 0; j < coords2.getSize(); j++) {
 	      if (coords2[j] < -0.00001) {
 		isPositive = false;
 		break;

AMDiS/src/ElInfo.h

@@ -215,8 +215,8 @@ namespace AMDiS {
     
     inline void getGrdLambda(mtl::dense2D<double>& grd_lam)
     {
-      grd_lam.change_dim(grdLambda.size(), Global::getGeo(WORLD));
-      for (size_t i = 0; i < static_cast<size_t>(grdLambda.size()); i++)
+      grd_lam.change_dim(grdLambda.getSize(), Global::getGeo(WORLD));
+      for (size_t i = 0; i < static_cast<size_t>(grdLambda.getSize()); i++)
 	for (size_t j = 0; j < static_cast<size_t>(Global::getGeo(WORLD)); j++)
 	  grd_lam(i,j) = grd_lam[i][j];
     }

AMDiS/src/FixVec.h

@@ -53,6 +53,8 @@ namespace AMDiS {
   template<typename T,GeoIndex d>
   class FixVec : public Vector<T>
   {
+    typedef FixVec     self;
+    typedef Vector<T>  super;
   public:
     
     /// Constructor without initialisation. initType must be NO_INIT. If dim is
@@ -80,6 +82,10 @@ namespace AMDiS {
       TEST_EXIT_DBG(initType == DEFAULT_VALUE)("wrong initType or wrong initializer\n");
       this->set(ini);
     }
+    
+    FixVec(self const& other)
+      : super(other)
+    { }
 
     /// Initialisation for dim.
     inline void init(int dim) 
@@ -94,10 +100,10 @@ namespace AMDiS {
     }
 
     /// Returns the \ref size_ of the FixVec.
-    inline int size() const 
-    { 
-      return this->getSize(); 
-    } 
+//     inline int size() const 
+//     { 
+//       return this->getSize(); 
+//     } 
 
   protected:
     /// Determines needed vector size.
@@ -178,7 +184,7 @@ namespace AMDiS {
     }
 
     /// Destructor
-    virtual ~VectorOfFixVecs()
+    ~VectorOfFixVecs()
     {
       for (int i = 0; i < size; i++)
 	delete vec[i];
@@ -274,7 +280,7 @@ namespace AMDiS {
     }
 
     /// destructor
-    virtual ~MatrixOfFixVecs()
+    ~MatrixOfFixVecs()
     {
       for (VectorOfFixVecs<FixVecType>** i = &vec[0]; i < &vec[rows]; i++)
 	delete *i;
@@ -417,18 +423,19 @@ namespace AMDiS {
     {}
     
     /// Copy assignement operator
-    self& operator=(self other) 
-    {
-      swap(*this, other);
-      return *this;
-    }
+//     self& operator=(self const& other)
+//     {
+//       assert( Global::getGeo(WORLD) == other.getSize() );
+//       this->setValues(other.getValArray());
+//       return *this;
+//     }
       
     
     /// Assignement operator
     template <typename S>
     self& operator=(const Vector<S>& other) 
     {
-      TEST_EXIT_DBG( other.getSize() == Global::getGeo(WORLD) )
+      TEST_EXIT_DBG( other.getSize() == this->size )
 		   ("Wrong dimensions in assignment.\n");
       this->setValues(other.getValArray());
       return *this;
@@ -503,26 +510,27 @@ namespace AMDiS {
     { }
     
     /// Copy assignment operator
-    self& operator=(self other) 
-    {
-      swap(*this, other);
-      return *this;
-    }
+//     self& operator=(self other) 
+//     {
+//       swap(*this, other);
+//       return *this;
+//     }
     
     /// Assignment operator
     template <typename S>
     self& operator=(const Matrix<S>& other) 
     {
-      TEST_EXIT_DBG( other.getNumOfRows() == Global::getGeo(WORLD) &&
-		     other.getNumOfCols() == Global::getGeo(WORLD) )
-		    ("Wrong dimensions in assignment.\n");
+      TEST_EXIT_DBG( this->size == other.getSize() && 
+		     this->rows == other.getNumRows() && 
+		     this->cols == other.getNumCols() )
+		   ("Wrong dimensions in assignment.\n");
       this->setValues(other.getValArray());
       return *this;
     }
     
     /// Assignment operator for scalars
     template <typename S>
-    typename enable_if<boost::is_convertible<S, T>, WorldMatrix<T> >::type &
+    typename enable_if<boost::is_convertible<S, T>, self>::type &
     operator=(S value) 
     {
       this->set(value);

AMDiS/src/FixVec.hh

@@ -87,8 +87,8 @@ namespace AMDiS {
   {
     FUNCNAME_DBG("WorldMatrix<T>::vecProduct()");
 
-    TEST_EXIT_DBG(v1.getSize() == v2.getSize())("invalid size 1\n");
-    TEST_EXIT_DBG(v1.getSize() == this->getSize())("invalid size 2\n");
+    TEST_EXIT_DBG(v1.getSize() == v2.getSize())("size(v1) != size(v2), %d != %d\n", v1.getSize(), v2.getSize());
+    TEST_EXIT_DBG(v1.getSize() == this->getNumRows())("size(v1) != num_rows(this), %d != %d\n", v1.getSize(), this->getNumRows());
 
     T *thisIt = this->begin();
 

AMDiS/src/Global.cc

@@ -61,7 +61,7 @@ namespace AMDiS {
     FUNCNAME("Msg::wait()");
 
     if (w) {
-      std::string line;
+      char line;
       MSG("wait for <enter> ...");
       std::cin >> line;
 //       char* result = fgets(line, 9, stdin);

AMDiS/src/Global.h

@@ -361,7 +361,7 @@ namespace AMDiS {
   /// function name nn for message output via MSG, WARNING, ...
 #define FUNCNAME(nn) const char *funcName; funcName = nn;
   
-#if (DEBUG == 0) 
+#ifdef NDEBUG
   #define FUNCNAME_DBG(nn)
 #else
   #define FUNCNAME_DBG(nn) const char *funcName; funcName = nn;
@@ -386,18 +386,18 @@ namespace AMDiS {
 #define TEST_EXIT(test) if ((test));else ERROR_EXIT
 
   /// In debug mode, it corresponds to ERROR_EXIT, otherwise it is noop.
-#if (DEBUG == 0) 
+#ifdef NDEBUG
   #define TEST_EXIT_DBG(test) if (false) Msg::catch_error_exit
   #define DBG_VAR(var)
 #else
-#define TEST_EXIT_DBG(test) if ((test));else ERROR_EXIT
+  #define TEST_EXIT_DBG(test) if ((test));else ERROR_EXIT
   #define DBG_VAR(var) var
 #endif
 
   /// prints a message
 #define MSG Msg::print_funcname(funcName), Msg::print
 
-#if (DEBUG == 0) 
+#ifdef NDEBUG
   #define MSG_DBG
 #else
   #define MSG_DBG Msg::print_funcname(funcName), Msg::print

AMDiS/src/MatrixVector.h

@@ -56,31 +56,11 @@ namespace AMDiS {
     {
       setValues(other.valArray);
     }
-    
-    /// Copy constructor for other of different value_type
-    // TODO: notwendig?
-//     template <typename S>
-//     Vector(Vector<S> const& rhs)
-//       : Serializable(),
-// 	size(rhs.getSize()),
-// 	valArray(size ? new T[size] : NULL)
-//     {
-//       setValues(rhs.valArray);
-//     }
-    
-// #if HAS_RVALUE_REFERENCES
-//     /// move constructor
-//     Vector(self&& other)
-//       : Vector()
-//     {
-//       swap(*this, other);
-//     }
-// #endif
 
     /// Destructor.
-    virtual ~Vector() 
+    ~Vector() 
     {
-      if (valArray != NULL) {
+      if (valArray) {
 	delete [] valArray; 
 	valArray = NULL;
       }
@@ -95,7 +75,7 @@ namespace AMDiS {
     inline void resize(int s) 
     {
       if (size != s && s > 0) {
-	if (valArray != NULL) 
+	if (valArray) 
 	  delete [] valArray;
 	valArray = new T[s];
 	size = s;
@@ -103,24 +83,9 @@ namespace AMDiS {
     }
 
     /// Assignement operator
-    template <typename S>
-    typename disable_if< boost::is_same<S, T>, self >::type &
-    operator=(Vector<S> const& other) 
-    {
-      assert( size == other.getSize() );
-//       resize(other.getSize());
-      setValues(other.getValArray());
-      return *this;
-    }
-    
-    /// copy assignment operator
-//     self& operator=(self other) 
-//     {
-//       swap(*this, other);
-
     self& operator=(self const& other)
     {
-      assert( size == other.getSize() );
+      resize( other.getSize() );
       setValues(other.getValArray());
       return *this;
     }
@@ -140,73 +105,33 @@ namespace AMDiS {
       setValues(vec);
       return *this;
     }
-  
-    /// A swap function for Vector, used in the Copy-and-swap idiom
-    // need non-templated arguments in order to eliminate a friend declaration warning in gcc
-    friend void swap(Vector& first, Vector& second)
-    {
-      using std::swap; // enable ADL
-      swap(first.size, second.size);
-      swap(first.valArray, second.valArray);
-    }
 
     /// Sets all entries to scal.
     template <typename S>
     typename enable_if< boost::is_convertible<S, T> >::type
     set(S value) 
     {
-      for (T *thisIt = this->begin(); thisIt != this->end(); ++thisIt)
-	*thisIt = value;
-//       std::fill(begin(), end(), value);
+      std::fill(begin(), end(), value);
     }
 
     /// Sets all entries.
     template <typename S>
     inline void setValues(const S* values) 
     {
-      T *thisIt;
-      const S *valuesIt;
-      for (thisIt = this->begin(), valuesIt = values; 
-	   thisIt != this->end(); 
-	   ++thisIt, ++valuesIt) 
-	*thisIt = *valuesIt;
-//       std::copy(values, values + size, valArray);
-    }
-
-    /// Comparison operator.
-    inline bool operator==(Vector<T> const& rhs) const 
-    {
-      if (size != rhs.size) 
-	return false;
-
-      T const* rhsIt;
-      T const* thisIt;
-      for (rhsIt = rhs.begin(), thisIt = this->begin();
-	   rhsIt != rhs.end();
-	   ++rhsIt, ++thisIt)
-	if (*thisIt != *rhsIt) 
-	  return false;
-
-      return true;
-    }
-
-    /// Comparison operator.
-    inline bool operator!=(Vector<T> const& rhs) const 
-    {
-      return !(*this==rhs);
+      std::copy(values, values + size, begin());
     }
 
     /// Access to the i-th vector element.
     inline T& operator[](int i) 
     {
-//       TEST_EXIT_DBG(i < size && i >= 0)("Invalid index %d!\n", i);
+      TEST_EXIT_DBG((unsigned)i < (unsigned)size)("Invalid index %d!\n", i);
       return valArray[i];
     }
 
     /// Access to the i-th vector element for const vectors.
     inline const T& operator[] (int i) const 
     {
-//       TEST_EXIT_DBG(i < size && i >= 0)("Invalid index %d!\n", i);
+      TEST_EXIT_DBG((unsigned)i < (unsigned)size)("Invalid index %d!\n", i);
       return valArray[i];
     }
 
@@ -288,15 +213,6 @@ namespace AMDiS {
 	rows(other.getNumRows()),
 	cols(other.getNumCols())
     { }
-    
-// #if HAS_RVALUE_REFERENCES
-//     /// move constructor
-//     Matrix(self&& other)
-//       : Matrix()
-//     {
-//       swap(*this, other);
-//     }
-// #endif
 
     /// Changes the size of the matrix to newRows x newCols.
     inline void resize(int newRows, int newCols) 
@@ -307,60 +223,22 @@ namespace AMDiS {
 	cols = newCols;
       }
     }
-    
-//     self& operator=(self other) 
-//     {
-//       swap(*this, other);
-//       return *this;
-//     }
-    
+        
     self& operator=(self const& other)
     {
-      assert( this->size == other.getSize() && rows == other.getNumRows() && cols == other.getNumCols() );
-      this->setValues(other.getValArray());
-      return *this;
-    }
-
-    /// Assignement operator
-    template <typename S>
-    inline self& operator=(const Matrix<S>& other) 
-    {
-      assert( rows == other.getNumRows() && cols == other.getNumCols() );
-//       resize(other.getNumRows(), other.getNumCols());
+      resize( other.getNumRows(), other.getNumCols() );
       this->setValues(other.getValArray());
       return *this;
     }
 
     /// Assignement operator for scalars.
-    self& operator=(T value) 
+    template <typename S>
+    typename enable_if< boost::is_convertible<S, T>, self>::type &
+    operator=(S value) 
     {
-      super::set(value);
+      this->set(value);
       return *this;
     }
-  
-    /// A swap function for Matrix, used in the Copy-and-swap idiom
-    // need non-templated arguments in order to eliminate a friend declaration warning in gcc
-    friend void swap(Matrix& first, Matrix& second)
-    {
-      using std::swap; // enable ADL
-      swap(first.rows, second.rows);
-      swap(first.cols, second.cols);
-      swap(static_cast<super&>(first), static_cast<super&>(second));
-    }
-
-    /// Comparison operator.
-    inline bool operator==(const self& rhs) const 
-    {
-      if (rows != rhs.getNumRows()) return false;
-      if (cols != rhs.getNumCols()) return false;
-      return super::operator==(rhs);
-    }
-
-    /// Comparison operator.
-    inline bool operator!=(const self& rhs) const 
-    {
-      return !(*this == rhs);
-    }
 
     /// Access to i-th matrix row.
     inline T *operator[](int i) 

AMDiS/src/MatrixVectorOperations.h

@@ -473,6 +473,20 @@ namespace AMDiS {
     return m1;
   }
 
+  /// matrix += matrix
+  template <typename T, typename S>
+  Matrix<T>& operator+=(Matrix<T>& m1, const Matrix<S>& m2)
+  {
+    T* m1It;
+    S const* m2It;
+    for (m1It = m1.begin(), m2It = m2.begin();
+	 m1It != m1.end(); 
+	 m1It++, m2It++)
+      *m1It += *m2It;
+
+    return m1;
+  }
+
   /// matrix := matrix + matrix
   template <typename T, typename S>
   WorldMatrix<T> operator+(WorldMatrix<T> M1, const WorldMatrix<S>& M2 )
@@ -480,6 +494,14 @@ namespace AMDiS {
     M1 += M2;
     return M1;
   }
+
+  /// matrix := matrix + matrix
+  template <typename T, typename S>
+  Matrix<T> operator+(Matrix<T> M1, const Matrix<S>& M2 )
+  {
+    M1 += M2;
+    return M1;
+  }
   
   /// matrix -= matrix
   template <typename T, typename S>
@@ -494,6 +516,20 @@ namespace AMDiS {
 
     return m1;
   }
+  
+  /// matrix -= matrix
+  template <typename T, typename S>
+  Matrix<T>& operator-=(Matrix<T>& m1, const Matrix<S>& m2)
+  {
+    T *m1It;
+    S const* m2It;
+    for (m1It = m1.begin(), m2It = m2.begin();
+	 m1It != m1.end(); 
+	 m1It++, m2It++)
+      *m1It -= *m2It;
+
+    return m1;
+  }
 
   /// matrix := matrix - matrix
   template <typename T, typename S>
@@ -502,6 +538,14 @@ namespace AMDiS {
     M1 -= M2;
     return M1;
   }
+
+  /// matrix := matrix - matrix
+  template <typename T, typename S>
+  Matrix<T> operator-(Matrix<T> M1, const Matrix<S>& M2 )
+  {
+    M1 -= M2;
+    return M1;
+  }
   
   // unary minus operators
   // ---------------------
@@ -548,6 +592,50 @@ namespace AMDiS {
     return true;
   }
   
+
+  /// Comparison operator for Vector<T>.
+  template <typename T>
+  inline bool operator==(Vector<T> const& lhs, Vector<T> const& rhs) 
+  {
+    if (lhs.getSize() != rhs.getSize()) 
+      return false;
+
+    T const* rhsIt;
+    T const* lhsIt;
+    for (rhsIt = rhs.begin(), lhsIt = lhs.begin();
+	  rhsIt != rhs.end();
+	  ++rhsIt, ++lhsIt)
+      if (*lhsIt != *rhsIt) 
+	return false;
+
+    return true;
+  }
+
+  /// Comparison operator for Vector<T>.
+  template <typename T>
+  inline bool operator!=(Vector<T> const& lhs, Vector<T> const& rhs) 
+  {
+    return !(lhs == rhs);
+  }
+  
+  
+
+  /// Comparison operator for Matrix<T>.
+  template <typename T>
+  inline bool operator==(Matrix<T> const& lhs, Matrix<T> const& rhs) 
+  {
+    if (lhs.getNumRows() != rhs.getNumRows()) return false;
+    if (lhs.getNumCols() != rhs.getNumCols()) return false;
+    return (static_cast<Vector<T> const&>(lhs) == static_cast<Vector<T> const&>(rhs));
+  }
+
+  /// Comparison operator for Matrix<T>.
+  template <typename T>
+  inline bool operator!=(Matrix<T> const& lhs, Matrix<T> const& rhs) 
+  {
+    return !(lhs == rhs);
+  }
+  
   // special operators
   // -----------------
   

AMDiS/src/Recovery.h

@@ -42,7 +42,7 @@ namespace AMDiS {
 	exponent(expon)
     {
       degree_ = exponent[0];
-      for (int i = 1; i<exponent.size(); i++)
+      for (int i = 1; i<exponent.getSize(); i++)
 	degree_ += exponent[i];
     }
     
@@ -52,7 +52,7 @@ namespace AMDiS {
 		      const WorldVector<double>& z) const
     {
       double result = std::pow(y[0] - z[0], double(exponent[0]));
-      for (int i = 1; i < exponent.size(); i++)
+      for (int i = 1; i < exponent.getSize(); i++)
 	result *= std::pow(y[i] - z[i], double(exponent[i]));
       
       return result;

AMDiS/src/VertexInfo.h

@@ -26,6 +26,7 @@
 #define AMDIS_VERTEXINFO_H
 
 #include "FixVec.h"
+#include "MatrixVectorOperations.h"
 
 namespace AMDiS {
   

AMDiS/src/expressions/value_expr.hpp

@@ -164,11 +164,10 @@ namespace AMDiS
     struct Eye : public LazyOperatorTermBase
     {
       typedef Matrix<int> value_type;
-      value_type M;
       size_t n;
+      value_type M;
       
-      Eye(int n) : M(n,n) {
-	assert((N >= 0 && n < 0) || (N < 0 && n >= 0));
+      Eye(int n_) : n(n_ < 0 ? N : n_), M(n,n) {
 	M = 0;
 	for (size_t i = 0; i < n; ++i)
 	  M[i][i] = 1;

AMDiS/src/io/VtkVectorWriter.cc

@@ -84,18 +84,18 @@ namespace AMDiS
       {
 	DOFVector<WorldVector<double> > *newValues = new DOFVector<WorldVector<double> >(values[0]->getFeSpace(), values[0]->getName());
 	WorldVector<DOFIterator<double>* > iterators;
-	for (size_t i = 0; i < static_cast<size_t>(values.getSize()); i++)
+	for (int i = 0; i < values.getSize(); i++)
 	  iterators[i] = new DOFIterator<double>(values[i],USED_DOFS);
-	for (size_t i = 0; i < static_cast<size_t>(iterators.getSize()); i++)
+	for (int i = 0; i < iterators.getSize(); i++)
 	  iterators[i]->reset();
 	DOFIterator<WorldVector<double> > resultIter(newValues, USED_DOFS);
 
 	for(resultIter.reset(); !resultIter.end(); resultIter++)
 	{
-	  for (size_t i = 0; i < static_cast<size_t>(iterators.getSize()); i++)