unit tests from complex dense matrix: Mult, MultAtB, Inverse, SystemMatrix

linalg/CMakeLists.txt

@@ -15,6 +15,7 @@ list(APPEND SRCS
   blockoperator.cpp
   blockvector.cpp
   complex_operator.cpp
+  complex_densemat.cpp
   constraints.cpp
   densemat.cpp
   symmat.cpp
@@ -34,6 +35,7 @@ list(APPEND HDRS
   blockoperator.hpp
   blockvector.hpp
   complex_operator.hpp
+  complex_densemat.hpp
   constraints.hpp
   densemat.hpp
   dinvariants.hpp

linalg/complex_densemat.cpp

@@ -0,0 +1,366 @@
+// Copyright (c) 2010-2021, 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_densemat.hpp"
+#include <complex>
+
+#ifdef MFEM_USE_LAPACK
+extern "C" void
+zgetrf_(int *, int *, std::complex<double> *, int *, int *, int *);
+extern "C" void
+zgetrs_(char *, int *, int *, std::complex<double> *, int *, int *,
+        std::complex<double> *, int *, int *);
+extern "C" void
+zgetri_(int *N, std::complex<double> *A, int *LDA, int *IPIV,
+        std::complex<double> *WORK,
+        int *LWORK, int *INFO);
+#endif
+
+namespace mfem
+{
+
+DenseMatrix & ComplexDenseMatrix::real()
+{
+   MFEM_ASSERT(Op_Real_, "ComplexDenseMatrix has no real part!");
+   return dynamic_cast<DenseMatrix &>(*Op_Real_);
+}
+
+DenseMatrix & ComplexDenseMatrix::imag()
+{
+   MFEM_ASSERT(Op_Imag_, "ComplexDenseMatrix has no imaginary part!");
+   return dynamic_cast<DenseMatrix &>(*Op_Imag_);
+}
+
+const DenseMatrix & ComplexDenseMatrix::real() const
+{
+   MFEM_ASSERT(Op_Real_, "ComplexDenseMatrix has no real part!");
+   return dynamic_cast<const DenseMatrix &>(*Op_Real_);
+}
+
+const DenseMatrix & ComplexDenseMatrix::imag() const
+{
+   MFEM_ASSERT(Op_Imag_, "ComplexDenseMatrix has no imaginary part!");
+   return dynamic_cast<const DenseMatrix &>(*Op_Imag_);
+}
+
+DenseMatrix * ComplexDenseMatrix::GetSystemMatrix() const
+{
+   int h = height/2;
+   int w = width/2;
+   DenseMatrix * A = new DenseMatrix(2*h,2*w);
+   double * data = A->Data();
+   double * data_r = nullptr;
+   double * data_i = nullptr;
+
+   // assuming Hermitian convension
+   *A = 0.;
+   if (hasRealPart())
+   {
+      data_r = real().Data();
+      for (int j = 0; j<w; j++)
+      {
+         for (int i = 0; i<h; i++)
+         {
+            data[i+j*height] = data_r[i+j*h];
+            data[i+h+(j+h)*height] = data_r[i+j*h];
+         }
+      }
+   }
+   if (hasImagPart())
+   {
+      data_i = imag().Data();
+      for (int j = 0; j<w; j++)
+      {
+         for (int i = 0; i<h; i++)
+         {
+            data[i+h+j*height] = data_i[i+j*h];
+            data[i+(j+h)*height] = -data_i[i+j*h];
+         }
+      }
+   }
+   return A;
+}
+
+ComplexDenseMatrix * ComplexDenseMatrix::ComputeInverse()
+{
+   MFEM_VERIFY(height == width, "Matrix has to be square");
+
+   // complex data
+   int h = height/2;
+   int w = width/2;
+   std::complex<double> * data = new std::complex<double>[h*w];
+
+   // copy data
+   if (hasRealPart() && hasImagPart())
+   {
+      double * data_r = real().Data();
+      double * data_i = imag().Data();
+      for (int i = 0; i < h*w; i++)
+      {
+         data[i] = std::complex<double> (data_r[i], data_i[i]);
+      }
+   }
+   else
+   {
+      MFEM_ABORT("ComplexDenseMatrix has either only real or imag part");
+   }
+
+#ifdef MFEM_USE_LAPACK
+   int   *ipiv = new int[w];
+   int    lwork = -1;
+   std::complex<double> qwork, *work;
+   int    info;
+
+   zgetrf_(&w, &w, data, &w, ipiv, &info);
+   if (info)
+   {
+      mfem_error("DenseMatrix::Invert() : Error in ZGETRF");
+   }
+
+   zgetri_(&w, data, &w, ipiv, &qwork, &lwork, &info);
+
+   lwork = (int) qwork.real();
+   work = new std::complex<double>[lwork];
+
+   zgetri_(&w, data, &w, ipiv, work, &lwork, &info);
+
+   if (info)
+   {
+      mfem_error("DenseMatrix::Invert() : Error in ZGETRI");
+   }
+
+   delete [] work;
+   delete [] ipiv;
+
+#else
+   // compiling without LAPACK
+   int c, i, j, n = w;
+   double a, b;
+   Array<int> piv(n);
+   std::complex<double> ac,bc;
+
+   for (c = 0; c < n; c++)
+   {
+      a = std::abs(data[c+c*h]);
+      i = c;
+      for (j = c + 1; j < n; j++)
+      {
+         b = std::abs(data[j+c*h]);
+         if (a < b)
+         {
+            a = b;
+            i = j;
+         }
+      }
+      if (a == 0.0)
+      {
+         mfem_error("DenseMatrix::Invert() : singular matrix");
+      }
+      piv[c] = i;
+      for (j = 0; j < n; j++)
+      {
+         mfem::Swap<std::complex<double>>(data[c+j*h], data[i+j*h]);
+      }
+
+      ac = data[c+c*h] = 1.0 / data[c+c*h];
+      for (j = 0; j < c; j++)
+      {
+         data[c+j*h] *= ac;
+      }
+      for (j++; j < n; j++)
+      {
+         data[c+j*h] *= ac;
+      }
+      for (i = 0; i < c; i++)
+      {
+         data[i+c*h] = ac * (bc = -data[i+c*h]);
+         for (j = 0; j < c; j++)
+         {
+            data[i+j*h] += bc * data[c+j*h];
+         }
+         for (j++; j < n; j++)
+         {
+            data[i+j*h] += bc * data[c+j*h];
+         }
+      }
+      for (i++; i < n; i++)
+      {
+         data[i+c*h] = ac * (bc = -data[i+c*h]);
+         for (j = 0; j < c; j++)
+         {
+            data[i+j*h] += bc * data[c+j*h];
+         }
+         for (j++; j < n; j++)
+         {
+            data[i+j*h] += bc * data[c+j*h];
+         }
+      }
+   }
+
+   for (c = n - 1; c >= 0; c--)
+   {
+      j = piv[c];
+      for (i = 0; i < n; i++)
+      {
+         mfem::Swap<std::complex<double>>(data[i+c*h], data[i+j*h]);
+      }
+   }
+
+#endif
+
+   DenseMatrix * C_r = new DenseMatrix(h);
+   DenseMatrix * C_i = new DenseMatrix(h);
+
+   double * datac_r = C_r->Data();
+   double * datac_i = C_i->Data();
+
+   for (int i = 0; i < h*w; i++)
+   {
+      datac_r[i] = data[i].real();
+      datac_i[i] = data[i].imag();
+   }
+
+   return new ComplexDenseMatrix(C_r,C_i,true,true);
+
+}
+
+ComplexDenseMatrix * Mult(const ComplexDenseMatrix &A,
+                          const ComplexDenseMatrix &B)
+{
+   // C = C_r + i C_i = (A_r + i * A_i) * (B_r + i * B_i)
+   //                 = A_r * B_r - A_i B_i + i (A_r * B_i + A_i * B_r)
+
+   int h = A.Height()/2;
+   int w = B.Width()/2;
+
+   MFEM_VERIFY(A.Width() == B.Height(), "Incompatible matrix dimenions");
+
+   //only real case (imag is null)
+   DenseMatrix * C_r = nullptr;
+   DenseMatrix * C_i = nullptr;
+   if ((A.hasRealPart() && B.hasRealPart()) ||
+       (A.hasImagPart() && B.hasImagPart()))
+   {
+      C_r = new DenseMatrix(h,w);
+   }
+   if ((A.hasRealPart() && B.hasImagPart()) ||
+       (A.hasImagPart() && B.hasRealPart()))
+   {
+      C_i = new DenseMatrix(h,w);
+   }
+
+   MFEM_VERIFY(C_r || C_i, "Both real and imag parts are null");
+
+   if (A.hasRealPart() && B.hasRealPart())
+   {
+      Mult(A.real(), B.real(),*C_r);
+   }
+   if (A.hasImagPart() && B.hasImagPart())
+   {
+      if (A.hasRealPart() && B.hasRealPart())
+      {
+         AddMult_a(-1.,A.imag(), B.imag(),*C_r);
+      }
+      else
+      {
+         Mult(A.imag(), B.imag(),*C_r);
+      }
+   }
+
+   if (A.hasRealPart() && B.hasImagPart())
+   {
+      Mult(A.real(), B.imag(),*C_i);
+   }
+
+   if (A.hasImagPart() && B.hasRealPart())
+   {
+      if (A.hasRealPart() && B.hasImagPart())
+      {
+         AddMult(A.imag(), B.real(),*C_i);
+      }
+      else
+      {
+         Mult(A.imag(), B.real(),*C_i);
+      }
+   }
+
+   return new ComplexDenseMatrix(C_r,C_i,true,true);
+}
+
+ComplexDenseMatrix * MultAtB(const ComplexDenseMatrix &A,
+                             const ComplexDenseMatrix &B)
+{
+   // C = C_r + i C_i = (A_r^t - i * A_i^t) * (B_r + i * B_i)
+   //                     = A_r^t * B_r + A_i^t * B_i + i (A_r^t * B_i - A_i^t * B_r)
+
+   int h = A.Width()/2;
+   int w = B.Width()/2;
+
+   MFEM_VERIFY(A.Height() == B.Height(), "Incompatible matrix dimenions");
+
+   //only real case (imag is null)
+   DenseMatrix * C_r = nullptr;
+   DenseMatrix * C_i = nullptr;
+   if ((A.hasRealPart() && B.hasRealPart()) ||
+       (A.hasImagPart() && B.hasImagPart()))
+   {
+      C_r = new DenseMatrix(h,w);
+   }
+   if ((A.hasRealPart() && B.hasImagPart()) ||
+       (A.hasImagPart() && B.hasRealPart()))
+   {
+      C_i = new DenseMatrix(h,w);
+   }
+
+   MFEM_VERIFY(C_r || C_i, "Both real and imag parts are null");
+
+   if (A.hasRealPart() && B.hasRealPart())
+   {
+      MultAtB(A.real(), B.real(),*C_r);
+   }
+   if (A.hasImagPart() && B.hasImagPart())
+   {
+      if (A.hasRealPart() && B.hasRealPart())
+      {
+         DenseMatrix tempC_r(h,w);
+         MultAtB(A.imag(), B.imag(),tempC_r);
+         (*C_r) += tempC_r;
+      }
+      else
+      {
+         MultAtB(A.imag(), B.imag(),*C_r);
+      }
+   }
+
+   if (A.hasRealPart() && B.hasImagPart())
+   {
+      MultAtB(A.real(), B.imag(),*C_i);
+   }
+
+   if (A.hasImagPart() && B.hasRealPart())
+   {
+      if (A.hasRealPart() && B.hasImagPart())
+      {
+         DenseMatrix tempC_i(h,w);
+         MultAtB(A.imag(), B.real(),tempC_i);
+         (*C_i) -= tempC_i;
+      }
+      else
+      {
+         MultAtB(A.imag(), B.real(),*C_i);
+      }
+   }
+
+   return new ComplexDenseMatrix(C_r,C_i,true,true);
+}
+
+
+} // mfem namespace

linalg/complex_densemat.hpp

@@ -0,0 +1,65 @@
+// Copyright (c) 2010-2021, 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_DENSEMAT
+#define MFEM_COMPLEX_DENSEMAT
+
+#include "complex_operator.hpp"
+
+namespace mfem
+{
+
+/** @brief Specialization of the ComplexOperator built from a pair of Dense
+    Matrices.
+
+    The purpose of this specialization is to support the inverse of a
+    ComplexDenseMatrix and various MatMat operations
+
+    See ComplexOperator documentation for more information.
+    Note: Only the Hermitian convention is supported
+ */
+class ComplexDenseMatrix : public ComplexOperator
+{
+
+public:
+   ComplexDenseMatrix(DenseMatrix * A_Real, DenseMatrix * A_Imag,
+                      bool ownReal, bool ownImag)
+      : ComplexOperator(A_Real, A_Imag, ownReal, ownImag)
+   { }
+
+   virtual DenseMatrix & real();
+   virtual DenseMatrix & imag();
+
+   virtual const DenseMatrix & real() const;
+   virtual const DenseMatrix & imag() const;
+
+   /** Combine the blocks making up this complex operator into a single
+       DenseMatrix. Note that this combined operator requires roughly
+       twice the memory of the block structured operator. */
+   DenseMatrix * GetSystemMatrix() const;
+
+   virtual Type GetType() const { return Complex_DenseMat; }
+
+   ComplexDenseMatrix * ComputeInverse();
+
+};
+
+/// Matrix matrix multiplication.  A = B * C.
+ComplexDenseMatrix * Mult(const ComplexDenseMatrix &B,
+                          const ComplexDenseMatrix &C);
+
+/// Multiply the Complex transpose of a matrix A with a matrix B:   Ah*B
+ComplexDenseMatrix * MultAtB(const ComplexDenseMatrix &A,
+                             const ComplexDenseMatrix &B);
+
+} // namespace mfem
+
+#endif // MFEM_COMPLEX_DENSEMAT

linalg/linalg.hpp

@@ -19,6 +19,7 @@
 #include "matrix.hpp"
 #include "sparsemat.hpp"
 #include "complex_operator.hpp"
+#include "complex_densemat.hpp"
 #include "blockvector.hpp"
 #include "blockmatrix.hpp"
 #include "blockoperator.hpp"

linalg/operator.hpp

@@ -265,7 +265,8 @@ public:
       PETSC_MATGENERIC, ///< ID for class PetscParMatrix, unspecified format.
       Complex_Operator, ///< ID for class ComplexOperator.
       MFEM_ComplexSparseMat, ///< ID for class ComplexSparseMatrix.
-      Complex_Hypre_ParCSR   ///< ID for class ComplexHypreParMatrix.
+      Complex_Hypre_ParCSR,   ///< ID for class ComplexHypreParMatrix.
+      Complex_DenseMat  ///< ID for class ComplexDenseMatrix
    };
 
    /// Return the type ID of the Operator class.

tests/unit/CMakeLists.txt

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

tests/unit/linalg/test_complex_dense_matrix.cpp

@@ -0,0 +1,277 @@
+// Copyright (c) 2010-2021, 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"
+
+using namespace mfem;
+
+TEST_CASE("ComplexDenseMatrix", "[ComplexDenseMatrix]")
+{
+   DenseMatrix A_r(
+   {
+      {
+         7.449631539398011e-01, 5.063269213774142e-01,
+         9.973798911162260e-01, 1.538014630136140e-01
+      },
+      {
+         3.242251927854619e-01, 5.449572274461251e-01,
+         3.627747214847648e-01, 1.751467921062383e-01
+      },
+      {
+         3.021051758964581e-02, 7.644365599277384e-01,
+         6.409585121280975e-01, 2.964283844937954e-01
+      },
+      {
+         8.934345605565150e-01, 5.936906950509614e-01,
+         9.667885787036209e-01, 3.997214406644853e-01
+      }
+   });
+
+   DenseMatrix A_i(
+   {
+      {
+         7.502971220456730e-01, 9.403353505608869e-01,
+         5.710855592871427e-01, 2.532954183189922e-02
+      },
+      {
+         5.439505663874815e-02, 1.994881941837258e-01,
+         7.235874076138660e-01, 2.960357099674006e-01
+      },
+      {
+         1.226216656446337e-02, 6.376122492902638e-01,
+         4.088277326796946e-01, 3.190335595202387e-01
+      },
+      {
+         3.399382051588566e-01, 5.349445908438923e-01,
+         5.141016802691041e-01, 3.124674268283065e-01
+      }
+   });
+
+   DenseMatrix B_r(
+   {
+      {
+         5.153034716511075e-01, 6.753590227871137e-01,
+         6.488329597542797e-01, 5.875278152498098e-01
+      },
+      {
+         4.817945488784101e-01, 1.506780316517289e-01,
+         3.272389525639043e-01, 9.870211717917942e-01
+      },
+      {
+         4.967654008687812e-01, 4.772880018003378e-01,
+         1.669243944018936e-01, 5.566959215458547e-01
+      },
+      {
+         3.884369981782157e-01, 6.518319777336401e-01,
+         4.941412579410056e-01, 6.560086428625783e-01
+      }
+   });
+
+   DenseMatrix B_i(
+   {
+      {
+         1.715019866173268e-01, 9.756201041995241e-01,
+         9.988839596883219e-01, 2.291758588510961e-01
+      },
+      {
+         3.445918467879527e-01, 8.640938244547280e-01,
+         7.006319538228559e-01, 6.955297132297025e-01
+      },
+      {
+         5.683873148693049e-01, 8.407962324993539e-01,
+         5.278065857458160e-01, 7.078222981496557e-01
+      },
+      {
+         5.352425442921088e-01, 7.445666042983262e-01,
+         9.595725554185789e-01, 5.915498198096532e-01
+      }
+   });
+
+   DenseMatrix AB_r(
+   {
+      {
+         3.921690721310207e-01, -8.878689447367504e-01,
+         -8.424837843181997e-01, 3.483824291867016e-01
+      },
+      {
+         3.007953729123389e-02, -4.658575034568339e-01,
+         -3.242816640008720e-01, 2.067214846847279e-01
+      },
+      {
+         1.924669164183410e-01, -5.094725252559890e-01,
+         -4.576737102182742e-01, 3.991530729424347e-01
+      },
+      {
+         6.798682275555494e-01, -4.396738816577034e-02,
+         -1.526715491253031e-01, 9.126227486839693e-01
+      }
+   });
+
+   DenseMatrix AB_i(
+   {
+      {
+         2.084671590130670, 3.054914385318618,
+         2.675264243507663, 3.023330671994920
+      },
+      {
+         1.141922343769045, 1.827762673469408,
+         1.432860194672228, 1.639604444532814
+      },
+      {
+         1.432106293784903, 1.957085447970639,
+         1.631010776832001, 2.241071404311889
+      },
+      {
+         1.930931621959122, 3.254384169173436,
+         2.838072204721421, 2.757357240594122
+      }
+   });
+
+   DenseMatrix AtB_r(
+   {
+      {
+         1.238483198389297e+00, 2.191183543124149e+00,
+         2.156219610054944e+00, 1.780175753880010e+00
+      },
+      {
+         2.012573711609136e+00, 3.200098959139964e+00,
+         2.856730798036862e+00, 2.772405085471168e+00
+      },
+      {
+         2.237504886974010e+00, 3.573290931326598e+00,
+         3.137085199309365e+00, 3.162745119391227e+00
+      },
+      {
+         9.210972068474134e-01, 1.313704573046731e+00,
+         1.105042827289885e+00, 1.312843091304551e+00
+      }
+   });
+
+   DenseMatrix AtB_i(
+   {
+      {
+         1.838895315375063e-01, 9.552331942138000e-01,
+         1.169912839045891e+00,  2.217923570878820e-01
+      },
+      {
+         -7.831791871315461e-02, 7.315134241118445e-01,
+            8.145681272268642e-01, -6.790194154090945e-02
+         },
+      {
+         1.321504844031072e-01, 1.520336262933564e+00,
+         1.586836956019061e+00, -1.080882911419373e-01
+      },
+      {
+         3.362593185692363e-02, 4.305904017992347e-01,
+         4.953952430293495e-01, -8.632019073781308e-02
+      }
+   });
+
+   DenseMatrix invA_r(
+   {
+      {
+         -6.213684345811854e-02,  2.582305494586467e-01,
+            -9.887388608096964e-01,  8.966161168842675e-01
+         },
+      {
+         2.869296542606237e-01,  7.437992704193498e-01,
+         4.018887102703056e-01, -7.468194529297197e-01
+      },
+      {
+         4.609546822804680e-01, -2.488378792722463e-01,
+         2.680431511842616e-01, -5.000096584597614e-01
+      },
+      {
+         -1.104494530992195e+00, -7.039157446455697e-01,
+            2.000991732708908e-01,  2.094621838961788e+00
+         }
+   });
+
+   DenseMatrix invA_i(
+   {
+      {
+         -9.514887849903558e-02,  6.972795844477648e-01,
+            1.368385088611629e-01, -4.694496308988390e-01
+         },
+      {
+         -3.084658924431760e-01,  7.692368880071571e-01,
+            -6.991229727229065e-01, -5.916993505502399e-03
+         },
+      {
+         -5.200646743549728e-01, -1.876793828449070e+00,
+            6.006690970332771e-01,  9.090251794168793e-01
+         },
+      {
+         1.868069474938261e+00,  1.493553266782004e+00,
+         -9.184332408291080e-01, -1.897901970739719e+00
+      }
+   });
+
+   ComplexDenseMatrix A(&A_r,&A_i,false,false);
+   ComplexDenseMatrix B(&B_r,&B_i,false,false);
+
+   SECTION("Mult")
+   {
+      ComplexDenseMatrix * AB = Mult(A,B);
+      AB_r -= AB->real();
+      AB_i -= AB->imag();
+
+      double norm_r = AB_r.MaxMaxNorm();
+      double norm_i = AB_i.MaxMaxNorm();
+
+      REQUIRE(norm_r == MFEM_Approx(0.));
+      REQUIRE(norm_i == MFEM_Approx(0.));
+      delete AB;
+   }
+   SECTION("MultAtB")
+   {
+      ComplexDenseMatrix * AtB = MultAtB(A,B);
+      AtB_r -= AtB->real();
+      AtB_i -= AtB->imag();
+
+      double norm_r = AtB_r.MaxMaxNorm();
+      double norm_i = AtB_i.MaxMaxNorm();
+
+      REQUIRE(norm_r == MFEM_Approx(0.));
+      REQUIRE(norm_i == MFEM_Approx(0.));
+      delete AtB;
+   }
+   SECTION("Inverse")
+   {
+      ComplexDenseMatrix * invA = A.ComputeInverse();
+      invA_r -= invA->real();
+      invA_i -= invA->imag();
+
+      double norm_r = invA_r.MaxMaxNorm();
+      double norm_i = invA_i.MaxMaxNorm();
+
+      REQUIRE(norm_r == MFEM_Approx(0.));
+      REQUIRE(norm_i == MFEM_Approx(0.));
+      delete invA;
+   }
+   SECTION("SystemMatrix")
+   {
+      DenseMatrix * sA = A.GetSystemMatrix();
+      sA->Invert();
+      ComplexDenseMatrix * invA = A.ComputeInverse();
+      DenseMatrix * sinvA = invA->GetSystemMatrix();
+
+      *sA-=*sinvA;
+      double norm = sA->MaxMaxNorm();
+      REQUIRE(norm == MFEM_Approx(0.));
+
+      delete sinvA;
+      delete sA;
+   }
+
+}
+