- add dot procut to dense vector

- update fvector test
(thanks to Matthias Wohlmuth)

[[Imported from SVN: r6863]]

dune/common/densevector.hh

@@ -8,7 +8,8 @@
 #include "genericiterator.hh"
 #include "ftraits.hh"
 #include "matvectraits.hh"
-
+#include "promotiontraits.hh"
+#include "dotproduct.hh"
 
 namespace Dune {
 
@@ -453,16 +454,39 @@ namespace Dune {
       return asImp();
     }
 
-    //===== Euclidean scalar product
+    /**
+     * \brief indefinite vector dot product \f$\left (x^T \cdot y \right)\f$ which corresponds to Petsc's VecTDot
+     *
+     * http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecTDot.html
+     * @param y other vector
+     * @return
+     */
+    template<class Other>
+    typename PromotionTraits<field_type,typename DenseVector<Other>::field_type>::PromotedType operator* (const DenseVector<Other>& y) const {
+      typedef typename PromotionTraits<field_type, typename DenseVector<Other>::field_type>::PromotedType PromotedType;
+      PromotedType result = PromotedType();
+      assert(y.size() == size());
+      for (size_type i=0; i<size(); i++) {
+        result += PromotedType((*this)[i]*y[i]);
+      }
+      return result;
+    }
 
-    //! scalar product (x^T y)
-    template <class Other>
-    value_type operator* (const DenseVector<Other>& y) const
-    {
+    /**
+     * @brief vector dot product \f$\left (x^H \cdot y \right)\f$ which corresponds to Petsc's VecDot
+     *
+     * http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecDot.html
+     * @param y other vector
+     * @return
+     */
+    template<class Other>
+    typename PromotionTraits<field_type,typename DenseVector<Other>::field_type>::PromotedType dot(const DenseVector<Other>& y) const {
+      typedef typename PromotionTraits<field_type, typename DenseVector<Other>::field_type>::PromotedType PromotedType;
+      PromotedType result = PromotedType();
       assert(y.size() == size());
-      value_type result( 0 );
-      for (size_type i=0; i<size(); i++)
-        result += (*this)[i]*y[i];
+      for (size_type i=0; i<size(); i++) {
+        result += Dune::dot((*this)[i],y[i]);
+      }
       return result;
     }
 

dune/common/test/fvectortest.cc

@@ -70,6 +70,7 @@ struct FieldVectorMainTest
 
     // test scalar product, axpy
     a = v * w;
+    a = v.dot(w);
     z = v.axpy(a,w);
 
     // test comparison
@@ -201,6 +202,76 @@ struct ScalarOrderingTest
   }
 };
 
+// scalar ordering doesn't work for complex numbers
+template <class rt, int d>
+struct DotProductTest
+{
+  DotProductTest() {
+    typedef std::complex<rt> ct;
+    const rt myEps(1e-6);
+
+    dune_static_assert(
+      ( Dune::is_same< typename Dune::FieldTraits<rt>::real_type, rt>::value ),
+      "DotProductTest requires real data type as template parameter!"
+      );
+
+    const ct I(0.,1.); // imaginary unit
+    const FieldVector<rt,d> one(1.); // vector filled with 1
+    const FieldVector<ct,d> iVec(ct(0.,1.)); // vector filled with I
+
+    std::cout << __func__ << "\t \t ( " << Dune::className(one) << " and " << Dune::className(iVec) << ")" << std::endl;
+
+    const bool isRealOne = Dune::is_same<typename Dune::FieldTraits<rt>::field_type,typename Dune::FieldTraits<rt>::real_type>::value;
+    const bool isRealIVec = Dune::is_same<typename Dune::FieldTraits<ct>::field_type,typename Dune::FieldTraits<ct>::real_type> ::value;
+    dune_static_assert(isRealOne,"1-vector expected to be real");
+    dune_static_assert(!isRealIVec,"i-vector expected to be complex");
+
+    ct result = ct();
+    ct length = ct(d);
+
+
+    // one^H*one should equal d
+    result = dot(one,one);
+    assert(std::abs(result-length)<= myEps);
+    result = one.dot(one);
+    assert(std::abs(result-length)<= myEps);
+
+
+    // iVec^H*iVec should equal d
+    result = dot(iVec,iVec);
+    assert(std::abs(result-length)<= myEps);
+    result = iVec.dot(iVec);
+    assert(std::abs(result-length)<= myEps);
+
+
+    // test that we do conjugate first argument
+    result = dot(one,iVec);
+    assert(std::abs(result-length*I)<= myEps);
+    result = dot(one,iVec);
+    assert(std::abs(result-length*I)<= myEps);
+
+
+    // test that we do not conjugate second argument
+    result = dot(iVec,one);
+    assert(std::abs(result+length*I)<= myEps);
+    result = iVec.dot(one);
+    assert(std::abs(result+length*I)<= myEps);
+
+
+    // test that dotT does not conjugate at all
+    result = dotT(one,one) + one*one;
+    assert(std::abs(result-ct(2)*length)<= myEps);
+    result = dotT(iVec,iVec) + iVec*iVec;
+    assert(std::abs(result+ct(2)*length)<= myEps);
+    result = dotT(one,iVec) + one*iVec;
+    assert(std::abs(result-ct(2)*length*I)<= myEps);
+    result = dotT(iVec,one) + iVec*one;
+    assert(std::abs(result-ct(2)*length*I)<= myEps);
+
+  }
+
+};
+
 template<class ft, int d>
 struct FieldVectorTest
 {
@@ -209,6 +280,7 @@ struct FieldVectorTest
     // --- test complex and real valued vectors
     FieldVectorMainTest<ft,ft,d>();
     FieldVectorMainTest<complex<ft>,ft,d>();
+    DotProductTest<ft,d>();
     // --- test next lower dimension
     FieldVectorTest<ft,d-1>();
   }
@@ -225,6 +297,7 @@ public:
     FieldVectorMainTest<ft,ft,1>();
     ScalarOperatorTest<ft>();
     ScalarOrderingTest<ft>();
+    DotProductTest<ft,1>();
     // --- complex valued
     FieldVectorMainTest<complex<ft>,ft,1>();
     ScalarOperatorTest< complex<ft> >();