Merge branch 'complex-coef-dev' of github.com:mfem/mfem into complex-coef-dev

examples/ex22.cpp

@@ -62,14 +62,9 @@ static real_t epsilon_ = 1.0;
 static real_t sigma_ = 20.0;
 static real_t omega_ = 10.0;
 
-real_t u0_real_exact(const Vector &);
-real_t u0_imag_exact(const Vector &);
-
-void u1_real_exact(const Vector &, Vector &);
-void u1_imag_exact(const Vector &, Vector &);
-
-void u2_real_exact(const Vector &, Vector &);
-void u2_imag_exact(const Vector &, Vector &);
+complex<real_t> u0_exact(const Vector &x);
+void u1_exact(const Vector &, ComplexVector &);
+void u2_exact(const Vector &, ComplexVector &);
 
 bool check_for_inline_mesh(const char * mesh_file);
 
@@ -215,54 +210,48 @@ int main(int argc, char *argv[])
    ComplexGridFunction * u_exact = NULL;
    if (exact_sol) { u_exact = new ComplexGridFunction(fespace); }
 
-   FunctionCoefficient u0_r(u0_real_exact);
-   FunctionCoefficient u0_i(u0_imag_exact);
-   VectorFunctionCoefficient u1_r(dim, u1_real_exact);
-   VectorFunctionCoefficient u1_i(dim, u1_imag_exact);
-   VectorFunctionCoefficient u2_r(dim, u2_real_exact);
-   VectorFunctionCoefficient u2_i(dim, u2_imag_exact);
+   ComplexFunctionCoefficient u0(u0_exact);
+   ComplexVectorFunctionCoefficient u1(dim, u1_exact);
+   ComplexVectorFunctionCoefficient u2(dim, u2_exact);
 
-   ConstantCoefficient zeroCoef(0.0);
-   ConstantCoefficient oneCoef(1.0);
+   ComplexConstantCoefficient oneCoef(1.0);
 
-   Vector zeroVec(dim); zeroVec = 0.0;
    Vector  oneVec(dim);  oneVec = 0.0; oneVec[(prob==2)?(dim-1):0] = 1.0;
