From aeec5e67ffe93a182050562b0155caf5727e22db Mon Sep 17 00:00:00 2001 From: Denis Davydov 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(*this); +#else + return PetscRealPart (static_cast(*this)); +#endif + } + + + + inline + PetscReal + VectorReference::imag () const + { +#ifndef PETSC_USE_COMPLEX + return PetscReal (0); +#else + return PetscImaginaryPart (static_cast(*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