From aeec5e67ffe93a182050562b0155caf5727e22db Mon Sep 17 00:00:00 2001
From: Denis Davydov <davydden@gmail.com>
Date: Tue, 16 Jun 2015 23:09:05 +0200
Subject: [PATCH] add real/imag to PETSc vector. Adjust description of
 operator* for complex case

---
 include/deal.II/lac/petsc_vector_base.h | 42 ++++++++++++++++++++++++-
 source/lac/petsc_vector_base.cc         |  6 +++-
 2 files changed, 46 insertions(+), 2 deletions(-)

diff --git a/include/deal.II/lac/petsc_vector_base.h b/include/deal.II/lac/petsc_vector_base.h
index 38412c37ab..6765f92117 100644
--- a/include/deal.II/lac/petsc_vector_base.h
+++ b/include/deal.II/lac/petsc_vector_base.h
@@ -136,6 +136,19 @@ namespace PETScWrappers
        */
       const VectorReference &operator /= (const PetscScalar &s) const;
 
+      /**
+       * Return the real part of the value of the referenced element.
+       */
+      PetscReal real () const;
+
+      /**
+       * Return the imaginary part of the value of the referenced element.
+       *
+       * @note This operation is not defined for real numbers and an
+       * exception is thrown.
+       */
+      PetscReal imag () const;
+
       /**
        * Convert the reference to an actual value, i.e. return the value of
        * the referenced element of the vector.
@@ -443,6 +456,8 @@ namespace PETScWrappers
     /**
      * Return the scalar product of two vectors. The vectors must have the
      * same size.
+     *
+     * For complex valued vector, this gives$\left(v^\ast,vec\right)$.
      */
     PetscScalar operator * (const VectorBase &vec) const;
 
@@ -1022,10 +1037,35 @@ namespace PETScWrappers
 
       return *this;
     }
-  }
 
 
 
+    inline
+    PetscReal
+    VectorReference::real () const
+    {
+#ifndef PETSC_USE_COMPLEX
+      return static_cast<PetscScalar>(*this);
+#else
+      return PetscRealPart (static_cast<PetscScalar>(*this));
+#endif
+    }
+
+
+
+    inline
+    PetscReal
+    VectorReference::imag () const
+    {
+#ifndef PETSC_USE_COMPLEX
+      return PetscReal (0);
+#else
+      return PetscImaginaryPart (static_cast<PetscScalar>(*this));
+#endif
+    }
+
+  } // namespace internal
+
   inline
   bool
   VectorBase::in_local_range (const size_type index) const
diff --git a/source/lac/petsc_vector_base.cc b/source/lac/petsc_vector_base.cc
index 0feeaef91a..e3191cd202 100644
--- a/source/lac/petsc_vector_base.cc
+++ b/source/lac/petsc_vector_base.cc
@@ -389,7 +389,11 @@ namespace PETScWrappers
 
     PetscScalar result;
 
-    const int ierr = VecDot (vector, vec.vector, &result);
+    //For complex vectors, VecDot() computes
+    //    val = (x,y) = y^H x,
+    //where y^H denotes the conjugate transpose of y.
+    //Note that this corresponds to the usual "mathematicians" complex inner product where the SECOND argument gets the complex conjugate.
+    const int ierr = VecDot (vec.vector, vector, &result);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
 
     return result;
-- 
2.39.5