-   VectorConstantCoefficient zeroVecCoef(zeroVec);
-   VectorConstantCoefficient oneVecCoef(oneVec);
+   ComplexVectorConstantCoefficient oneVecCoef(oneVec);
 
    switch (prob)
    {
       case 0:
          if (exact_sol)
          {
-            u.ProjectBdrCoefficient(u0_r, u0_i, ess_bdr);
-            u_exact->ProjectCoefficient(u0_r, u0_i);
+            u.ProjectBdrCoefficient(u0, ess_bdr);
+            u_exact->ProjectCoefficient(u0);
          }
          else
          {
-            u.ProjectBdrCoefficient(oneCoef, zeroCoef, ess_bdr);
+            u.ProjectBdrCoefficient(oneCoef, ess_bdr);
          }
          break;
       case 1:
          if (exact_sol)
          {
-            u.ProjectBdrCoefficientTangent(u1_r, u1_i, ess_bdr);
-            u_exact->ProjectCoefficient(u1_r, u1_i);
+            u.ProjectBdrCoefficientTangent(u1, ess_bdr);
+            u_exact->ProjectCoefficient(u1);
          }
          else
          {
-            u.ProjectBdrCoefficientTangent(oneVecCoef, zeroVecCoef, ess_bdr);
+            u.ProjectBdrCoefficientTangent(oneVecCoef, ess_bdr);
          }
          break;
       case 2:
          if (exact_sol)
          {
-            u.ProjectBdrCoefficientNormal(u2_r, u2_i, ess_bdr);
-            u_exact->ProjectCoefficient(u2_r, u2_i);
+            u.ProjectBdrCoefficientNormal(u2, ess_bdr);
+            u_exact->ProjectCoefficient(u2);
          }
          else
          {
-            u.ProjectBdrCoefficientNormal(oneVecCoef, zeroVecCoef, ess_bdr);
+            u.ProjectBdrCoefficientNormal(oneVecCoef, ess_bdr);
          }
          break;
       default: break; // This should be unreachable
@@ -300,27 +289,24 @@ int main(int argc, char *argv[])
    ConstantCoefficient lossCoef(omega_ * sigma_);
    ConstantCoefficient negMassCoef(omega_ * omega_ * epsilon_);
 
+   ComplexConstantCoefficient complexMassCoef(-omega_ * omega_ * epsilon_,
+                                              omega_ * sigma_);
+
    SesquilinearForm *a = new SesquilinearForm(fespace, conv);
    if (pa) { a->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
    switch (prob)
    {
       case 0:
-         a->AddDomainIntegrator(new DiffusionIntegrator(stiffnessCoef),
-                                NULL);
-         a->AddDomainIntegrator(new MassIntegrator(massCoef),
-                                new MassIntegrator(lossCoef));
+         a->AddDomainIntegrator<DiffusionIntegrator>(stiffnessCoef);
+         a->AddDomainIntegrator<MassIntegrator>(complexMassCoef);
          break;
       case 1:
-         a->AddDomainIntegrator(new CurlCurlIntegrator(stiffnessCoef),
-                                NULL);
-         a->AddDomainIntegrator(new VectorFEMassIntegrator(massCoef),
-                                new VectorFEMassIntegrator(lossCoef));
+         a->AddDomainIntegrator<CurlCurlIntegrator>(stiffnessCoef);
+         a->AddDomainIntegrator<VectorFEMassIntegrator>(complexMassCoef);
          break;
       case 2:
-         a->AddDomainIntegrator(new DivDivIntegrator(stiffnessCoef),
-                                NULL);
-         a->AddDomainIntegrator(new VectorFEMassIntegrator(massCoef),
-                                new VectorFEMassIntegrator(lossCoef));
+         a->AddDomainIntegrator<DivDivIntegrator>(stiffnessCoef);
+         a->AddDomainIntegrator<VectorFEMassIntegrator>(complexMassCoef);
          break;
       default: break; // This should be unreachable
    }
@@ -436,29 +422,24 @@ int main(int argc, char *argv[])
 
    if (exact_sol)
    {
-      real_t err_r = -1.0;
-      real_t err_i = -1.0;
+      real_t err_u = -1.0;
 
       switch (prob)
       {
          case 0:
-            err_r = u.real().ComputeL2Error(u0_r);
-            err_i = u.imag().ComputeL2Error(u0_i);
+            err_u = u.ComputeL2Error(u0);
             break;
          case 1:
-            err_r = u.real().ComputeL2Error(u1_r);
-            err_i = u.imag().ComputeL2Error(u1_i);
+            err_u = u.ComputeL2Error(u1);
             break;
          case 2:
-            err_r = u.real().ComputeL2Error(u2_r);
-            err_i = u.imag().ComputeL2Error(u2_i);
+            err_u = u.ComputeL2Error(u2);
             break;
          default: break; // This should be unreachable
       }
 
       cout << endl;
-      cout << "|| Re (u_h - u) ||_{L^2} = " << err_r << endl;
-      cout << "|| Im (u_h - u) ||_{L^2} = " << err_i << endl;
+      cout << "|| u_h - u ||_{L^2} = " << err_u << endl;
       cout << endl;
    }
 
@@ -564,36 +545,14 @@ complex<real_t> u0_exact(const Vector &x)
    return std::exp(-i * kappa * x[dim - 1]);
 }
 
-real_t u0_real_exact(const Vector &x)
-{
-   return u0_exact(x).real();
-}
-
-real_t u0_imag_exact(const Vector &x)
-{
-   return u0_exact(x).imag();
-}
-
-void u1_real_exact(const Vector &x, Vector &v)
+void u1_exact(const Vector &x, ComplexVector &v)
 {
    int dim = x.Size();
-   v.SetSize(dim); v = 0.0; v[0] = u0_real_exact(x);
+   v.SetSize(dim); v = 0.0; v[0] = u0_exact(x);
 }
 
-void u1_imag_exact(const Vector &x, Vector &v)
+void u2_exact(const Vector &x, ComplexVector &v)
 {
    int dim = x.Size();
-   v.SetSize(dim); v = 0.0; v[0] = u0_imag_exact(x);
-}
-
-void u2_real_exact(const Vector &x, Vector &v)
-{
-   int dim = x.Size();
-   v.SetSize(dim); v = 0.0; v[dim-1] = u0_real_exact(x);
-}
-
-void u2_imag_exact(const Vector &x, Vector &v)
-{
-   int dim = x.Size();
-   v.SetSize(dim); v = 0.0; v[dim-1] = u0_imag_exact(x);
+   v.SetSize(dim); v = 0.0; v[dim-1] = u0_exact(x);
 }

examples/ex22p.cpp

@@ -62,6 +62,10 @@ static real_t epsilon_ = 1.0;
 static real_t sigma_ = 20.0;
 static real_t omega_ = 10.0;
 
+complex<real_t> u0_exact(const Vector &x);
+void u1_exact(const Vector &, ComplexVector &);
+void u2_exact(const Vector &, ComplexVector &);
+
 real_t u0_real_exact(const Vector &);
 real_t u0_imag_exact(const Vector &);
 
@@ -244,13 +248,22 @@ int main(int argc, char *argv[])
    ParComplexGridFunction * u_exact = NULL;
    if (exact_sol) { u_exact = new ParComplexGridFunction(fespace); }
 
+   ComplexFunctionCoefficient u0(u0_exact);
+   ComplexVectorFunctionCoefficient u1(dim, u1_exact);
+   ComplexVectorFunctionCoefficient u2(dim, u2_exact);
+
+   ComplexConstantCoefficient oneCoef(1.0);
+
+   Vector  oneVec(dim);  oneVec = 0.0; oneVec[(prob==2)?(dim-1):0] = 1.0;
+   ComplexVectorConstantCoefficient oneVecCoef(oneVec);
+
    FunctionCoefficient u0_r(u0_real_exact);
    FunctionCoefficient u0_i(u0_imag_exact);
    VectorFunctionCoefficient u1_r(dim, u1_real_exact);
    VectorFunctionCoefficient u1_i(dim, u1_imag_exact);
    VectorFunctionCoefficient u2_r(dim, u2_real_exact);
    VectorFunctionCoefficient u2_i(dim, u2_imag_exact);
-
+   /*
    ConstantCoefficient zeroCoef(0.0);
    ConstantCoefficient oneCoef(1.0);
 
@@ -258,40 +271,40 @@ int main(int argc, char *argv[])
    Vector  oneVec(dim);  oneVec = 0.0; oneVec[(prob==2)?(dim-1):0] = 1.0;
    VectorConstantCoefficient zeroVecCoef(zeroVec);
    VectorConstantCoefficient oneVecCoef(oneVec);
-
+   */
    switch (prob)
    {
       case 0:
          if (exact_sol)
          {
-            u.ProjectBdrCoefficient(u0_r, u0_i, ess_bdr);
-            u_exact->ProjectCoefficient(u0_r, u0_i);
+            u.ProjectBdrCoefficient(u0, ess_bdr);
+            u_exact->ProjectCoefficient(u0);
          }
          else
          {
-            u.ProjectBdrCoefficient(oneCoef, zeroCoef, ess_bdr);
+            u.ProjectBdrCoefficient(oneCoef, ess_bdr);
          }
          break;
       case 1:
          if (exact_sol)
          {
-            u.ProjectBdrCoefficientTangent(u1_r, u1_i, ess_bdr);
-            u_exact->ProjectCoefficient(u1_r, u1_i);
+            u.ProjectBdrCoefficientTangent(u1, ess_bdr);
+            u_exact->ProjectCoefficient(u1);
          }
          else
          {
-            u.ProjectBdrCoefficientTangent(oneVecCoef, zeroVecCoef, ess_bdr);
+            u.ProjectBdrCoefficientTangent(oneVecCoef, ess_bdr);
          }
          break;
       case 2:
          if (exact_sol)
          {
-            u.ProjectBdrCoefficientNormal(u2_r, u2_i, ess_bdr);
-            u_exact->ProjectCoefficient(u2_r, u2_i);
+            u.ProjectBdrCoefficientNormal(u2, ess_bdr);
+            u_exact->ProjectCoefficient(u2);
          }
          else
          {
-            u.ProjectBdrCoefficientNormal(oneVecCoef, zeroVecCoef, ess_bdr);
+            u.ProjectBdrCoefficientNormal(oneVecCoef, ess_bdr);
          }
          break;
       default: break; // This should be unreachable
@@ -331,27 +344,24 @@ int main(int argc, char *argv[])
    ConstantCoefficient lossCoef(omega_ * sigma_);
    ConstantCoefficient negMassCoef(omega_ * omega_ * epsilon_);
 
+   ComplexConstantCoefficient complexMassCoef(-omega_ * omega_ * epsilon_,
+                                              omega_ * sigma_);
+
    ParSesquilinearForm *a = new ParSesquilinearForm(fespace, conv);
    if (pa) { a->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
    switch (prob)
    {
       case 0:
-         a->AddDomainIntegrator(new DiffusionIntegrator(stiffnessCoef),
-                                NULL);
-         a->AddDomainIntegrator(new MassIntegrator(massCoef),
-                                new MassIntegrator(lossCoef));
+         a->AddDomainIntegrator<DiffusionIntegrator>(stiffnessCoef);
+         a->AddDomainIntegrator<MassIntegrator>(complexMassCoef);
          break;
       case 1:
-         a->AddDomainIntegrator(new CurlCurlIntegrator(stiffnessCoef),
-                                NULL);
-         a->AddDomainIntegrator(new VectorFEMassIntegrator(massCoef),
-                                new VectorFEMassIntegrator(lossCoef));
+         a->AddDomainIntegrator<CurlCurlIntegrator>(stiffnessCoef);
+         a->AddDomainIntegrator<VectorFEMassIntegrator>(complexMassCoef);
          break;
       case 2:
-         a->AddDomainIntegrator(new DivDivIntegrator(stiffnessCoef),
-                                NULL);
-         a->AddDomainIntegrator(new VectorFEMassIntegrator(massCoef),
-                                new VectorFEMassIntegrator(lossCoef));
+         a->AddDomainIntegrator<DivDivIntegrator>(stiffnessCoef);
+         a->AddDomainIntegrator<VectorFEMassIntegrator>(complexMassCoef);
          break;
       default: break; // This should be unreachable
    }
@@ -475,22 +485,18 @@ int main(int argc, char *argv[])
 
    if (exact_sol)
    {
-      real_t err_r = -1.0;
-      real_t err_i = -1.0;
+      real_t err_u = -1.0;
 
       switch (prob)
       {
          case 0:
-            err_r = u.real().ComputeL2Error(u0_r);
-            err_i = u.imag().ComputeL2Error(u0_i);
+            err_u = u.ComputeL2Error(u0);
             break;
          case 1:
-            err_r = u.real().ComputeL2Error(u1_r);
-            err_i = u.imag().ComputeL2Error(u1_i);
+            err_u = u.ComputeL2Error(u1);
             break;
          case 2:
-            err_r = u.real().ComputeL2Error(u2_r);
-            err_i = u.imag().ComputeL2Error(u2_i);
+            err_u = u.ComputeL2Error(u2);
             break;
          default: break; // This should be unreachable
       }
@@ -498,8 +504,7 @@ int main(int argc, char *argv[])
       if ( myid == 0 )
       {
          cout << endl;
-         cout << "|| Re (u_h - u) ||_{L^2} = " << err_r << endl;
-         cout << "|| Im (u_h - u) ||_{L^2} = " << err_i << endl;
+         cout << "|| u_h - u ||_{L^2} = " << err_u << endl;
          cout << endl;
       }
    }
@@ -627,6 +632,12 @@ real_t u0_imag_exact(const Vector &x)
    return u0_exact(x).imag();
 }
 
+void u1_exact(const Vector &x, ComplexVector &v)
+{
+   int dim = x.Size();
+   v.SetSize(dim); v = 0.0; v[0] = u0_exact(x);
+}
+
 void u1_real_exact(const Vector &x, Vector &v)
 {
    int dim = x.Size();
@@ -639,6 +650,12 @@ void u1_imag_exact(const Vector &x, Vector &v)
    v.SetSize(dim); v = 0.0; v[0] = u0_imag_exact(x);
 }
 
+void u2_exact(const Vector &x, ComplexVector &v)
+{
+   int dim = x.Size();
+   v.SetSize(dim); v = 0.0; v[dim-1] = u0_exact(x);
+}
+
 void u2_real_exact(const Vector &x, Vector &v)
 {
    int dim = x.Size();

fem/CMakeLists.txt

@@ -59,6 +59,7 @@ set(SRCS
   integ/nonlininteg_vecconvection_pa.cpp
   integ/nonlininteg_vecconvection_mf.cpp
   coefficient.cpp
+  complex_coefficient.cpp
   complex_fem.cpp
   convergence.cpp
   datacollection.cpp
@@ -176,6 +177,7 @@ set(HDRS
   integ/bilininteg_hcurlhdiv_kernels.hpp
   integ/bilininteg_mass_kernels.hpp
   coefficient.hpp
+  complex_coefficient.hpp
   complex_fem.hpp
   convergence.hpp
   datacollection.hpp

fem/complex_coefficient.cpp

@@ -0,0 +1,217 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "complex_fem.hpp"
+#include "../general/forall.hpp"
+
+using namespace std;
+
+namespace mfem
+{
+
+real_t
+RealPartCoefficient::Eval(ElementTransformation &T,
+                          const IntegrationPoint &ip)
+{
+   complex_t val = complex_coef_.Eval(T, ip);
+   return val.real();
+}
+
+real_t
+ImagPartCoefficient::Eval(ElementTransformation &T,
+                          const IntegrationPoint &ip)
+{
+   complex_t val = complex_coef_.Eval(T, ip);
+   return val.imag();
+}
+
+RealPartVectorCoefficient::RealPartVectorCoefficient(ComplexVectorCoefficient &
+                                                     complex_vcoef)
+   : VectorCoefficient(complex_vcoef.GetVDim()),
+     complex_vcoef_(complex_vcoef),
+     val_(vdim)
+{}
+
+void
+RealPartVectorCoefficient::Eval(Vector &V, ElementTransformation &T,
+                                const IntegrationPoint &ip)
+{
+   complex_vcoef_.Eval(val_, T, ip);
+   V = val_.real();
+}
+
+ImagPartVectorCoefficient::ImagPartVectorCoefficient(ComplexVectorCoefficient &
+                                                     complex_vcoef)
+   : VectorCoefficient(complex_vcoef.GetVDim()),
+     complex_vcoef_(complex_vcoef),
+     val_(vdim)
+{}
+
+void
+ImagPartVectorCoefficient::Eval(Vector &V, ElementTransformation &T,
+                                const IntegrationPoint &ip)
+{
+   complex_vcoef_.Eval(val_, T, ip);
+   V = val_.imag();
+}
+
+RealPartMatrixCoefficient::RealPartMatrixCoefficient(ComplexMatrixCoefficient &
+                                                     complex_mcoef)
+   : MatrixCoefficient(complex_mcoef.GetHeight(), complex_mcoef.GetWidth()),
+     complex_mcoef_(complex_mcoef),
+     val_(height, width)
+{}
+
+void
+RealPartMatrixCoefficient::Eval(DenseMatrix &M, ElementTransformation &T,
+                                const IntegrationPoint &ip)
+{
+   complex_mcoef_.Eval(val_, T, ip);
+   M = val_.real();
+}
+
+ImagPartMatrixCoefficient::ImagPartMatrixCoefficient(ComplexMatrixCoefficient &
+                                                     complex_mcoef)
+   : MatrixCoefficient(complex_mcoef.GetHeight(), complex_mcoef.GetWidth()),
+     complex_mcoef_(complex_mcoef),
+     val_(height, width)
+{}
+
+void
+ImagPartMatrixCoefficient::Eval(DenseMatrix &M, ElementTransformation &T,
+                                const IntegrationPoint &ip)
+{
+   complex_mcoef_.Eval(val_, T, ip);
+   M = val_.imag();
+}
+
+ComplexCoefficient::ComplexCoefficient()
+   : time(0.),
+     re_part_coef_(*this), im_part_coef_(*this),
+     real_coef_(re_part_coef_), imag_coef_(im_part_coef_)
+{ }
+
+ComplexCoefficient::ComplexCoefficient(Coefficient &c_r,
+                                       Coefficient &c_i)
+   : time(c_r.GetTime()),
+     re_part_coef_(*this), im_part_coef_(*this),
+     real_coef_(c_r), imag_coef_(c_i)
+{
+   c_i.SetTime(time);
+}
+
+complex_t
+ComplexCoefficient::Eval(ElementTransformation &T,
+                         const IntegrationPoint &ip)
+{
+   // Avoid circular dependency
+   MFEM_VERIFY(std::addressof(real_coef_) != std::addressof(re_part_coef_) &&
+               std::addressof(imag_coef_) != std::addressof(im_part_coef_),
+               "Classes dervied from ComplexCoefficient must either "
+               "implement an Eval method or supply Coefficients "
+               "for both the real and imaginary parts of the field.");
+
+   return complex_t(real_coef_.Eval(T, ip), imag_coef_.Eval(T, ip));
+}
+
+ComplexVectorCoefficient::ComplexVectorCoefficient(VectorCoefficient &v_r,
+                                                   VectorCoefficient &v_i)
+   : vdim(v_r.GetVDim()), time(v_r.GetTime()),
+     re_part_vcoef_(*this), im_part_vcoef_(*this),
+     real_vcoef_(v_r), imag_vcoef_(v_i)
+{
+   MFEM_ASSERT(v_r.GetVDim() == v_i.GetVDim(), "ComplexVectorCoefficient"
+               " - incompatible vector dimensions of real and imaginary parts.");
+
+   v_i.SetTime(time);
+}
+
+void ComplexVectorCoefficient::Eval(ComplexVector &V, ElementTransformation &T,
+                                    const IntegrationPoint &ip)
+{
+   // Avoid circular dependency
+   MFEM_VERIFY(std::addressof(real_vcoef_) != std::addressof(re_part_vcoef_) &&
+               std::addressof(imag_vcoef_) != std::addressof(im_part_vcoef_),
+               "Classes dervied from ComplexVectorCoefficient must either "
+               "implement an Eval method or supply VectorCoefficients "
+               "for both the real and imaginary parts of the field.");
+
+   V_r_.SetSize(vdim);
+   V_i_.SetSize(vdim);
+
+   real_vcoef_.Eval(V_r_, T, ip);
+   imag_vcoef_.Eval(V_i_, T, ip);
+
+   V.Set(V_r_, V_i_);
+}
+
+ComplexConstantCoefficient::ComplexConstantCoefficient(
+   const complex_t z)
+   : val(z), real_coef(z.real()), imag_coef(z.imag())
+{
+   real_coef_ = real_coef;
+   imag_coef_ = imag_coef;
+}
+
+ComplexConstantCoefficient::ComplexConstantCoefficient(
+   real_t z_r, real_t z_i)
+   : real_coef(z_r), imag_coef(z_i)
+{
+   val = complex_t(z_r, z_i);
+
+   real_coef_ = real_coef;
+   imag_coef_ = imag_coef;
+}
+
+complex_t ComplexFunctionCoefficient::Eval(ElementTransformation & T,
+                                           const IntegrationPoint & ip)
+{
+   real_t x[3];
+   Vector transip(x, 3);
+
+   T.Transform(ip, transip);
+
+   if (Function)
+   {
+      return Function(transip);
+   }
+   else
+   {
+      return TDFunction(transip, GetTime());
+   }
+}
+
+void ComplexVectorFunctionCoefficient::Eval(ComplexVector &V,
+                                            ElementTransformation &T,
+                                            const IntegrationPoint &ip)
+{
+   real_t x[3];
+   Vector transip(x, 3);
+
+   T.Transform(ip, transip);
+
+   V.SetSize(vdim);
+   if (Function)
+   {
+      Function(transip, V);
+   }
+   else
+   {
+      TDFunction(transip, GetTime(), V);
+   }
+   if (Q)
+   {
+      V *= Q->Eval(T, ip, GetTime());
+   }
+}
+
+} // end namespace mfem
+

fem/complex_coefficient.hpp

@@ -0,0 +1,523 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_COMPLEX_COEFFICIENT
+#define MFEM_COMPLEX_COEFFICIENT
+
+#include "../config/config.hpp"
+#include "../linalg/linalg.hpp"
+#include "coefficient.hpp"
+#include "intrules.hpp"
+#include "eltrans.hpp"
+
+namespace mfem
+{
+
+class ComplexCoefficient;
+class ComplexVectorCoefficient;
+class ComplexMatrixCoefficient;
+
+/// Standard Coefficient which returns the real part of a ComplexCoefficient
+class RealPartCoefficient : public Coefficient
+{
+private:
+   ComplexCoefficient &complex_coef_;
+
+public:
+   RealPartCoefficient(ComplexCoefficient & complex_coef)
+      : complex_coef_(complex_coef) {}
+
+   real_t Eval(ElementTransformation &T,
+               const IntegrationPoint &ip);
+};
+
+/// Standard Coefficient which returns the imaginary part of a
+/// ComplexCoefficient
+class ImagPartCoefficient : public Coefficient
+{
+private:
+   ComplexCoefficient &complex_coef_;
+
+public:
+   ImagPartCoefficient(ComplexCoefficient & complex_coef)
+      : complex_coef_(complex_coef) {}
+
+   real_t Eval(ElementTransformation &T,
+               const IntegrationPoint &ip);
+};
+
+typedef ImagPartCoefficient ImaginaryPartCoefficient;
+
+class RealPartVectorCoefficient : public VectorCoefficient
+{
+private:
+   ComplexVectorCoefficient &complex_vcoef_;
+   mutable ComplexVector val_;
+
+public:
+   RealPartVectorCoefficient(ComplexVectorCoefficient & complex_vcoef);
+
+   void Eval(Vector &V, ElementTransformation &T,
+             const IntegrationPoint &ip);
+};
+
+class ImagPartVectorCoefficient : public VectorCoefficient
+{
+private:
+   ComplexVectorCoefficient &complex_vcoef_;
+   mutable ComplexVector val_;
+
+public:
+   ImagPartVectorCoefficient(ComplexVectorCoefficient & complex_vcoef);
+
+   void Eval(Vector &V, ElementTransformation &T,
+             const IntegrationPoint &ip);
+};
+
+typedef ImagPartVectorCoefficient ImaginaryPartVectorCoefficient;
+
+class RealPartMatrixCoefficient : public MatrixCoefficient
+{
+private:
+   ComplexMatrixCoefficient &complex_mcoef_;
+   mutable ComplexTypeDenseMatrix val_;
+
+public:
+   RealPartMatrixCoefficient(ComplexMatrixCoefficient & complex_mcoef);
+
+   void Eval(DenseMatrix &M, ElementTransformation &T,
+             const IntegrationPoint &ip);
+};
+
+class ImagPartMatrixCoefficient : public MatrixCoefficient
+{
+private:
+   ComplexMatrixCoefficient &complex_mcoef_;
+   mutable ComplexTypeDenseMatrix val_;
+
+public:
+   ImagPartMatrixCoefficient(ComplexMatrixCoefficient & complex_mcoef);
+
+   void Eval(DenseMatrix &V, ElementTransformation &T,
+             const IntegrationPoint &ip);
+};
+
+typedef ImagPartMatrixCoefficient ImaginaryPartMatrixCoefficient;
+
+/** @brief Base class ComplexCoefficients that optionally depend on space and
+    time. These are used by the SesquilinearForm, ComplexLinearForm, and
+    ComplexGridFunction classes to represent the physical coefficients in
+    the PDEs that are being discretized. This class can also be used in a more
+    general way to represent functions that don't necessarily belong to a FE
+    space, e.g., to project onto ComplexGridFunctions to use as initial
+    conditions, exact solutions, etc. See, e.g., ex22 for these uses. */
+class ComplexCoefficient
+{
+protected:
+   real_t time;
+
+private:
+   RealPartCoefficient re_part_coef_;
+   ImagPartCoefficient im_part_coef_;
+
+protected:
+   Coefficient &real_coef_;
+   Coefficient &imag_coef_;
+
+public:
+
+   ComplexCoefficient();
+   ComplexCoefficient(Coefficient &c_r, Coefficient &c_i);
+
+   /// Set the time for time dependent coefficients
+   virtual void SetTime(real_t t)
+   { time = t; real_coef_.SetTime(t); imag_coef_.SetTime(t); }
+
+   /// Get the time for time dependent coefficients
+   real_t GetTime() { return time; }
+
+   /** @brief Evaluate the coefficient in the element described by @a T at the
+       point @a ip. */
+   /** @note When this method is called, the caller must make sure that the
+       IntegrationPoint associated with @a T is the same as @a ip. This can be
+       achieved by calling T.SetIntPoint(&ip). */
+   virtual complex_t Eval(ElementTransformation &T,
+                          const IntegrationPoint &ip);
+
+   /** @brief Evaluate the coefficient in the element described by @a T at the
+       point @a ip at time @a t. */
+   /** @note When this method is called, the caller must make sure that the
+       IntegrationPoint associated with @a T is the same as @a ip. This can be
+       achieved by calling T.SetIntPoint(&ip). */
+   complex_t Eval(ElementTransformation &T,
+                  const IntegrationPoint &ip, real_t t)
+   {
+      SetTime(t);
+      return Eval(T, ip);
+   }
+
+   /** @brief Access a standard Coefficient object reproducing the real part of
+       the complex-valued field */
+   /** @note By default this method returns an internal object which
+       computes the complex value using the above Eval method and
+       returns its real part. Custom implementations may choose to
+       override this method with a more efficient real-valued
+       coefficient. */
+   virtual Coefficient & real() { return real_coef_; }
+
+   /** @brief Access a standard Coefficient object reproducing the imaginary
+       part of the complex-valued field */
+   /** @note By default this method returns an internal object which
+       computes the complex value using the above Eval method and
+       returns its imaginary part. Custom implementations may choose to
+       override this method with a more efficient real-valued
+       coefficient. */
+   virtual Coefficient & imag() { return imag_coef_; }
+
+   virtual ~ComplexCoefficient() { }
+};
+
+/** @brief Base class ComplexVectorCoefficients that optionally depend
+    on space and time. These are used by the SesquilinearForm,
+    ComplexLinearForm, and ComplexGridFunction classes to represent
+    the physical vector-valued coefficients in the PDEs that are being
+    discretized. This class can also be used in a more general way to
+    represent functions that don't necessarily belong to a FE space,
+    e.g., to project onto ComplexGridFunctions to use as initial
+    conditions, exact solutions, etc. See, e.g., ex22 for these
+    uses. */
+class ComplexVectorCoefficient
+{
+protected:
+   int vdim;
+   real_t time;
+
+private:
+   RealPartVectorCoefficient re_part_vcoef_;
+   ImagPartVectorCoefficient im_part_vcoef_;
+
+protected:
+   VectorCoefficient &real_vcoef_;
+   VectorCoefficient &imag_vcoef_;
+
+   mutable Vector V_r_;
+   mutable Vector V_i_;
+
+public:
+   ComplexVectorCoefficient(int vd)
+      : vdim(vd), time(0.),
+        re_part_vcoef_(*this), im_part_vcoef_(*this),
+        real_vcoef_(re_part_vcoef_), imag_vcoef_(im_part_vcoef_)
+   { }
+
+   ComplexVectorCoefficient(VectorCoefficient &v_r, VectorCoefficient &v_i);
+
+
+   /// Set the time for time dependent coefficients
+   virtual void SetTime(real_t t)
+   { time = t; real_vcoef_.SetTime(t); imag_vcoef_.SetTime(t); }
+
+   /// Get the time for time dependent coefficients
+   real_t GetTime() { return time; }
+
+   /// Returns dimension of the vector.
+   int GetVDim() { return vdim; }
+
+   /** @brief Evaluate the vector coefficient in the element described by @a T
+       at the point @a ip, storing the result in @a V. */
+   /** @note When this method is called, the caller must make sure that the
+       IntegrationPoint associated with @a T is the same as @a ip. This can be
+       achieved by calling T.SetIntPoint(&ip). */
+   virtual void Eval(ComplexVector &V, ElementTransformation &T,
+                     const IntegrationPoint &ip);
+
+   /** @brief Evaluate the vector coefficient in the element described by @a T
+       at the point @a ip at time @a t, storing the result in @a V. */
+   /** @note When this method is called, the caller must make sure that the
+       IntegrationPoint associated with @a T is the same as @a ip. This can be
+       achieved by calling T.SetIntPoint(&ip). */
+   void Eval(ComplexVector &V, ElementTransformation &T,
+             const IntegrationPoint &ip, real_t t)
+   {
+      SetTime(t);
+      Eval(V, T, ip);
+   }
+
+   /** @brief Access a standard Coefficient object reproducing the real part of
+       the complex-valued field */
+   /** @note By default this method returns an internal object which
+       computes the complex value using the above Eval method and
+       returns its real part. Custom implementations may choose to
+       override this method with a more efficient real-valued
+       coefficient. */
+   virtual VectorCoefficient & real() { return real_vcoef_; }
+
+   /** @brief Access a standard Coefficient object reproducing the imaginary
+       part of the complex-valued field */
+   /** @note By default this method returns an internal object which
+       computes the complex value using the above Eval method and
+       returns its imaginary part. Custom implementations may choose to
+       override this method with a more efficient real-valued
+       coefficient. */
+   virtual VectorCoefficient & imag() { return imag_vcoef_; }
+
+   virtual ~ComplexVectorCoefficient() { }
+};
+
+/** @brief Base class ComplexMatrixCoefficients that optionally depend
+    on space and time. These are used by the SesquilinearForm,
+    ComplexLinearForm, and ComplexGridFunction classes to represent
+    the physical matrix-valued coefficients in the PDEs that are being
+    discretized. This class can also be used in a more general way to
+    represent functions that don't necessarily belong to a FE space.
+    See, e.g., ex22 for these uses. */
+class ComplexMatrixCoefficient
+{
+protected:
+   int height, width;
+   real_t time;
+
+private:
+   RealPartMatrixCoefficient re_part_mcoef_;
+   ImagPartMatrixCoefficient im_part_mcoef_;
+
+protected:
+   MatrixCoefficient &real_mcoef_;
+   MatrixCoefficient &imag_mcoef_;
+
+   mutable DenseMatrix M_r_;
+   mutable DenseMatrix M_i_;
+
+public:
+   /// Construct a dim x dim matrix coefficient.
+   explicit ComplexMatrixCoefficient(int dim)
+      : height(dim), width(dim), time(0.),
+        re_part_mcoef_(*this), im_part_mcoef_(*this),
+        real_mcoef_(re_part_mcoef_), imag_mcoef_(im_part_mcoef_)
+   { }
+
+   /// Construct a h x w matrix coefficient.
+   ComplexMatrixCoefficient(int h, int w) :
+      height(h), width(w), time(0.),
+      re_part_mcoef_(*this), im_part_mcoef_(*this),
+      real_mcoef_(re_part_mcoef_), imag_mcoef_(im_part_mcoef_)
+   { }
+
+   /// Set the time for time dependent coefficients
+   virtual void SetTime(real_t t) { time = t; }
+
+   /// Get the time for time dependent coefficients
+   real_t GetTime() { return time; }
+
+   /// Get the height of the matrix.
+   int GetHeight() const { return height; }
+
+   /// Get the width of the matrix.
+   int GetWidth() const { return width; }
+
+   /// For backward compatibility get the width of the matrix.
+   int GetVDim() const { return width; }
+
+   /** @brief Evaluate the matrix coefficient in the element described by @a T
+       at the point @a ip, storing the result in @a K. */
+   /** @note When this method is called, the caller must make sure that the
+       IntegrationPoint associated with @a T is the same as @a ip. This can be
+       achieved by calling T.SetIntPoint(&ip). */
+   virtual void Eval(ComplexTypeDenseMatrix &K, ElementTransformation &T,
+                     const IntegrationPoint &ip) = 0;
+
+   /** @brief Access a standard Coefficient object reproducing the real part of
+       the complex-valued field */
+   /** @note By default this method returns an internal object which
+       computes the complex value using the above Eval method and
+       returns its real part. Custom implementations may choose to
+       override this method with a more efficient real-valued
+       coefficient. */
+   virtual MatrixCoefficient & real() { return real_mcoef_; }
+
+   /** @brief Access a standard Coefficient object reproducing the imaginary
+       part of the complex-valued field */
+   /** @note By default this method returns an internal object which
+       computes the complex value using the above Eval method and
+       returns its imaginary part. Custom implementations may choose to
+       override this method with a more efficient real-valued
+       coefficient. */
+   virtual MatrixCoefficient & imag() { return imag_mcoef_; }
+
+   virtual ~ComplexMatrixCoefficient() { }
+};
+
+/// A complex-valued coefficient that is constant across space and time
+class ComplexConstantCoefficient : public ComplexCoefficient
+{
+private:
+   complex_t val;
+
+   ConstantCoefficient real_coef;
+   ConstantCoefficient imag_coef;
+
+public:
+   ComplexConstantCoefficient(const complex_t z);
+
+   ComplexConstantCoefficient(real_t z_r, real_t z_i = 0.);
+
+   complex_t Eval(ElementTransformation &T,
+                  const IntegrationPoint &ip) { return val; }
+};
+
+/// Complex-valued vector coefficient that is constant in space and time.
+class ComplexVectorConstantCoefficient : public ComplexVectorCoefficient
+{
+private:
+   ComplexVector vec;
+
+public:
+   /// Construct the coefficient with constant vector @a v.
+   ComplexVectorConstantCoefficient(const ComplexVector &v)
+      : ComplexVectorCoefficient(v.Size()), vec(v) { }
+
+   /// Construct the coefficient with constant vector @a v.
+   ComplexVectorConstantCoefficient(const Vector &v)
+      : ComplexVectorCoefficient(v.Size()), vec(v) { }
+
+   using ComplexVectorCoefficient::Eval;
+
+   ///  Evaluate the vector coefficient at @a ip.
+   void Eval(ComplexVector &V, ElementTransformation &T,
+             const IntegrationPoint &ip) override { V = vec; }
+
+   /// Return a reference to the constant vector in this class.
+   const ComplexVector& GetVec() const { return vec; }
+};
+
+/// Complex-valued vector coefficient that is constant in space and time.
+class ComplexMatrixConstantCoefficient : public ComplexMatrixCoefficient
+{
+private:
+   ComplexTypeDenseMatrix mat;
+
+public:
+   /// Construct the coefficient with constant vector @a v.
+   ComplexMatrixConstantCoefficient(const ComplexTypeDenseMatrix &m)
+      : ComplexMatrixCoefficient(m.Height(), m.Width()), mat(m) { }
+
+   /// Construct the coefficient with constant vector @a v.
+   ComplexMatrixConstantCoefficient(const DenseMatrix &m)
+      : ComplexMatrixCoefficient(m.Height(), m.Width()), mat(m) { }
+
+   using ComplexMatrixCoefficient::Eval;
+
+   ///  Evaluate the matrix coefficient at @a ip.
+   void Eval(ComplexTypeDenseMatrix &M, ElementTransformation &T,
+             const IntegrationPoint &ip) override { M = mat; }
+
+   /// Return a reference to the constant matrix in this class.
+   const ComplexTypeDenseMatrix& GetMat() const { return mat; }
+};
+
+/// A general complex-valued function coefficient
+class ComplexFunctionCoefficient : public ComplexCoefficient
+{
+protected:
+   std::function<complex_t(const Vector &)> Function;
+   std::function<complex_t(const Vector &, real_t)> TDFunction;
+
+public:
+   /// Define a time-independent coefficient from a std function
+   /** \param F time-independent std::function */
+   ComplexFunctionCoefficient(std::function<complex_t
+                              (const Vector &)> F)
+      : Function(std::move(F))
+   { }
+
+   /// Define a time-dependent coefficient from a std function
+   /** \param TDF time-dependent function */
+   ComplexFunctionCoefficient(std::function<complex_t
+                              (const Vector &, real_t)> TDF)
+      : TDFunction(std::move(TDF))
+   { }
+
+   /// (DEPRECATED) Define a time-independent coefficient from a C-function
+   /** @deprecated Use the method where the C-function, @a f, uses a const
+       Vector argument instead of Vector. */
+   MFEM_DEPRECATED ComplexFunctionCoefficient(complex_t
+                                              (*f)(Vector &))
+   {
+      // Cast first to (void*) to suppress a warning from newer version of
+      // Clang when using -Wextra.
+      Function = reinterpret_cast<complex_t(*)
+                 (const Vector&)>((void*)f);
+      TDFunction = NULL;
+   }
+
+   /// (DEPRECATED) Define a time-dependent coefficient from a C-function
+   /** @deprecated Use the method where the C-function, @a tdf, uses a const
+       Vector argument instead of Vector. */
+   MFEM_DEPRECATED ComplexFunctionCoefficient(complex_t
+                                              (*tdf)(Vector &, real_t))
+   {
+      Function = NULL;
+      // Cast first to (void*) to suppress a warning from newer version of
+      // Clang when using -Wextra.
+      TDFunction =
+         reinterpret_cast<complex_t(*)(const Vector&,
+                                       real_t)>((void*)tdf);
+   }
+
+   /// Evaluate the coefficient at @a ip.
+   complex_t Eval(ElementTransformation &T,
+                  const IntegrationPoint &ip) override;
+};
+
+/// A general vector function coefficient
+class ComplexVectorFunctionCoefficient : public ComplexVectorCoefficient
+{
+private:
+   std::function<void(const Vector &, ComplexVector &)> Function;
+   std::function<void(const Vector &, real_t, ComplexVector &)> TDFunction;
+   ComplexCoefficient *Q;
+
+public:
+   /// Define a time-independent complex-valued vector coefficient
+   /// from a std function
+   /** \param dim - the size of the vector
+       \param F - time-independent function
+       \param q - optional scalar Coefficient to scale the vector coefficient */
+   ComplexVectorFunctionCoefficient(int dim,
+                                    std::function<void(const Vector &,
+                                                       ComplexVector &)> F,
+                                    ComplexCoefficient *q = nullptr)
+      : ComplexVectorCoefficient(dim), Function(std::move(F)), Q(q)
+   { }
+
+   /// Define a time-dependent complex-valued vector coefficient from
+   /// a std function
+   /** \param dim - the size of the vector
+       \param TDF - time-dependent function
+       \param q - optional scalar ComplexCoefficient to scale the vector coefficient */
+   ComplexVectorFunctionCoefficient(int dim,
+                                    std::function<void(const Vector &, real_t,
+                                                       ComplexVector &)> TDF,
+                                    ComplexCoefficient *q = nullptr)
+      : ComplexVectorCoefficient(dim), TDFunction(std::move(TDF)), Q(q)
+   { }
+
+   using ComplexVectorCoefficient::Eval;
+   /// Evaluate the vector coefficient at @a ip.
+   void Eval(ComplexVector &V, ElementTransformation &T,
+             const IntegrationPoint &ip) override;
+
+   virtual ~ComplexVectorFunctionCoefficient() { }
+};
+
+} // end namespace mfem
+
+#endif

fem/complex_fem.cpp

@@ -96,6 +96,23 @@ ComplexGridFunction::ProjectCoefficient(Coefficient &real_coeff,
    gfi->SyncAliasMemory(*this);
 }
 
+void
+ComplexGridFunction::ProjectCoefficient(Coefficient &real_coeff)
+{
+   gfr->SyncMemory(*this);
+   gfi->SyncMemory(*this);
+   gfr->ProjectCoefficient(real_coeff);
+   *gfi = 0.0;
+   gfr->SyncAliasMemory(*this);
+   gfi->SyncAliasMemory(*this);
+}
+
+void
+ComplexGridFunction::ProjectCoefficient(ComplexCoefficient &coeff)
+{
+   this->ProjectCoefficient(coeff.real(), coeff.imag());
+}
+
 void
 ComplexGridFunction::ProjectCoefficient(VectorCoefficient &real_vcoeff,
                                         VectorCoefficient &imag_vcoeff)
@@ -108,6 +125,23 @@ ComplexGridFunction::ProjectCoefficient(VectorCoefficient &real_vcoeff,
    gfi->SyncAliasMemory(*this);
 }
 
+void
+ComplexGridFunction::ProjectCoefficient(VectorCoefficient &real_vcoeff)
+{
+   gfr->SyncMemory(*this);
+   gfi->SyncMemory(*this);
+   gfr->ProjectCoefficient(real_vcoeff);
+   *gfi = 0.0;
+   gfr->SyncAliasMemory(*this);
+   gfi->SyncAliasMemory(*this);
+}
+
+void
+ComplexGridFunction::ProjectCoefficient(ComplexVectorCoefficient &vcoeff)
+{
+   this->ProjectCoefficient(vcoeff.real(), vcoeff.imag());
+}
+
 void
 ComplexGridFunction::ProjectBdrCoefficient(Coefficient &real_coeff,
                                            Coefficient &imag_coeff,
@@ -121,6 +155,26 @@ ComplexGridFunction::ProjectBdrCoefficient(Coefficient &real_coeff,
    gfi->SyncAliasMemory(*this);
 }
 
+void
+ComplexGridFunction::ProjectBdrCoefficient(Coefficient &real_coeff,
+                                           Array<int> &attr)
+{
+   ConstantCoefficient zero_coeff(0.0);
+   gfr->SyncMemory(*this);
+   gfi->SyncMemory(*this);
+   gfr->ProjectBdrCoefficient(real_coeff, attr);
+   gfi->ProjectBdrCoefficient(zero_coeff, attr);
+   gfr->SyncAliasMemory(*this);
+   gfi->SyncAliasMemory(*this);
+}
+
+void
+ComplexGridFunction::ProjectBdrCoefficient(ComplexCoefficient &coeff,
+                                           Array<int> &attr)
+{
+   this->ProjectBdrCoefficient(coeff.real(), coeff.imag(), attr);
+}
+
 void
 ComplexGridFunction::ProjectBdrCoefficientNormal(VectorCoefficient &real_vcoeff,
                                                  VectorCoefficient &imag_vcoeff,
@@ -134,6 +188,28 @@ ComplexGridFunction::ProjectBdrCoefficientNormal(VectorCoefficient &real_vcoeff,
    gfi->SyncAliasMemory(*this);
 }
 
+void
+ComplexGridFunction::ProjectBdrCoefficientNormal(VectorCoefficient &real_vcoeff,
+                                                 Array<int> &attr)
+{
+   Vector zero_vec(real_vcoeff.GetVDim()); zero_vec = 0.;
+   VectorConstantCoefficient zero_vcoeff(zero_vec);
+   gfr->SyncMemory(*this);
+   gfi->SyncMemory(*this);
+   gfr->ProjectBdrCoefficientNormal(real_vcoeff, attr);
+   gfi->ProjectBdrCoefficientNormal(zero_vcoeff, attr);
+   gfr->SyncAliasMemory(*this);
+   gfi->SyncAliasMemory(*this);
+}
+
+void
+ComplexGridFunction::ProjectBdrCoefficientNormal(
+   ComplexVectorCoefficient &vcoeff,
+   Array<int> &attr)
+{
+   this->ProjectBdrCoefficientNormal(vcoeff.real(), vcoeff.imag(), attr);
+}
+
 void
 ComplexGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient
                                                   &real_vcoeff,
@@ -149,6 +225,80 @@ ComplexGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient
    gfi->SyncAliasMemory(*this);
 }
 
+void
+ComplexGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient
+                                                  &real_vcoeff,
+                                                  Array<int> &attr)
+{
+   Vector zero_vec(real_vcoeff.GetVDim()); zero_vec = 0.;
+   VectorConstantCoefficient zero_vcoeff(zero_vec);
+   gfr->SyncMemory(*this);
+   gfi->SyncMemory(*this);
+   gfr->ProjectBdrCoefficientTangent(real_vcoeff, attr);
+   gfi->ProjectBdrCoefficientTangent(zero_vcoeff, attr);
+   gfr->SyncAliasMemory(*this);
+   gfi->SyncAliasMemory(*this);
+}
+
+void
+ComplexGridFunction::ProjectBdrCoefficientTangent(
+   ComplexVectorCoefficient &vcoeff,
+   Array<int> &attr)
+{
+   this->ProjectBdrCoefficientTangent(vcoeff.real(), vcoeff.imag(), attr);
+}
+
+real_t
+ComplexGridFunction::ComputeL2Error(Coefficient &re_exsol,
+                                    Coefficient &im_exsol,
+                                    const IntegrationRule *irs[],
+                                    const Array<int> *elems) const
+{
+   real_t err_r = gfr->ComputeL2Error(re_exsol, irs, elems);
+   real_t err_i = gfi->ComputeL2Error(im_exsol, irs, elems);
+
+   return sqrt(err_r * err_r + err_i * err_i);
+}
+
+real_t
+ComplexGridFunction::ComputeL2Error(Coefficient &re_exsol,
+                                    const IntegrationRule *irs[],
+                                    const Array<int> *elems) const
+{
+   ConstantCoefficient zero_coef(0.0);
+
+   real_t err_r = gfr->ComputeL2Error(re_exsol, irs, elems);
+   real_t err_i = gfi->ComputeL2Error(zero_coef, irs, elems);
+
+   return sqrt(err_r * err_r + err_i * err_i);
+}
+
+real_t
+ComplexGridFunction::ComputeL2Error(VectorCoefficient &re_exsol,
+                                    VectorCoefficient &im_exsol,
+                                    const IntegrationRule *irs[],
+                                    const Array<int> *elems) const
+{
+   real_t err_r = gfr->ComputeL2Error(re_exsol, irs, elems);
+   real_t err_i = gfi->ComputeL2Error(im_exsol, irs, elems);
+
+   return sqrt(err_r * err_r + err_i * err_i);
+}
+
+real_t
+ComplexGridFunction::ComputeL2Error(VectorCoefficient &re_exsol,
+                                    const IntegrationRule *irs[],
+                                    const Array<int> *elems) const
+{
+   Vector zero_vec(re_exsol.GetVDim()); zero_vec = 0.0;
+   VectorConstantCoefficient zero_coef(zero_vec);
+
+   real_t err_r = gfr->ComputeL2Error(re_exsol, irs, elems);
+   real_t err_i = gfi->ComputeL2Error(zero_coef, irs, elems);
+
+   return sqrt(err_r * err_r + err_i * err_i);
+}
+
 
 ComplexLinearForm::ComplexLinearForm(FiniteElementSpace *fes,
                                      ComplexOperator::Convention convention)
@@ -731,6 +881,17 @@ ParComplexGridFunction::ProjectCoefficient(Coefficient &real_coeff,
    pgfi->SyncAliasMemory(*this);
 }
 
+void
+ParComplexGridFunction::ProjectCoefficient(Coefficient &real_coeff)
+{
+   pgfr->SyncMemory(*this);
+   pgfi->SyncMemory(*this);
+   pgfr->ProjectCoefficient(real_coeff);
+   *pgfi = 0.0;
+   pgfr->SyncAliasMemory(*this);
+   pgfi->SyncAliasMemory(*this);
+}
+
 void
 ParComplexGridFunction::ProjectCoefficient(VectorCoefficient &real_vcoeff,
                                            VectorCoefficient &imag_vcoeff)
@@ -743,6 +904,17 @@ ParComplexGridFunction::ProjectCoefficient(VectorCoefficient &real_vcoeff,
    pgfi->SyncAliasMemory(*this);
 }
 
+void
+ParComplexGridFunction::ProjectCoefficient(VectorCoefficient &real_vcoeff)
+{
+   pgfr->SyncMemory(*this);
+   pgfi->SyncMemory(*this);
+   pgfr->ProjectCoefficient(real_vcoeff);
+   *pgfi = 0.0;
+   pgfr->SyncAliasMemory(*this);
+   pgfi->SyncAliasMemory(*this);
+}
+
 void
 ParComplexGridFunction::ProjectBdrCoefficient(Coefficient &real_coeff,
                                               Coefficient &imag_coeff,
@@ -756,6 +928,19 @@ ParComplexGridFunction::ProjectBdrCoefficient(Coefficient &real_coeff,
    pgfi->SyncAliasMemory(*this);
 }
 
+void
+ParComplexGridFunction::ProjectBdrCoefficient(Coefficient &real_coeff,
+                                              Array<int> &attr)
+{
+   ConstantCoefficient zero_coeff(0.0);
+   pgfr->SyncMemory(*this);
+   pgfi->SyncMemory(*this);
+   pgfr->ProjectBdrCoefficient(real_coeff, attr);
+   pgfi->ProjectBdrCoefficient(zero_coeff, attr);
+   pgfr->SyncAliasMemory(*this);
+   pgfi->SyncAliasMemory(*this);
+}
+
 void
 ParComplexGridFunction::ProjectBdrCoefficientNormal(VectorCoefficient
                                                     &real_vcoeff,
@@ -771,6 +956,21 @@ ParComplexGridFunction::ProjectBdrCoefficientNormal(VectorCoefficient
    pgfi->SyncAliasMemory(*this);
 }
 
+void
+ParComplexGridFunction::ProjectBdrCoefficientNormal(VectorCoefficient
+                                                    &real_vcoeff,
+                                                    Array<int> &attr)
+{
+   Vector zero_vec(real_vcoeff.GetVDim()); zero_vec = 0.;
+   VectorConstantCoefficient zero_vcoeff(zero_vec);
+   pgfr->SyncMemory(*this);
+   pgfi->SyncMemory(*this);
+   pgfr->ProjectBdrCoefficientNormal(real_vcoeff, attr);
+   pgfi->ProjectBdrCoefficientNormal(zero_vcoeff, attr);
+   pgfr->SyncAliasMemory(*this);
+   pgfi->SyncAliasMemory(*this);
+}
+
 void
 ParComplexGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient
                                                      &real_vcoeff,
@@ -786,6 +986,21 @@ ParComplexGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient
    pgfi->SyncAliasMemory(*this);
 }
 
+void
+ParComplexGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient
+                                                     &real_vcoeff,
+                                                     Array<int> &attr)
+{
+   Vector zero_vec(real_vcoeff.GetVDim()); zero_vec = 0.;
+   VectorConstantCoefficient zero_vcoeff(zero_vec);
+   pgfr->SyncMemory(*this);
+   pgfi->SyncMemory(*this);
+   pgfr->ProjectBdrCoefficientTangent(real_vcoeff, attr);
+   pgfi->ProjectBdrCoefficientTangent(zero_vcoeff, attr);
+   pgfr->SyncAliasMemory(*this);
+   pgfi->SyncAliasMemory(*this);
+}
+
 void
 ParComplexGridFunction::Distribute(const Vector *tv)
 {
@@ -825,6 +1040,31 @@ ParComplexGridFunction::ParallelProject(Vector &tv) const
    tvi.SyncAliasMemory(tv);
 }
 
+real_t
+ParComplexGridFunction::ComputeL2Error(Coefficient &exsolr,
+                                       const IntegrationRule *irs[],
+                                       Array<int> *elems) const
+{
+   ConstantCoefficient zeroCoef(0.0);
+
+   real_t err_r = pgfr->ComputeL2Error(exsolr, irs, elems);
+   real_t err_i = pgfi->ComputeL2Error(zeroCoef, irs, elems);
+   return sqrt(err_r * err_r + err_i * err_i);
+}
+
+real_t
+ParComplexGridFunction::ComputeL2Error(VectorCoefficient &exsolr,
+                                       const IntegrationRule *irs[],
+                                       Array<int> *elems) const
+{
+   Vector zeroVec(exsolr.GetVDim()); zeroVec = 0.0;
+   VectorConstantCoefficient zeroCoef(zeroVec);
+
+   real_t err_r = pgfr->ComputeL2Error(exsolr, irs, elems);
+   real_t err_i = pgfi->ComputeL2Error(zeroCoef, irs, elems);
+   return sqrt(err_r * err_r + err_i * err_i);
+}
+
 
 ParComplexLinearForm::ParComplexLinearForm(ParFiniteElementSpace *pfes,
                                            ComplexOperator::Convention

general/CMakeLists.txt

@@ -39,6 +39,7 @@ list(APPEND HDRS
   arrays_by_name.hpp
   backends.hpp
   binaryio.hpp
+  complex_type.hpp
   cuda.hpp
   device.hpp
   error.hpp

general/complex_type.hpp

@@ -0,0 +1,125 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_COMPLEX_TYPE
+#define MFEM_COMPLEX_TYPE
+
+#include "../config/config.hpp"
+
+#if !(defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP))
+#include <complex>
+#include <utility>
+#endif
+
+#if defined(MFEM_USE_CUDA)
+#include <cuComplex.h>
+#endif
+
+#if defined(MFEM_USE_HIP)
+#include <hip/hip_complex.h>
+#endif
+
+
+namespace mfem
+{
+
+/// @brief Complex number type for device.
+#if !(defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP))
+
+#define zAbs std::abs
+#define zExp std::exp
+#define zNorm std::norm
+using complex_t = std::complex<real_t>;
+
+#else // CUDA or HIP
+
+#if defined(MFEM_USE_CUDA)
+using DoubleComplex_t = cuDoubleComplex;
+#endif
+
+#if defined(MFEM_USE_HIP)
+using DoubleComplex_t = hipDoubleComplex;
+#endif
+
+struct Complex : public DoubleComplex_t
+{
+   MFEM_HOST_DEVICE Complex() = default;
+   MFEM_HOST_DEVICE Complex(real_t r) { x = r, y = 0.0; }
+   MFEM_HOST_DEVICE Complex(real_t r, real_t i) { x = r, y = i; }
+   MFEM_HOST_DEVICE real_t real() const { return x; }
+   MFEM_HOST_DEVICE void real(real_t r) { x = r; }
+   MFEM_HOST_DEVICE real_t imag() const { return y; }
+   MFEM_HOST_DEVICE void imag(real_t i) { y = i; }
+
+   template <typename U>
+   MFEM_HOST_DEVICE inline Complex &operator*=(const U &z)
+   {
+      return *this = *this * z, *this;
+   }
+
+   template <typename U>
+   MFEM_HOST_DEVICE inline Complex &operator/=(const U &z)
+   {
+      return *this = *this / z, *this;
+   }
+};
+
+MFEM_HOST_DEVICE inline Complex operator*(const Complex &x, const real_t &y)
+{
+   return Complex(x.real() * y, x.imag() * y);
+}
+
+MFEM_HOST_DEVICE inline Complex operator+(const Complex &a, const Complex &b)
+{
+   return Complex(a.real() + b.real(), a.imag() + b.imag());
+}
+
+MFEM_HOST_DEVICE inline Complex operator*(const real_t d, const Complex &z)
+{
+   return Complex(z.real() * d, z.imag() * d);
+}
+
+MFEM_HOST_DEVICE inline Complex operator*(const Complex &a, const Complex &b)
+{
+   return Complex(a.real() * b.real() - a.imag() * b.imag(),
+                  a.real() * b.imag() + a.imag() * b.real());
+}
+
+MFEM_HOST_DEVICE inline Complex operator/(const Complex &z, const real_t &d)
+{
+   return Complex(z.real() / d, z.imag() / d);
+}
+
+MFEM_HOST_DEVICE inline real_t zAbs(const Complex &z)
+{
+   return std::hypot(z.real(), z.imag());
+}
+
+MFEM_HOST_DEVICE inline Complex zExp(const Complex &q)
+{
+   Complex z;
+   real_t s, c, e = std::exp(q.real());
+   sincos(q.imag(), &s, &c);
+   z.real(c * e), z.imag(s * e);
+   return z;
+}
+
+MFEM_HOST_DEVICE inline real_t zNorm(const Complex &z)
+{
+   return z.real() * z.real() + z.imag() * z.imag();
+}
+
+using complex_t = Complex;
+#endif // MFEM_USE_CUDA || MFEM_USE_HIP
+
+} // namespace mfem
+
+#endif // MFEM_COMPLEX_TYPE

linalg/CMakeLists.txt

@@ -21,6 +21,7 @@ list(APPEND SRCS
   blockvector.cpp
   complex_densemat.cpp
   complex_operator.cpp
+  complex_vector.cpp
   constraints.cpp
   densemat.cpp
   symmat.cpp
@@ -47,6 +48,7 @@ list(APPEND HDRS
   blockvector.hpp
   complex_densemat.hpp
   complex_operator.hpp
+  complex_vector.hpp
   constraints.hpp
   densemat.hpp
   dinvariants.hpp

linalg/complex_densemat.cpp

@@ -9,6 +9,7 @@
 // terms of the BSD-3 license. We welcome feedback and contributions, see file
 // CONTRIBUTING.md for details.
 
+#include "../general/forall.hpp"
 #include "complex_densemat.hpp"
 #include "lapack.hpp"
 #include <complex>
@@ -16,6 +17,8 @@
 namespace mfem
 {
 
+using namespace std;
+
 DenseMatrix & ComplexDenseMatrix::real()
 {
    MFEM_ASSERT(Op_Real_, "ComplexDenseMatrix has no real part!");
@@ -1017,4 +1020,303 @@ void ComplexCholeskyFactors::GetInverseMatrix(int m, real_t * X_r,
    delete [] X;
 }
 
+ComplexTypeDenseMatrix::ComplexTypeDenseMatrix()
+   : height(0), width(0)
+{}
+
+ComplexTypeDenseMatrix::ComplexTypeDenseMatrix(const ComplexTypeDenseMatrix &m)
+   : height(m.Height()), width(m.Width())
+{
+   const int hw = height * width;
+   if (hw > 0)
+   {
+      MFEM_ASSERT(m.data, "invalid source matrix");
+      data.New(hw);
+      std::memcpy(data, m.data, sizeof(complex_t)*hw);
+   }
+}
+
+ComplexTypeDenseMatrix::ComplexTypeDenseMatrix(const DenseMatrix &m)
+   : height(m.Height()), width(m.Width())
+{
+   const int hw = height * width;
+   if (hw > 0)
+   {
+      MFEM_ASSERT(m.data, "invalid source matrix");
+      data.New(hw);
+
+      for (int i = 0; i < hw; i++)
+      {
+         data[i] = m.data[i];
+      }
+   }
+}
+
+ComplexTypeDenseMatrix::ComplexTypeDenseMatrix(int s)
+   : height(s), width(s)
+{
+   MFEM_ASSERT(s >= 0, "invalid DenseMatrix size: " << s);
+   if (s > 0)
+   {
+      data.New(s*s);
+      *this = 0.0; // init with zeroes
+   }
+}
+
+ComplexTypeDenseMatrix::ComplexTypeDenseMatrix(int m, int n)
+   : height(m), width(n)
+{
+   MFEM_ASSERT(m >= 0 && n >= 0,
+               "invalid DenseMatrix size: " << m << " x " << n);
+   const int capacity = m*n;
+   if (capacity > 0)
+   {
+      data.New(capacity);
+      *this = 0.0; // init with zeroes
+   }
+}
+
+void ComplexTypeDenseMatrix::SetSize(int h, int w)
+{
+   MFEM_ASSERT(h >= 0 && w >= 0,
+               "invalid ComplexTypeDenseMatrix size: " << h << " x " << w);
+   if (Height() == h && Width() == w)
+   {
+      return;
+   }
+   height = h;
+   width = w;
+   const int hw = h*w;
+   if (hw > data.Capacity())
+   {
+      data.Delete();
+      data.New(hw);
+      *this = 0.0; // init with zeroes
+   }
+}
+
+/// Returns reference to a_{ij}.
+complex_t &ComplexTypeDenseMatrix::Elem(int i, int j)
+{
+   return (*this)(i,j);
+}
+
+/// Returns constant reference to a_{ij}.
+const complex_t &ComplexTypeDenseMatrix::Elem(int i, int j) const
+{
+   return (*this)(i,j);
+}
+
+
+ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator=(real_t c)
+{
+   const int s = Height()*Width();
+   for (int i = 0; i < s; i++)
+   {
+      data[i] = c;
+   }
+   return *this;
+}
+
+ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator=(complex_t c)
+{
+   const int s = Height()*Width();
+   for (int i = 0; i < s; i++)
+   {
+      data[i] = c;
+   }
+   return *this;
+}
+
+/// Copy the matrix entries from the given array
+ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator=(const real_t *d)
+{
+   const int s = Height()*Width();
+   for (int i = 0; i < s; i++)
+   {
+      data[i] = d[i];
+   }
+   return *this;
+}
+
+ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator=
+(const complex_t *d)
+{
+   const int s = Height()*Width();
+   for (int i = 0; i < s; i++)
+   {
+      data[i] = d[i];
+   }
+   return *this;
+}
+
+/// Sets the matrix size and elements equal to those of m
+ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator=(const DenseMatrix &m)
+{
+   SetSize(m.height, m.width);
+
+   const int hw = height * width;
+   for (int i = 0; i < hw; i++)
+   {
+      data[i] = m.data[i];
+   }
+
+   return *this;
+}
+
+ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator=
+(const ComplexTypeDenseMatrix &m)
+{
+   SetSize(m.height, m.width);
+
+   const int hw = height * width;
+   for (int i = 0; i < hw; i++)
+   {
+      data[i] = m.data[i];
+   }
+
+   return *this;
+}
+
+ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator+=(const real_t *m)
+{
+   const int s = Height()*Width();
+   for (int i = 0; i < s; i++)
+   {
+      data[i] += m[i];
+   }
+   return *this;
+}
+
+ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator+=
+(const complex_t *m)
+{
+   const int s = Height()*Width();
+   for (int i = 0; i < s; i++)
+   {
+      data[i] += m[i];
+   }
+   return *this;
+}
+
+ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator+=(const DenseMatrix &m)
+{
+   const int hw = height * width;
+   for (int i = 0; i < hw; i++)
+   {
+      data[i] += m.data[i];
+   }
+
+   return *this;
+}
+
+ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator+=
+(const ComplexTypeDenseMatrix &m)
+{
+   const int hw = height * width;
+   for (int i = 0; i < hw; i++)
+   {
+      data[i] += m.data[i];
+   }
+
+   return *this;
+}
+
+ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator-=(const DenseMatrix &m)
+{
+   const int hw = height * width;
+   for (int i = 0; i < hw; i++)
+   {
+      data[i] -= m.data[i];
+   }
+
+   return *this;
+}
+
+ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator-=
+(const ComplexTypeDenseMatrix &m)
+{
+   const int hw = height * width;
+   for (int i = 0; i < hw; i++)
+   {
+      data[i] -= m.data[i];
+   }
+
+   return *this;
+}
+
+ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator*=(real_t c)
+{
+   const int hw = height * width;
+   for (int i = 0; i < hw; i++)
+   {
+      data[i] *= c;
+   }
+
+   return *this;
+}
+
+ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator*=(complex_t c)
+{
+   const int hw = height * width;
+   for (int i = 0; i < hw; i++)
+   {
+      data[i] *= c;
+   }
+
+   return *this;
+}
+
+ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::Set(const DenseMatrix &Mr,
+                                                    const DenseMatrix &Mi)
+{
+   MFEM_ASSERT(height == Mr.Height() && height == Mi.Height() &&
+               width == Mr.Width() && width == Mi.Width(),
+               "incompatible Matrices!");
+
+   const int hw = height * width;
+   for (int i = 0; i < hw; i++)
+   {
+      data[i] = complex_t(Mr.data[i], Mi.data[i]);
+   }
+
+   return *this;
+}
+
+void ComplexTypeDenseMatrix::Swap(ComplexTypeDenseMatrix &other)
+{
+   mfem::Swap(width, other.width);
+   mfem::Swap(height, other.height);
+   mfem::Swap(data, other.data);
+}
+
+ComplexTypeDenseMatrix::~ComplexTypeDenseMatrix()
+{
+   data.Delete();
+}
+
+const DenseMatrix &ComplexTypeDenseMatrix::real() const
+{
+   re_part.SetSize(height, width);
+   const int hw = height * width;
+   for (int i = 0; i < hw; i++)
+   {
+      re_part.data[i] = data[i].real();
+   }
+
+   return re_part;
+}
+
+const DenseMatrix &ComplexTypeDenseMatrix::imag() const
+{
+   im_part.SetSize(height, width);
+   const int hw = height * width;
+   for (int i = 0; i < hw; i++)
+   {
+      im_part.data[i] = data[i].imag();
+   }
+
+   return im_part;
+}
+
 } // mfem namespace

linalg/complex_densemat.hpp

@@ -13,6 +13,7 @@
 #define MFEM_COMPLEX_DENSEMAT
 
 #include "complex_operator.hpp"
+#include "../general/complex_type.hpp"
 #include <complex>
 
 namespace mfem
@@ -241,6 +242,220 @@ public:
 
 };
 
+class ComplexTypeDenseMatrix
+{
+protected:
+   int height; ///< Dimension of the output / number of rows in the matrix.
+   int width;  ///< Dimension of the input / number of columns in the matrix.
+
+private:
+   Memory<complex_t > data;
+
+   mutable DenseMatrix re_part;
+   mutable DenseMatrix im_part;
+
+public:
+   /** Default constructor for DenseMatrix.
+       Sets data = NULL and height = width = 0. */
+   ComplexTypeDenseMatrix();
+
+   /// Copy constructor
+   ComplexTypeDenseMatrix(const ComplexTypeDenseMatrix &);
+   ComplexTypeDenseMatrix(const DenseMatrix &);
+
+   /// Creates square matrix of size s.
+   explicit ComplexTypeDenseMatrix(int s);
+
+   /// Creates rectangular matrix of size m x n.
+   ComplexTypeDenseMatrix(int m, int n);
+
+   /// Construct a ComplexTypeDenseMatrix using an existing data array.
+   /** The ComplexTypeDenseMatrix does not assume ownership of the data array,
+       i.e. it will not delete the array. */
+   ComplexTypeDenseMatrix(complex_t *d, int h, int w)
+      : height(h), width(w) { UseExternalData(d, h, w); }
+
+   /// Create a dense matrix using a braced initializer list
+   /// The inner lists correspond to rows of the matrix
+   template <int M, int N, typename T = real_t>
+   explicit ComplexTypeDenseMatrix(const T (&values)[M][N]) :
+      ComplexTypeDenseMatrix(
+         M, N)
+   {
+      // DenseMatrix is column-major so copies have to be element-wise
+      for (int i = 0; i < M; i++)
+      {
+         for (int j = 0; j < N; j++)
+         {
+            (*this)(i,j) = values[i][j];
+         }
+      }
+   }
+
+   /// Change the data array and the size of the DenseMatrix.
+   /** The DenseMatrix does not assume ownership of the data array, i.e. it will
+       not delete the data array @a d. This method should not be used with
+       DenseMatrix that owns its current data array. */
+   void UseExternalData(complex_t *d, int h, int w)
+   {
+      data.Wrap(d, h*w, false);
+      height = h; width = w;
+   }
+
+   /// Change the data array and the size of the DenseMatrix.
+   /** The DenseMatrix does not assume ownership of the data array, i.e. it will
+       not delete the new array @a d. This method will delete the current data
+       array, if owned. */
+   void Reset(complex_t *d, int h, int w)
+   { if (OwnsData()) { data.Delete(); } UseExternalData(d, h, w); }
+
+   /** Clear the data array and the dimensions of the DenseMatrix. This method
+       should not be used with DenseMatrix that owns its current data array. */
+   void ClearExternalData() { data.Reset(); height = width = 0; }
+
+   /// Delete the matrix data array (if owned) and reset the matrix state.
+   void Clear()
+   { if (OwnsData()) { data.Delete(); } ClearExternalData(); }
+
+   /// Get the height (size of output) of the Operator. Synonym with NumRows().
+   inline int Height() const { return height; }
+   /** @brief Get the number of rows (size of output) of the Operator. Synonym
+       with Height(). */
+   inline int NumRows() const { return height; }
+
+   /// Get the width (size of input) of the Operator. Synonym with NumCols().
+   inline int Width() const { return width; }
+   /** @brief Get the number of columns (size of input) of the Operator. Synonym
+       with Width(). */
+   inline int NumCols() const { return width; }
+
+   /// For backward compatibility define Size to be synonym of Width()
+   int Size() const { return Width(); }
+
+   // Total size = width*height
+   int TotalSize() const { return width*height; }
+
+   /// Change the size of the DenseMatrix to s x s.
+   void SetSize(int s) { SetSize(s, s); }
+
+   /// Change the size of the DenseMatrix to h x w.
+   void SetSize(int h, int w);
+
+   /// Returns the matrix data array.
+   inline complex_t *Data() const
+   {
+      return const_cast<complex_t*>
+             ((const complex_t*)data);
+   }
+
+   /// Returns the matrix data array.
+   inline complex_t *GetData() const { return Data(); }
+
+   Memory<complex_t > &GetMemory() { return data; }
+   const Memory<complex_t > &GetMemory() const { return data; }
+
+   /// Return the DenseMatrix data (host pointer) ownership flag.
+   inline bool OwnsData() const { return data.OwnsHostPtr(); }
+
+   /// Returns reference to a_{ij}.
+   inline complex_t &operator()(int i, int j);
+
+   /// Returns constant reference to a_{ij}.
+   inline const complex_t &operator()(int i, int j) const;
+
+   /// Returns reference to a_{ij}.
+   complex_t &Elem(int i, int j);
+
+   /// Returns constant reference to a_{ij}.
+   const complex_t &Elem(int i, int j) const;
+
+   /// Sets the matrix elements equal to constant c
+   ComplexTypeDenseMatrix &operator=(real_t c);
+   ComplexTypeDenseMatrix &operator=(complex_t c);
+
+   /// Copy the matrix entries from the given array
+   ComplexTypeDenseMatrix &operator=(const real_t *d);
+   ComplexTypeDenseMatrix &operator=(const complex_t *d);
+
+   /// Sets the matrix size and elements equal to those of m
+   ComplexTypeDenseMatrix &operator=(const DenseMatrix &m);
+   ComplexTypeDenseMatrix &operator=(const ComplexTypeDenseMatrix &m);
+
+   ComplexTypeDenseMatrix &operator+=(const real_t *m);
+   ComplexTypeDenseMatrix &operator+=(const complex_t *m);
+   ComplexTypeDenseMatrix &operator+=(const DenseMatrix &m);
+   ComplexTypeDenseMatrix &operator+=(const ComplexTypeDenseMatrix &m);
+
+   ComplexTypeDenseMatrix &operator-=(const DenseMatrix &m);
+   ComplexTypeDenseMatrix &operator-=(const ComplexTypeDenseMatrix &m);
+
+   ComplexTypeDenseMatrix &operator*=(real_t c);
+   ComplexTypeDenseMatrix &operator*=(complex_t c);
+
+   /// (*this) = x + i * y
+   ComplexTypeDenseMatrix &Set(const DenseMatrix &x, const DenseMatrix &y);
+
+   std::size_t MemoryUsage() const
+   { return data.Capacity() * sizeof(complex_t); }
+
+   /// Shortcut for mfem::Read( GetMemory(), TotalSize(), on_dev).
+   const complex_t *Read(bool on_dev = true) const
+   { return mfem::Read(data, Height()*Width(), on_dev); }
+
+   /// Shortcut for mfem::Read(GetMemory(), TotalSize(), false).
+   const complex_t *HostRead() const
+   { return mfem::Read(data, Height()*Width(), false); }
+
+   /// Shortcut for mfem::Write(GetMemory(), TotalSize(), on_dev).
+   complex_t *Write(bool on_dev = true)
+   { return mfem::Write(data, Height()*Width(), on_dev); }
+
+   /// Shortcut for mfem::Write(GetMemory(), TotalSize(), false).
+   complex_t *HostWrite()
+   { return mfem::Write(data, Height()*Width(), false); }
+
+   /// Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), on_dev).
+   complex_t *ReadWrite(bool on_dev = true)
+   { return mfem::ReadWrite(data, Height()*Width(), on_dev); }
+
+   /// Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), false).
+   complex_t *HostReadWrite()
+   { return mfem::ReadWrite(data, Height()*Width(), false); }
+
+   void Swap(ComplexTypeDenseMatrix &other);
+
+   /// Return a reference to the real part of this matrix
+   const DenseMatrix &real() const;
+
+   /// Return a reference to the imaginary part of this matrix
+   const DenseMatrix &imag() const;
+
+   /// Destroys dense matrix.
+   virtual ~ComplexTypeDenseMatrix();
+};
+
+/// Specialization of the template function Swap<> for class ComplexTypeDenseMatrix
+template<> inline void Swap<ComplexTypeDenseMatrix>(ComplexTypeDenseMatrix &a,
+                                                    ComplexTypeDenseMatrix &b)
+{
+   a.Swap(b);
+}
+
+// Inline methods
+
+inline complex_t &ComplexTypeDenseMatrix::operator()(int i, int j)
+{
+   MFEM_ASSERT(data && i >= 0 && i < height && j >= 0 && j < width, "");
+   return data[i+j*height];
+}
+
+inline const complex_t &ComplexTypeDenseMatrix::operator()
+(int i, int j) const
+{
+   MFEM_ASSERT(data && i >= 0 && i < height && j >= 0 && j < width, "");
+   return data[i+j*height];
+}
+
 } // namespace mfem
 
 #endif // MFEM_COMPLEX_DENSEMAT

linalg/complex_vector.cpp

@@ -0,0 +1,424 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "../general/forall.hpp"
+#include "../general/reducers.hpp"
+#include "complex_vector.hpp"
+
+using namespace std;
+
+namespace mfem
+{
+
+ComplexVector::ComplexVector(const ComplexVector &v)
+{
+   const int s = v.Size();
+   size = s;
+   if (s > 0)
+   {
+      MFEM_ASSERT(!v.data.Empty(), "invalid source vector");
+      data.New(s, v.data.GetMemoryType());
+      data.CopyFrom(v.data, s);
+   }
+   UseDevice(v.UseDevice());
+}
+
+ComplexVector::ComplexVector(const Vector &v)
+{
+   const int s = v.Size();
+   size = s;
+   if (s > 0)
+   {
+      MFEM_ASSERT(!v.data.Empty(), "invalid source vector");
+      data.New(s, v.data.GetMemoryType());
+      MFEM_FORALL(i, size, data[i] = v.data[i]; );
+   }
+   UseDevice(v.UseDevice());
+}
+
+ComplexVector::ComplexVector(ComplexVector &&v)
+{
+   *this = std::move(v);
+}
+
+complex_t &ComplexVector::Elem(int i)
+{
+   return operator()(i);
+}
+
+const complex_t &ComplexVector::Elem(int i) const
+{
+   return operator()(i);
+}
+
+complex_t ComplexVector::operator*(const complex_t *v) const
+{
+   HostRead();
+   complex_t dot = 0.0;
+#ifdef MFEM_USE_LEGACY_OPENMP
+   #pragma omp parallel for reduction(+:dot)
+#endif
+   for (int i = 0; i < size; i++)
+   {
+      dot += data[i] * v[i];
+   }
+   return dot;
+}
+
+complex_t ComplexVector::operator*(const real_t *v) const
+{
+   HostRead();
+   complex_t dot = 0.0;
+#ifdef MFEM_USE_LEGACY_OPENMP
+   #pragma omp parallel for reduction(+:dot)
+#endif
+   for (int i = 0; i < size; i++)
+   {
+      dot += data[i] * v[i];
+   }
+   return dot;
+}
+
+complex_t ComplexVector::operator*(const ComplexVector &v) const
+{
+   MFEM_ASSERT(size == v.size, "incompatible Vectors!");
+
+   if (size == 0) { return 0.0; }
+
+   const bool use_dev = UseDevice() || v.UseDevice();
+   const auto m_data = Read(use_dev), v_data = v.Read(use_dev);
+
+   // The standard way of computing the dot product is non-deterministic
+   complex_t prod = 0.0;
+   for (int i = 0; i < size; i++)
+   {
+      prod += m_data[i] * v_data[i];
+   }
+   return prod;
+}
+
+complex_t ComplexVector::operator*(const Vector &v) const
+{
+   MFEM_ASSERT(size == v.size, "incompatible Vectors!");
+
+   if (size == 0) { return 0.0; }
+
+   const bool use_dev = UseDevice() || v.UseDevice();
+   const auto m_data = Read(use_dev);
+   const auto v_data = v.Read(use_dev);
+
+   // The standard way of computing the dot product is non-deterministic
+   complex_t prod = 0.0;
+   for (int i = 0; i < size; i++)
+   {
+      prod += m_data[i] * v_data[i];
+   }
+   return prod;
+}
+
+ComplexVector &ComplexVector::operator=(const complex_t *v)
+{
+   HostRead();
+   MFEM_FORALL(i, size, data[i] = v[i]; );
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator=(const real_t *v)
+{
+   HostRead();
+   MFEM_FORALL(i, size, data[i] = v[i]; );
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator=(const ComplexVector &v)
+{
+#if 0
+   SetSize(v.Size(), v.data.GetMemoryType());
+   data.CopyFrom(v.data, v.Size());
+   UseDevice(v.UseDevice());
+#else
+   SetSize(v.Size());
+   const bool vuse = v.UseDevice();
+   const bool use_dev = UseDevice() || vuse;
+   v.UseDevice(use_dev);
+   // keep 'data' where it is, unless 'use_dev' is true
+   if (use_dev) { Write(); }
+   data.CopyFrom(v.data, v.Size());
+   v.UseDevice(vuse);
+#endif
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator=(const Vector &v)
+{
+   SetSize(v.Size());
+   const bool vuse = v.UseDevice();
+   const bool use_dev = UseDevice() || vuse;
+   v.UseDevice(use_dev);
+   // keep 'data' where it is, unless 'use_dev' is true
+   if (use_dev) { Write(); }
+   MFEM_FORALL(i, size, data[i] = v[i]; );
+   v.UseDevice(vuse);
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator=(ComplexVector &&v)
+{
+   v.Swap(*this);
+   if (this != &v) { v.Destroy(); }
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator=(complex_t value)
+{
+   const bool use_dev = UseDevice();
+   const int N = size;
+   auto y = Write(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] = value; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator=(real_t value)
+{
+   const bool use_dev = UseDevice();
+   const int N = size;
+   auto y = Write(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] = value; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator*=(complex_t c)
+{
+   const bool use_dev = UseDevice();
+   const int N = size;
+   auto y = ReadWrite(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] *= c; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator*=(real_t c)
+{
+   const bool use_dev = UseDevice();
+   const int N = size;
+   auto y = ReadWrite(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] *= c; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator*=(const ComplexVector &v)
+{
+   MFEM_ASSERT(size == v.size, "incompatible Vectors!");
+
+   const bool use_dev = UseDevice() || v.UseDevice();
+   const int N = size;
+   const auto x = v.Read(use_dev);
+   auto y = ReadWrite(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] *= x[i]; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator*=(const Vector &v)
+{
+   MFEM_ASSERT(size == v.size, "incompatible Vectors!");
+
+   const bool use_dev = UseDevice() || v.UseDevice();
+   const int N = size;
+   const auto x = v.Read(use_dev);
+   auto y = ReadWrite(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] *= x[i]; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator/=(complex_t c)
+{
+   const bool use_dev = UseDevice();
+   const int N = size;
+   const complex_t m = conj(c) / norm(c);
+   auto y = ReadWrite(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] *= m; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator/=(real_t c)
+{
+   const bool use_dev = UseDevice();
+   const int N = size;
+   const real_t m = 1.0/c;
+   auto y = ReadWrite(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] *= m; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator/=(const ComplexVector &v)
+{
+   MFEM_ASSERT(size == v.size, "incompatible Vectors!");
+
+   const bool use_dev = UseDevice() || v.UseDevice();
+   const int N = size;
+   const auto x = v.Read(use_dev);
+   auto y = ReadWrite(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] /= x[i]; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator/=(const Vector &v)
+{
+   MFEM_ASSERT(size == v.size, "incompatible Vectors!");
+
+   const bool use_dev = UseDevice() || v.UseDevice();
+   const int N = size;
+   const auto x = v.Read(use_dev);
+   auto y = ReadWrite(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] /= x[i]; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator-=(complex_t c)
+{
+   const bool use_dev = UseDevice();
+   const int N = size;
+   auto y = ReadWrite(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] -= c; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator-=(real_t c)
+{
+   const bool use_dev = UseDevice();
+   const int N = size;
+   auto y = ReadWrite(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] -= c; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator-=(const ComplexVector &v)
+{
+   MFEM_ASSERT(size == v.size, "incompatible Vectors!");
+
+   const bool use_dev = UseDevice() || v.UseDevice();
+   const int N = size;
+   const auto x = v.Read(use_dev);
+   auto y = ReadWrite(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] -= x[i]; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator-=(const Vector &v)
+{
+   MFEM_ASSERT(size == v.size, "incompatible Vectors!");
+
+   const bool use_dev = UseDevice() || v.UseDevice();
+   const int N = size;
+   const auto x = v.Read(use_dev);
+   auto y = ReadWrite(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] -= x[i]; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator+=(complex_t c)
+{
+   const bool use_dev = UseDevice();
+   const int N = size;
+   auto y = ReadWrite(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] += c; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator+=(real_t c)
+{
+   const bool use_dev = UseDevice();
+   const int N = size;
+   auto y = ReadWrite(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] += c; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator+=(const ComplexVector &v)
+{
+   MFEM_ASSERT(size == v.size, "incompatible Vectors!");
+
+   const bool use_dev = UseDevice() || v.UseDevice();
+   const int N = size;
+   const auto x = v.Read(use_dev);
+   auto y = ReadWrite(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] += x[i]; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator+=(const Vector &v)
+{
+   MFEM_ASSERT(size == v.size, "incompatible Vectors!");
+
+   const bool use_dev = UseDevice() || v.UseDevice();
+   const int N = size;
+   const auto x = v.Read(use_dev);
+   auto y = ReadWrite(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] += x[i]; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::Set(const Vector &Vr, const Vector &Vi)
+{
+   MFEM_ASSERT(size == Vr.size && size == Vi.size, "incompatible Vectors!");
+
+   const bool use_dev = UseDevice() || Vr.UseDevice() || Vi.UseDevice();
+   const int N = size;
+   const auto x = Vr.Read(use_dev);
+   const auto y = Vi.Read(use_dev);
+   auto z = Write(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { z[i] = complex_t(x[i], y[i]); });
+   return *this;
+}
+
+const Vector &ComplexVector::real() const
+{
+   re_part.SetSize(size);
+   const bool use_dev = UseDevice();
+   const int N = size;
+   const auto z = Read(use_dev);
+   auto x = re_part.Write(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { x[i] = z[i].real(); });
+   return re_part;
+}
+
+const Vector &ComplexVector::imag() const
+{
+   im_part.SetSize(size);
+   const bool use_dev = UseDevice();
+   const int N = size;
+   const auto z = Read(use_dev);
+   auto y = im_part.Write(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] = z[i].imag(); });
+   return im_part;
+}
+
+}

linalg/complex_vector.hpp

@@ -0,0 +1,479 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_COMPLEX_VECTOR
+#define MFEM_COMPLEX_VECTOR
+
+#include "vector.hpp"
+#include "../general/complex_type.hpp"
+
+namespace mfem
+{
+
+class ComplexVector
+{
+private:
+
+   Memory<complex_t > data;
+   int size;
+
+   mutable Vector re_part;
+   mutable Vector im_part;
+
+public:
+
+   /// Default constructor for ComplexVector. Sets size = 0
+   ComplexVector() : size(0) { }
+
+   /// Copy constructor. Allocates a new data array and copies the data.
+   ComplexVector(const ComplexVector &);
+
+   /// Copy constructor. Allocates a new data array and copies the
+   /// data into real part of this vector.
+   ComplexVector(const Vector &);
+
+   /// Move constructor. "Steals" data from its argument.
+   ComplexVector(ComplexVector&& v);
+
+   /// @brief Creates vector of size s.
+   /// @warning Entries are not initialized to zero!
+   explicit ComplexVector(int s);
+
+   /// Creates a vector referencing an array of complex<doubles>,
+   /// owned by someone else.
+   /// The pointer @a data_ can be NULL. The data array can be replaced later
+   /// with SetData().
+   ComplexVector(complex_t *data_, int size_)
+   { data.Wrap(data_, size_, false); size = size_; }
+
+   /// @brief Create a ComplexVector referencing a sub-vector of the
+   //  ComplexVector @a base starting at the given offset, @a
+   //  base_offset, and size @a size_.
+   ComplexVector(ComplexVector &base, int base_offset, int size_)
+      : data(base.data, base_offset, size_), size(size_) { }
+
+   /// Create a ComplexVector of size @a size_ using MemoryType @a mt.
+   ComplexVector(int size_, MemoryType mt)
+      : data(size_, mt), size(size_) { }
+
+   /// @brief Create a ComplexVector of size @a size_ using host
+   /// MemoryType @a h_mt and device MemoryType @a d_mt.
+   ComplexVector(int size_, MemoryType h_mt, MemoryType d_mt)
+      : data(size_, h_mt, d_mt), size(size_) { }
+
+   /// Create a vector from a statically sized C-style array of convertible type
+   template <typename CT, int N>
+   explicit ComplexVector(const CT (&values)[N]) : ComplexVector(N)
+   { std::copy(values, values + N, begin()); }
+
+   /// Create a vector using a braced initializer list
+   template <typename CT, typename std::enable_if<
+                std::is_convertible<CT,complex_t >::value,bool>::type = true>
+   explicit ComplexVector(std::initializer_list<CT> values) : ComplexVector(
+         values.size())
+   { std::copy(values.begin(), values.end(), begin()); }
+
+   /// Enable execution of Vector operations using the mfem::Device.
+   /// The default is to use Backend::CPU (serial execution on each MPI rank),
+   /// regardless of the mfem::Device configuration.
+   ///
+   /// When appropriate, MFEM functions and class methods will enable the use
+   /// of the mfem::Device for their Vector parameters.
+   ///
+   /// Some derived classes, e.g. GridFunction, enable the use of the
+   /// mfem::Device by default.
+   virtual void UseDevice(bool use_dev) const { data.UseDevice(use_dev); }
+
+   /// Return the device flag of the Memory object used by the Vector
+   virtual bool UseDevice() const { return data.UseDevice(); }
+
+   /// @brief Resize the vector to size @a s.
+   /// If the new size is less than or equal to Capacity() then the internal
+   /// data array remains the same. Otherwise, the old array is deleted, if
+   /// owned, and a new array of size @a s is allocated without copying the
+   /// previous content of the ComplexVector.
+   /// @warning In the second case above (new size greater than current one),
+   /// the vector will allocate new data array, even if it did not own the
+   /// original data! Also, new entries are not initialized!
+   void SetSize(int s);
+
+   /// Resize the vector to size @a s using MemoryType @a mt.
+   void SetSize(int s, MemoryType mt);
+
+   /// Resize the vector to size @a s using the MemoryType of @a v.
+   void SetSize(int s, const ComplexVector &v)
+   { SetSize(s, v.GetMemory().GetMemoryType()); }
+
+   /// Resize the vector to size @a s using the MemoryType of @a v.
+   void SetSize(int s, const Vector &v)
+   { SetSize(s, v.GetMemory().GetMemoryType()); }
+
+   /// Set the Vector data.
+   /// @warning This method should be called only when OwnsData() is false.
+   void SetData(complex_t *d)
+   { data.Wrap(d, data.Capacity(), false); }
+
+   /// Set the Vector data and size.
+   /// The Vector does not assume ownership of the new data. The new size is
+   /// also used as the new Capacity().
+   /// @warning This method should be called only when OwnsData() is false.
+   /// @sa NewDataAndSize().
+   void SetDataAndSize(complex_t *d, int s)
+   { data.Wrap(d, s, false); size = s; }
+
+   /// Set the Vector data and size, deleting the old data, if owned.
+   /// The Vector does not assume ownership of the new data. The new size is
+   /// also used as the new Capacity().
+   /// @sa SetDataAndSize().
+   void NewDataAndSize(complex_t *d, int s)
+   {
+      data.Delete();
+      SetDataAndSize(d, s);
+   }
+
+   /// Reset the Vector to use the given external Memory @a mem and size @a s.
+   /// If @a own_mem is false, the Vector will not own any of the pointers of
+   /// @a mem.
+   ///
+   /// Note that when @a own_mem is true, the @a mem object can be destroyed
+   /// immediately by the caller but `mem.Delete()` should NOT be called since
+   /// the Vector object takes ownership of all pointers owned by @a mem.
+   ///
+   /// @sa NewDataAndSize().
+   inline void NewMemoryAndSize(const Memory<complex_t > &mem,
+                                int s, bool own_mem);
+
+   /// Reset the Vector to be a reference to a sub-vector of @a base.
+   inline void MakeRef(ComplexVector &base, int offset, int size);
+
+   /// @brief Reset the Vector to be a reference to a sub-vector of @a base
+   /// without changing its current size.
+   inline void MakeRef(ComplexVector &base, int offset);
+
+   /// Set the Vector data (host pointer) ownership flag.
+   void MakeDataOwner() const { data.SetHostPtrOwner(true); }
+
+   /// Destroy a vector
+   void Destroy();
+
+   /// @brief Delete the device pointer, if owned. If @a copy_to_host is true
+   /// and the data is valid only on device, move it to host before deleting.
+   /// Invalidates the device memory.
+   void DeleteDevice(bool copy_to_host = true)
+   { data.DeleteDevice(copy_to_host); }
+
+   /// Returns the size of the vector.
+   inline int Size() const { return size; }
+
+   /// Return the size of the currently allocated data array.
+   /// It is always true that Capacity() >= Size().
+   inline int Capacity() const { return data.Capacity(); }
+
+   /// Return a pointer to the beginning of the ComplexVector data.
+   /// @warning This method should be used with caution as it gives write access
+   /// to the data of const-qualified ComplexVector%s.
+   inline complex_t *GetData() const
+   { return const_cast<complex_t*>((const complex_t*)data); }
+
+   /// STL-like begin.
+   inline complex_t *begin() { return data; }
+
+   /// STL-like end.
+   inline complex_t *end() { return data + size; }
+
+   /// STL-like begin (const version).
+   inline const complex_t *begin() const { return data; }
+
+   /// STL-like end (const version).
+   inline const complex_t *end() const { return data + size; }
+
+   /// Return a reference to the Memory object used by the Vector.
+   Memory<complex_t > &GetMemory() { return data; }
+
+   /// @brief Return a reference to the Memory object used by the
+   /// ComplexVector, const version.
+   const Memory<complex_t > &GetMemory() const { return data; }
+
+   /// Update the memory location of the vector to match @a v.
+   void SyncMemory(const ComplexVector &v) const
+   { GetMemory().Sync(v.GetMemory()); }
+
+   /// Update the alias memory location of the vector to match @a v.
+   void SyncAliasMemory(const ComplexVector &v) const
+   { GetMemory().SyncAlias(v.GetMemory(),Size()); }
+
+   /// Read the Vector data (host pointer) ownership flag.
+   inline bool OwnsData() const { return data.OwnsHostPtr(); }
+
+   /// Changes the ownership of the data; after the call the Vector is empty
+   inline void StealData(complex_t **p)
+   { *p = data; data.Reset(); size = 0; }
+
+   /// Changes the ownership of the data; after the call the Vector is empty
+   inline complex_t *StealData()
+   { complex_t *p; StealData(&p); return p; }
+
+   /// Access Vector entries. Index i = 0 .. size-1.
+   complex_t &Elem(int i);
+
+   /// Read only access to Vector entries. Index i = 0 .. size-1.
+   const complex_t &Elem(int i) const;
+
+   /// Access Vector entries using () for 0-based indexing.
+   /// @note If MFEM_DEBUG is enabled, bounds checking is performed.
+   inline complex_t &operator()(int i);
+
+   /// Read only access to Vector entries using () for 0-based indexing.
+   /// @note If MFEM_DEBUG is enabled, bounds checking is performed.
+   inline const complex_t &operator()(int i) const;
+
+   /// Access Vector entries using [] for 0-based indexing.
+   /// @note If MFEM_DEBUG is enabled, bounds checking is performed.
+   inline complex_t &operator[](int i) { return (*this)(i); }
+
+   /// Read only access to Vector entries using [] for 0-based indexing.
+   /// @note If MFEM_DEBUG is enabled, bounds checking is performed.
+   inline const complex_t &operator[](int i) const
+   { return (*this)(i); }
+
+   /// Dot product with a `complex<double> *` array.
+   /// @note No complex conjugate is performed
+   complex_t operator*(const complex_t *v) const;
+   complex_t operator*(const real_t *v) const;
+
+   /// Return the inner-product.
+   /// @note No complex conjugate is performed
+   complex_t operator*(const ComplexVector &v) const;
+   complex_t operator*(const Vector &v) const;
+
+   /// Copy Size() entries from @a v.
+   ComplexVector &operator=(const complex_t *v);
+   ComplexVector &operator=(const real_t *v);
+
+   /// Copy assignment.
+   /// @note Defining this method overwrites the implicitly defined copy
+   /// assignment operator.
+   ComplexVector &operator=(const ComplexVector &v);
+   ComplexVector &operator=(const Vector &v);
+
+   /// Move assignment
+   ComplexVector &operator=(ComplexVector&& v);
+
+   /// Redefine '=' for vector = constant.
+   ComplexVector &operator=(complex_t value);
+   ComplexVector &operator=(real_t value);
+
+   /// Scale vector by a constant
+   ComplexVector &operator*=(complex_t c);
+   ComplexVector &operator*=(real_t c);
+
+   /// Component-wise scaling: (*this)(i) *= v(i)
+   ComplexVector &operator*=(const ComplexVector &v);
+   ComplexVector &operator*=(const Vector &v);
+
+   /// Divide vector by a consant
+   ComplexVector &operator/=(complex_t c);
+   ComplexVector &operator/=(real_t c);
+
+   /// Component-wise division: (*this)(i) /= v(i)
+   ComplexVector &operator/=(const ComplexVector &v);
+   ComplexVector &operator/=(const Vector &v);
+
+   /// Subtract a constant from this vector
+   ComplexVector &operator-=(complex_t c);
+   ComplexVector &operator-=(real_t c);
+
+   /// Subtract a vector from this vector
+   ComplexVector &operator-=(const ComplexVector &v);
+   ComplexVector &operator-=(const Vector &v);
+
+   /// Add a constant to this vector
+   ComplexVector &operator+=(complex_t c);
+   ComplexVector &operator+=(real_t c);
+
+   /// Add a vector to this vector
+   ComplexVector &operator+=(const ComplexVector &v);
+   ComplexVector &operator+=(const Vector &v);
+
+   /// (*this) = x + i * y
+   ComplexVector &Set(const Vector &x, const Vector &y);
+
+   /// Swap the contents of two Vectors
+   inline void Swap(ComplexVector &other);
+
+   /// Return a reference to the real part of this vector
+   const Vector &real() const;
+
+   /// Return a reference to the imaginary part of this vector
+   const Vector &imag() const;
+
+   /// Destroys vector.
+   virtual ~ComplexVector();
+
+   /// Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
+   virtual const complex_t *Read(bool on_dev = true) const
+   { return mfem::Read(data, size, on_dev); }
+
+   /// Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
+   virtual const complex_t *HostRead() const
+   { return mfem::Read(data, size, false); }
+
+   /// Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
+   virtual complex_t *Write(bool on_dev = true)
+   { return mfem::Write(data, size, on_dev); }
+
+   /// Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), false).
+   virtual complex_t *HostWrite()
+   { return mfem::Write(data, size, false); }
+
+   /// Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
+   virtual complex_t *ReadWrite(bool on_dev = true)
+   { return mfem::ReadWrite(data, size, on_dev); }
+
+   /// Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
+   virtual complex_t *HostReadWrite()
+   { return mfem::ReadWrite(data, size, false); }
+};
+
+inline ComplexVector::ComplexVector(int s)
+{
+   MFEM_ASSERT(s>=0,"Unexpected negative size.");
+   size = s;
+   if (s > 0)
+   {
+      data.New(s);
+   }
+}
+
+inline void ComplexVector::SetSize(int s)
+{
+   if (s == size)
+   {
+      return;
+   }
+   if (s <= data.Capacity())
+   {
+      size = s;
+      return;
+   }
+   // preserve a valid MemoryType and device flag
+   const MemoryType mt = data.GetMemoryType();
+   const bool use_dev = data.UseDevice();
+   data.Delete();
+   size = s;
+   data.New(s, mt);
+   data.UseDevice(use_dev);
+}
+
+inline void ComplexVector::SetSize(int s, MemoryType mt)
+{
+   if (mt == data.GetMemoryType())
+   {
+      if (s == size)
+      {
+         return;
+      }
+      if (s <= data.Capacity())
+      {
+         size = s;
+         return;
+      }
+   }
+   const bool use_dev = data.UseDevice();
+   data.Delete();
+   if (s > 0)
+   {
+      data.New(s, mt);
+      size = s;
+   }
+   else
+   {
+      data.Reset();
+      size = 0;
+   }
+   data.UseDevice(use_dev);
+}
+
+inline void ComplexVector::NewMemoryAndSize(
+   const Memory<complex_t > &mem,
+   int s,
+   bool own_mem)
+{
+   data.Delete();
+   size = s;
+   if (own_mem)
+   {
+      data = mem;
+   }
+   else
+   {
+      data.MakeAlias(mem, 0, s);
+   }
+}
+
+inline void ComplexVector::MakeRef(ComplexVector &base, int offset, int s)
+{
+   data.Delete();
+   size = s;
+   data.MakeAlias(base.GetMemory(), offset, s);
+}
+
+inline void ComplexVector::MakeRef(ComplexVector &base, int offset)
+{
+   data.Delete();
+   data.MakeAlias(base.GetMemory(), offset, size);
+}
+
+inline void ComplexVector::Destroy()
+{
+   const bool use_dev = data.UseDevice();
+   data.Delete();
+   size = 0;
+   data.Reset();
+   data.UseDevice(use_dev);
+}
+
+inline complex_t &ComplexVector::operator()(int i)
+{
+   MFEM_ASSERT(data && i >= 0 && i < size,
+               "index [" << i << "] is out of range [0," << size << ")");
+
+   return data[i];
+}
+
+inline const complex_t &ComplexVector::operator()(int i) const
+{
+   MFEM_ASSERT(data && i >= 0 && i < size,
+               "index [" << i << "] is out of range [0," << size << ")");
+
+   return data[i];
+}
+
+inline void ComplexVector::Swap(ComplexVector &other)
+{
+   mfem::Swap(data, other.data);
+   mfem::Swap(size, other.size);
+}
+
+/// Specialization of the template function Swap<> for class ComplexVector
+template<> inline void Swap<ComplexVector>(ComplexVector &a, ComplexVector &b)
+{
+   a.Swap(b);
+}
+
+inline ComplexVector::~ComplexVector()
+{
+   data.Delete();
+}
+
+} // namespace mfem
+
+#endif

linalg/densemat.hpp

@@ -24,6 +24,7 @@ class DenseMatrix : public Matrix
 {
    friend class DenseTensor;
    friend class DenseMatrixInverse;
+   friend class ComplexTypeDenseMatrix;
 
 private:
    Memory<real_t> data;

linalg/vector.hpp

@@ -80,6 +80,8 @@ inline real_t rand_real()
 /// Vector data type.
 class Vector
 {
+   friend class ComplexVector;
+
 protected:
 
    Memory<real_t> data;

tests/unit/CMakeLists.txt

@@ -32,6 +32,7 @@ set(UNIT_TESTS_SRCS
   linalg/test_chebyshev.cpp
   linalg/test_complex_dense_matrix.cpp
   linalg/test_complex_operator.cpp
+  linalg/test_complex_vector.cpp
   linalg/test_constrainedsolver.cpp
   linalg/test_direct_solvers.cpp
   linalg/test_hypre_ilu.cpp

tests/unit/linalg/test_complex_vector.cpp

@@ -0,0 +1,145 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+#include <numeric>
+
+using namespace mfem;
+
+TEST_CASE("Complex Vector init-list and C-style array constructors",
+          "[ComplexVector]")
+{
+   std::complex<real_t> ContigData[6] = {std::complex<real_t>(6.0,1.0),
+                                         std::complex<real_t>(5.0,2.0),
+                                         std::complex<real_t>(4.0,3.0),
+                                         std::complex<real_t>(3.0,4.0),
+                                         std::complex<real_t>(2.0,5.0),
+                                         std::complex<real_t>(1.0,6.0)
+                                        };
+   // Point and size constructor
+   ComplexVector a(ContigData, 6);
+   // Braced-list constructor
+   ComplexVector b({std::complex<real_t>(6.0,1.0),
+                    std::complex<real_t>(5.0,2.0),
+                    std::complex<real_t>(4.0,3.0),
+                    std::complex<real_t>(3.0,4.0),
+                    std::complex<real_t>(2.0,5.0),
+                    std::complex<real_t>(1.0,6.0)});
+   // Statically sized C-style array constructor
+   ComplexVector c(ContigData);
+
+   for (int i = 0; i < a.Size(); i++)
+   {
+      REQUIRE(a[i] == b[i]);
+      REQUIRE(a[i] == c[i]);
+   }
+}
+
+TEST_CASE("Complex Vector Move Constructor", "[ComplexVector]")
+{
+   constexpr int N = 6;
+   std::complex<real_t> ContigData[6] = {std::complex<real_t>(6.0,1.0),
+                                         std::complex<real_t>(5.0,2.0),
+                                         std::complex<real_t>(4.0,3.0),
+                                         std::complex<real_t>(3.0,4.0),
+                                         std::complex<real_t>(2.0,5.0),
+                                         std::complex<real_t>(1.0,6.0)
+                                        };
+   ComplexVector a(ContigData, N);
+   ComplexVector b(N);
+   for (int i = 0; i < N; i++)
+   {
+      b(i) = std::complex<real_t>(N - i, i + 1);
+   }
+
+   std::complex<real_t>* a_data = a.GetData();
+   std::complex<real_t>* b_data = b.GetData();
+
+   ComplexVector move_non_owning(std::move(a));
+   ComplexVector move_owning(std::move(b));
+
+   REQUIRE(a.Size() == 0);
+   REQUIRE(a.GetData() == nullptr);
+   REQUIRE(b.Size() == 0);
+   REQUIRE(b.GetData() == nullptr);
+
+   // Should both be no-ops
+   a.Destroy();
+   b.Destroy();
+
+   REQUIRE(move_non_owning.OwnsData() == false);
+   REQUIRE(move_owning.OwnsData() == true);
+
+   REQUIRE(move_non_owning.Size() == N);
+   REQUIRE(move_owning.Size() == N);
+
+   // Make sure that the pointers were reused
+   REQUIRE(move_non_owning.GetData() == a_data);
+   REQUIRE(move_owning.GetData() == b_data);
+
+   for (int i = 0; i < N; i++)
+   {
+      REQUIRE(move_non_owning(i) == std::complex<real_t>(N - i, i + 1));
+      REQUIRE(move_owning(i) == std::complex<real_t>(N - i, i + 1));
+   }
+}
+
+TEST_CASE("Complex Vector Move Assignment", "[ComplexVector]")
+{
+   constexpr int N = 6;
+   std::complex<real_t> ContigData[6] = {std::complex<real_t>(6.0,1.0),
+                                         std::complex<real_t>(5.0,2.0),
+                                         std::complex<real_t>(4.0,3.0),
+                                         std::complex<real_t>(3.0,4.0),
+                                         std::complex<real_t>(2.0,5.0),
+                                         std::complex<real_t>(1.0,6.0)
+                                        };
+   ComplexVector a(ContigData, N);
+   ComplexVector b(N);
+   for (int i = 0; i < N; i++)
+   {
+      b(i) = std::complex<real_t>(N - i, i + 1);
+   }
+
+   std::complex<real_t>* a_data = a.GetData();
+   std::complex<real_t>* b_data = b.GetData();
+
+   ComplexVector move_non_owning;
+   move_non_owning = std::move(a);
+   ComplexVector move_owning;
+   move_owning = std::move(b);
+
+   REQUIRE(a.Size() == 0);
+   REQUIRE(a.GetData() == nullptr);
+   REQUIRE(b.Size() == 0);
+   REQUIRE(b.GetData() == nullptr);
+
+   // Should both be no-ops
+   a.Destroy();
+   b.Destroy();
+
+   REQUIRE(move_non_owning.OwnsData() == false);
+   REQUIRE(move_owning.OwnsData() == true);
+
+   REQUIRE(move_non_owning.Size() == N);
+   REQUIRE(move_owning.Size() == N);
+
+   // Make sure that the pointers were reused
+   REQUIRE(move_non_owning.GetData() == a_data);
+   REQUIRE(move_owning.GetData() == b_data);
+
+   for (int i = 0; i < N; i++)
+   {
+      REQUIRE(move_non_owning(i) == std::complex<real_t>(N - i, i + 1));
+      REQUIRE(move_owning(i) == std::complex<real_t>(N - i, i + 1));
+   }
+}