unit KronMult methods

linalg/CMakeLists.txt

@@ -17,6 +17,7 @@ list(APPEND SRCS
   complex_operator.cpp
   constraints.cpp
   densemat.cpp
+  fdsolver.cpp
   symmat.cpp
   handle.cpp
   matrix.cpp
@@ -39,6 +40,7 @@ list(APPEND HDRS
   dinvariants.hpp
   symmat.hpp
   dtensor.hpp
+  fdsolver.hpp
   handle.hpp
   invariants.hpp
   kernels.hpp

linalg/densemat.cpp

@@ -50,6 +50,10 @@ dsyevr_(char *JOBZ, char *RANGE, char *UPLO, int *N, double *A, int *LDA,
         double *W, double *Z, int *LDZ, int *ISUPPZ, double *WORK, int *LWORK,
         int *IWORK, int *LIWORK, int *INFO);
 extern "C" void
+dgeev_(const char * jobvl, const char * jobvr, int *n, double * A, int * lda,
+       double * wr, double * wl, double * vl, int * ldvl, double * vr, int * ldvr,
+       double * work, int * lwork, int * info);
+extern "C" void
 dsyev_(char *JOBZ, char *UPLO, int *N, double *A, int *LDA, double *W,
        double *WORK, int *LWORK, int *INFO);
 extern "C" void
@@ -174,7 +178,14 @@ const double &DenseMatrix::Elem(int i, int j) const
 
 void DenseMatrix::Mult(const double *x, double *y) const
 {
-   kernels::Mult(height, width, Data(), x, y);
+   const double *data = Read();
+   const int h = height;
+   const int w = width;
+
+   MFEM_FORALL(i, 1,
+   {
+      kernels::Mult(h, w, data, x, y);
+   });
 }
 
 void DenseMatrix::Mult(const Vector &x, Vector &y) const
@@ -182,7 +193,9 @@ void DenseMatrix::Mult(const Vector &x, Vector &y) const
    MFEM_ASSERT(height == y.Size() && width == x.Size(),
                "incompatible dimensions");
 
-   Mult((const double *)x, (double *)y);
+   const double *dx = x.Read();
+   double *dy = y.ReadWrite();
+   Mult(dx, dy);
 }
 
 double DenseMatrix::operator *(const DenseMatrix &m) const
@@ -2003,7 +2016,7 @@ void Mult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a)
    MFEM_ASSERT(a.Height() == b.Height() && a.Width() == c.Width() &&
                b.Width() == c.Height(), "incompatible dimensions");
 
-#ifdef MFEM_USE_LAPACK
+#if defined(MFEM_USE_LAPACK) && !defined(MFEM_USE_CUDA) && !defined(MFEM_USE_HIP)
    static char transa = 'N', transb = 'N';
    static double alpha = 1.0, beta = 0.0;
    int m = b.Height(), n = c.Width(), k = b.Width();
@@ -2014,10 +2027,13 @@ void Mult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a)
    const int ah = a.Height();
    const int aw = a.Width();
    const int bw = b.Width();
-   double *ad = a.Data();
-   const double *bd = b.Data();
-   const double *cd = c.Data();
-   kernels::Mult(ah,aw,bw,bd,cd,ad);
+   double *ad = a.ReadWrite();
+   const double *bd = b.Read();
+   const double *cd = c.Read();
+   MFEM_FORALL(i, 1,
+   {
+      kernels::Mult(ah, aw, bw, bd, cd, ad);
+   });
 #endif
 }
 
@@ -2859,6 +2875,155 @@ void AddMult_a_VVt(const double a, const Vector &v, DenseMatrix &VVt)
    }
 }
 
+void KronProd(const DenseMatrix & A, const DenseMatrix & B, DenseMatrix & C)
+{
+   const int ah = A.Height();
+   const int aw = A.Width();
+   const int bh = B.Height();
+   const int bw = B.Width();
+
+   C.SetSize(ah*bh,aw*bw);
+   const double * ad = A.Read();
+   const double * bd = B.Read();
+   double * cd = C.ReadWrite();
+
+   MFEM_FORALL(i, 1,
+   {
+      for (int ja = 0; ja<aw; ++ja)
+         for (int jb = 0; jb<bw; ++jb)
+            for (int ia = 0; ia<ah; ++ia)
+               for (int ib = 0; ib<bh; ++ib)
+                  cd[bh*ia + ib + ah*bh*(bw*ja + jb)]
+                     = ad[ia + ja * ah] * bd[ib + jb*bh];
+   });
+}
+
+#if 0 // this is a finer level parallel KronProd
+void KronProd2(const DenseMatrix & A, const DenseMatrix & B, DenseMatrix & C)
+{
+   const int ah = A.Height();
+   const int aw = A.Width();
+   const int bh = B.Height();
+   const int bw = B.Width();
+   const int ch = ah * bh;
+   const int cw = aw * bw;
+
+   C.SetSize(ch, cw);
+   const double *ad = A.Read();
+   const double *bd = B.Read();
+   double *cd = C.ReadWrite();
+
+   MFEM_FORALL(i, ch * cw,
+   {
+      const int jc = i / ch;
+      const int ic = i - jc * ch;
+      const int ja = jc / bw;
+      const int jb = jc - ja * bw;
+      const int ia = ic / bh;
+      const int ib = ic - ia * bh;
+      cd[jc * ch + ic] = ad[ja * ah + ia] * bd[jb * bh + ib];
+   });
+}
+#endif
+
+void KronMult(const DenseMatrix &A, const DenseMatrix &B, const Vector &r,
+              Vector &z)
+{
+   const int nA = A.Height();
+   const int mA = A.Width();
+   const int nB = B.Height();
+   const int mB = B.Width();
+   const int nr = r.Size();
+   MFEM_VERIFY(nr == mA*mB, "Wrong size of Vector r");
+   z.SetSize(nA*nB);
+#if !defined(MFEM_USE_CUDA)
+   DenseMatrix R(r.GetData(),mB,mA);
+   DenseMatrix X(nB,mA);
+   DenseMatrix Y(z.GetData(),nB,nA);
+   Mult(B,R,X);
+   MultABt(X,A,Y);
+#else
+   const double *ad = A.Read();
+   const double *bd = B.Read();
+   const double *rd = r.Read();
+   double *zd = z.Write();
+   MFEM_FORALL(i, 1,
+   {
+      kernels::KronMult(nA, mA, ad, nB, mB, bd, rd, zd);
+   });
+#endif
+}
+
+void KronMult(const DenseMatrix &A, const DenseMatrix &B, const DenseMatrix &R,
+              DenseMatrix & Z)
+{
+   const int nA = A.Height();
+   const int nB = B.Height();
+   const int nR = R.Height();
+   const int mR = R.Width();
+   Z.SetSize(nA*nB,mR);
+
+   Vector r,z;
+   double * dataR = R.Data();
+   for (int i = 0; i<mR; i++)
+   {
+      r.SetDataAndSize(&dataR[i*nR],nR);
+      KronMult(A,B,r,z);
+      Z.SetCol(i,z);
+   }
+}
+
+void KronMult(const DenseMatrix &A, const DenseMatrix &B, const DenseMatrix &C,
+              const Vector &r, Vector &z)
+{
+   const int nA = A.Height();
+   const int mA = A.Width();
+   const int nB = B.Height();
+   const int mB = B.Width();
+   const int nC = C.Height();
+   const int mC = C.Width();
+   const int nr = r.Size();
+   MFEM_VERIFY(nr == mA*mB*mC, "Wrong size of Vector r");
+   z.SetSize(nA*nB*nC);
+
+#if !defined(MFEM_USE_CUDA)
+   double * dataR = r.GetData();
+   DenseMatrix R(dataR,mC,mA*mB);
+   DenseMatrix X(nC,mA*mB);
+   Mult(C,R,X);
+   X.Transpose();
+   DenseMatrix Z(z.GetData(),mA*mB,nC);
+   KronMult(A,B,X,Z);
+   Z.Transpose();
+#else
+   const double *ad = A.Read();
+   const double *bd = B.Read();
+   const double *cd = C.Read();
+   const double *rd = r.Read();
+   double *zd = z.Write();
+   MFEM_FORALL(i, 1,
+   {
+      kernels::KronMult(nA, mA, ad, nB, mB, bd, nC, mC, cd, rd, zd);
+   });
+#endif
+}
+
+void KronMult(const Array<DenseMatrix *> & A, const Vector & r, Vector & z)
+{
+   int dim = A.Size();
+   if (dim == 2)
+   {
+      KronMult(*A[0],*A[1],r,z);
+   }
+   else if (dim == 3)
+   {
+      KronMult(*A[0],*A[1],*A[2], r,z);
+   }
+   else
+   {
+      MFEM_ABORT("KronMult::Wrong dimension");
+   }
+}
 
 bool LUFactors::Factor(int m, double TOL)
 {
@@ -3310,23 +3475,105 @@ DenseMatrixInverse::~DenseMatrixInverse()
    delete [] lu.ipiv;
 }
 
+void KronMult(const DenseMatrixInverse &A, const DenseMatrixInverse &B,
+              const Vector &r, Vector & z)
+{
+   // A and B are square matrices
+   z.SetSize(r.Size());
+   int nA = A.Height();
+   int nB = B.Height();
+   DenseMatrix R(r.GetData(),nB,nA);
+   DenseMatrix X(nB,nA);
+   B.Mult(R,X);
+   X.Transpose();
+   DenseMatrix Y(z.GetData(),nA,nB);
+   A.Mult(X,Y);
+   Y.Transpose();
+}
 
-DenseMatrixEigensystem::DenseMatrixEigensystem(DenseMatrix &m)
+void KronMult(const DenseMatrixInverse &A, const DenseMatrixInverse &B,
+              const DenseMatrix &R, DenseMatrix & Z)
+{
+   // A and B are square matrices
+   int nR = R.Height();
+   int mR = R.Width();
+   Z.SetSize(nR,mR);
+   Vector r(nR);
+   Vector z(nR);
+   double * dataR = R.GetData();
+   double * dataZ = Z.GetData();
+   for (int i = 0; i<mR; i++)
+   {
+      r.SetData(&dataR[i*nR]);
+      z.SetData(&dataZ[i*nR]);
+      KronMult(A,B,r,z);
+   }
+}
+
+void KronMult(const DenseMatrixInverse &A, const DenseMatrixInverse &B,
+              const DenseMatrixInverse &C, const Vector &r, Vector & z)
+{
+   // A, B and C are square matrices
+   int n = r.Size();
+   z.SetSize(n);
+   int nA = A.Height();
+   int nB = B.Height();
+   int nC = C.Height();
+   double * dataR = r.GetData();
+   DenseMatrix R(dataR,nC,nA*nB);
+   DenseMatrix X(nC,nA*nB);
+   C.Mult(R,X);
+   X.Transpose();
+   DenseMatrix Z(z.GetData(),nC,nA*nB);
+   KronMult(A,B,X,Z);
+   Z.Transpose();
+}
+
+void KronMult(const Array<DenseMatrixInverse *> & A, const Vector & r,
+              Vector & z)
+{
+   int dim = A.Size();
+   if (dim == 2)
+   {
+      KronMult(*A[0],*A[1],r,z);
+   }
+   else if (dim == 3)
+   {
+      KronMult(*A[0],*A[1],*A[2], r,z);
+   }
+   else
+   {
+      MFEM_ABORT("KronMult::Wrong dimension");
+   }
+}
+
+DenseMatrixEigensystem::DenseMatrixEigensystem(DenseMatrix &m, bool sym_)
    : mat(m)
 {
    n = mat.Width();
    EVal.SetSize(n);
+   EVali.SetSize(n);
    EVect.SetSize(n);
    ev.SetDataAndSize(NULL, n);
 
 #ifdef MFEM_USE_LAPACK
+   sym = sym_;
    jobz = 'V';
-   uplo = 'U';
    lwork = -1;
    double qwork;
-   dsyev_(&jobz, &uplo, &n, EVect.Data(), &n, EVal.GetData(),
-          &qwork, &lwork, &info);
-
+   if (sym)
+   {
+      uplo = 'U';
+      dsyev_(&jobz, &uplo, &n, EVect.Data(), &n, EVal.GetData(),
+             &qwork, &lwork, &info);
+   }
+   else
+   {
+      char jobvl = 'N';
+      int ldvl = 1;
+      dgeev_(&jobvl,&jobz,&n, mat.GetData(), &n, EVal.GetData(), EVali.GetData(),
+             nullptr, &ldvl, EVect.GetData(), &n, &qwork, &lwork, &info);
+   }
    lwork = (int) qwork;
    work = new double[lwork];
 #endif
@@ -3338,6 +3585,7 @@ DenseMatrixEigensystem::DenseMatrixEigensystem(
      n(other.n)
 {
 #ifdef MFEM_USE_LAPACK
+   sym = other.sym;
    jobz = other.jobz;
    uplo = other.uplo;
    lwork = other.lwork;
@@ -3356,13 +3604,24 @@ void DenseMatrixEigensystem::Eval()
 #endif
 
 #ifdef MFEM_USE_LAPACK
-   EVect = mat;
-   dsyev_(&jobz, &uplo, &n, EVect.Data(), &n, EVal.GetData(),
-          work, &lwork, &info);
-
+   if (sym)
+   {
+      EVect = mat;
+      dsyev_(&jobz, &uplo, &n, EVect.Data(), &n, EVal.GetData(),
+             work, &lwork, &info);
+   }
+   else
+   {
+      char jobvl = 'N';
+      int ldvl = 1;
+      DenseMatrix T = mat; // mat is overwritten by dgeev
+      dgeev_(&jobvl,&jobz,&n, T.GetData(), &n, EVal.GetData(), EVali.GetData(),
+             nullptr, &ldvl, EVect.GetData(), &n, work, &lwork, &info);
+   }
    if (info != 0)
    {
-      mfem::err << "DenseMatrixEigensystem::Eval(): DSYEV error code: "
+      string lpck = (sym) ? "DSYEV" : "DGEEV";
+      mfem::err << "DenseMatrixEigensystem::Eval(): " << lpck << "error code: "
                 << info << endl;
       mfem_error();
    }

linalg/densemat.hpp

@@ -523,6 +523,23 @@ void AddMult_a_VWt(const double a, const Vector &v, const Vector &w,
 /// VVt += a * v v^t
 void AddMult_a_VVt(const double a, const Vector &v, DenseMatrix &VVt);
 
+/// C = A ⊗ B
+void KronProd(const DenseMatrix & A, const DenseMatrix & B, DenseMatrix & C);
+void KronProd2(const DenseMatrix & A, const DenseMatrix & B, DenseMatrix & C);
+
+/// z = (A ⊗ B) r  = vec(B R A^T), where R := vec^-1 (r)
+void KronMult(const DenseMatrix &A, const DenseMatrix &B, const Vector &r,
+              Vector & z);
+
+/// z = (A ⊗ B) R
+void KronMult(const DenseMatrix &A, const DenseMatrix &B, const DenseMatrix &R,
+              DenseMatrix & Z);
+
+/// z = ( A ⊗ B ⊗ C ) r
+void KronMult(const DenseMatrix &A, const DenseMatrix &B, const DenseMatrix &C,
+              const Vector &r, Vector & z);
+
+void KronMult(const Array<DenseMatrix *> & A, const Vector & r, Vector & z);
 
 /** Class that can compute LU factorization of external data and perform various
     operations with the factored data. */
@@ -685,16 +702,33 @@ public:
    virtual ~DenseMatrixInverse();
 };
 
+/// z = (A^-1 ⊗ B^-1) r  = vec(B^-1 R A^-T), where R := vec^-1 (r)
+void KronMult(const DenseMatrixInverse &Ainv, const DenseMatrixInverse &Binv,
+              const Vector &r, Vector & z);
+
+/// z = (A^-1 ⊗ B^-1) R
+void KronMult(const DenseMatrixInverse &Ainv, const DenseMatrixInverse &Binv,
+              const DenseMatrix &R, DenseMatrix & Z);
+
+/// z = ( A^-1 ⊗ B^-1 ⊗ C^-1 ) r
+void KronMult(const DenseMatrixInverse &Ainv, const DenseMatrixInverse &Binv,
+              const DenseMatrixInverse &Cinv, const Vector &r, Vector & z);
+
+void KronMult(const Array<DenseMatrixInverse *> & A, const Vector & r,
+              Vector & z);
 
 class DenseMatrixEigensystem
 {
    DenseMatrix &mat;
    Vector      EVal;
+   // Possible non zero imaginary part of Eigenvalues
+   Vector      EVali;
    DenseMatrix EVect;
    Vector ev;
    int n;
 
 #ifdef MFEM_USE_LAPACK
+   bool sym;
    double *work;
    char jobz, uplo;
    int lwork, info;
@@ -702,10 +736,10 @@ class DenseMatrixEigensystem
 
 public:
 
-   DenseMatrixEigensystem(DenseMatrix &m);
+   DenseMatrixEigensystem(DenseMatrix &m, bool sym_ = false);
    DenseMatrixEigensystem(const DenseMatrixEigensystem &other);
    void Eval();
-   Vector &Eigenvalues() { return EVal; }
+   Vector &Eigenvalues(bool imag = false) { return imag ? EVali : EVal; }
    DenseMatrix &Eigenvectors() { return EVect; }
    double Eigenvalue(int i) { return EVal(i); }
    const Vector &Eigenvector(int i)

linalg/fdsolver.cpp

@@ -0,0 +1,143 @@
+// 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 "linalg.hpp"
+
+namespace mfem
+{
+
+void KronProdInvDiag(const Vector & a, const Vector & b, Vector & dinv)
+{
+   int n = a.Size(), m = b.Size();
+   dinv.SetSize(n*m);
+
+   for (int j = 0; j<m; j++)
+      for (int i = 0; i<n; i++)
+      {
+         dinv(i*m+j) = 1./(a(i) + b(j));
+      }
+}
+
+void KronProdInvDiag(const Vector & a, const Vector & b,
+                     const Vector & c, Vector & dinv)
+{
+   int n = a.Size(), m = b.Size(), l = c.Size();
+   dinv.SetSize(n*m*l);
+
+   for (int k = 0; k<l; k++)
+      for (int j = 0; j<m; j++)
+         for (int i = 0; i<n; i++)
+         {
+            dinv(i*m*l+j*l+k) = 1./(a(i) + b(j) + c(k));
+         }
+}
+
+void KronProdInvDiag(const Array<Vector *> & X, Vector & dinv)
+{
+   int dim = X.Size();
+   if (dim == 1)
+   {
+      int n = X[0]->Size();
+      dinv.SetSize(n);
+      for (int i = 0; i<n; i++) { dinv(i) = 1./(*X[0])(i); }
+   }
+   else if (dim == 2)
+   {
+      KronProdInvDiag(*X[0], *X[1], dinv);
+   }
+   else if (dim == 3)
+   {
+      KronProdInvDiag(*X[0], *X[1], *X[2], dinv);
+   }
+   else
+   {
+      MFEM_ABORT("KronProdInvDiag::Wrong dimension");
+   }
+}
+
+#ifdef MFEM_USE_LAPACK
+
+FDSolver::FDSolver(const Array<DenseMatrix *> & A,
+                   const Array<DenseMatrix *> & B)
+{
+   MFEM_ASSERT(A.Size() == B.Size(), "DenseFDSolver: Incompatible Dimensions");
+   dim = A.Size();
+
+   int solver_size = 1;
+   for (int i = 0; i<dim; i++)
+   {
+      MFEM_ASSERT(A[i]->Height() == A[i]->Width(),
+                  "DenseFDSolver: Matrix is not square");
+      MFEM_ASSERT(B[i]->Height() == B[i]->Width(),
+                  "DenseFDSolver: Matrix is not square");
+      MFEM_ASSERT(A[i]->Height() == B[i]->Height(),
+                  "DenseFDSolver: Matrices A and B have incompatible size");
+      solver_size *= A[i]->Height();
+   }
+   this->height = solver_size;
+   this->width  = solver_size;
+   if (solver_size) { Setup(A,B); }
+}
+
+void FDSolver::Setup(const Array<DenseMatrix *> & A,
+                     const Array<DenseMatrix *> & B)
+{
+   EigSystem.SetSize(dim);
+   eigv.SetSize(dim);
+   Array<Vector *> evalues(dim);
+   SQ.SetSize(dim);
+   DenseMatrix D;
+   for (int i = 0; i<dim; i++)
+   {
+      DenseMatrixInverse Minv(*B[i]);
+      Minv.Mult(*A[i],D);
+      EigSystem[i] = new DenseMatrixEigensystem(D);
+      EigSystem[i]->Eval();
+      evalues[i] = &EigSystem[i]->Eigenvalues();
+      eigv[i] = &EigSystem[i]->Eigenvectors();
+      DenseMatrixInverse Qinv(*eigv[i]);
+      DenseMatrix Sdinv;
+      Minv.GetInverseMatrix(Sdinv);
+      SQ[i] = new DenseMatrix;
+      Qinv.Mult(Sdinv,*SQ[i]);
+   }
+   KronProdInvDiag(evalues,dinv);
+}
+
+
+void FDSolver::Mult(const Vector & r,Vector & z) const
+{
+   MFEM_ASSERT(height == r.Size(),
+               "DenseFDSolver::Mult: Inconsistent vector size");
+   if (r.Size() == 0) { return; }
+   Vector rtemp;
+   KronMult(SQ,r,rtemp);
+   // 2. Diagonal solve;
+   rtemp *= dinv;
+   // 3. Modify RHS; z <-- (Q1 x Q2) rtemp
+   KronMult(eigv,rtemp,z);
+}
+
+FDSolver::~FDSolver()
+{
+   if (height)
+   {
+      for (int i=0; i<dim; i++)
+      {
+         delete SQ[i];
+         delete EigSystem[i];
+      }
+   }
+}
+
+#endif // MFEM_USE_LAPACK
+
+} // namespace mfem

linalg/fdsolver.hpp

@@ -0,0 +1,60 @@
+// 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_FDSOLVER
+#define MFEM_FDSOLVER
+
+#include "../config/config.hpp"
+#include "densemat.hpp"
+
+namespace mfem
+{
+
+/// Computes the inverse diagonal dinv = (a⊗I + I⊗b)^-1
+/// where a, b are diagonal matrices and I is the identity of the
+/// appropriate size
+void KronProdInvDiag(const Vector & a, const Vector & b, Vector & dinv);
+
+/// Computes the inverse diagonal dinv = (a⊗I⊗I + I⊗b⊗I + I⊗I⊗c)^-1
+/// where a, b, c are diagonal matrices and I is the identity of the
+/// appropriate size
+void KronProdInvDiag(const Vector & a, const Vector & b,
+                     const Vector & c, Vector & dinv);
+
+void KronProdInvDiag(const Array<Vector *> & X, Vector & dinv);
+
+#ifdef MFEM_USE_LAPACK
+
+/// In 2D it solves the system (A_0 ⊗ B_1 + B_0 ⊗ A_1) z = r
+/// In 3D it solves the system
+///     (A_0 ⊗ B_1 ⊗ B_2 + B_0 ⊗ A_1 ⊗ B_2 + B_0 ⊗ B_1 ⊗ A_2) z = r
+class FDSolver: public Solver
+{
+private:
+   int dim = 2;
+   Array<DenseMatrixEigensystem *> EigSystem;
+   Array<DenseMatrix *> eigv; // eigenvectors
+   Array<DenseMatrix *> SQ;
+   mutable Vector dinv;
+   void Setup(const Array<DenseMatrix *> & A, const Array<DenseMatrix *> & B);
+public:
+   FDSolver(const Array<DenseMatrix *> & A, const Array<DenseMatrix *> & B);
+   virtual void SetOperator(const Operator &op) {}
+   virtual void Mult(const Vector &r, Vector &z) const;
+   virtual ~FDSolver();
+};
+
+#endif // MFEM_USE_LAPACK
+
+} // mfem name space
+
+
+#endif // MFEM_FDSOLVER

linalg/kernels.hpp

@@ -160,7 +160,7 @@ double Norml2(const int size, const T *data)
     data of the input and output vectors. */
 template<typename TA, typename TX, typename TY>
 MFEM_HOST_DEVICE inline
-void Mult(const int height, const int width, TA *data, const TX *x, TY *y)
+void Mult(const int height, const int width, const TA *data, const TX *x, TY *y)
 {
    if (width == 0)
    {
@@ -170,7 +170,8 @@ void Mult(const int height, const int width, TA *data, const TX *x, TY *y)
       }
       return;
    }
-   TA *d_col = data;
+
+   TA *d_col = (TA *) data;
    TX x_col = x[0];
    for (int row = 0; row < height; row++)
    {
@@ -188,6 +189,52 @@ void Mult(const int height, const int width, TA *data, const TX *x, TY *y)
    }
 }
 
+template<typename TA, typename TB, typename TR, typename TZ>
+MFEM_HOST_DEVICE inline
+void KronMult(const int ah, const int aw, const TA *ad, const int bh, const int bw, const TB *bd, TR *r, TZ *z)
+{
+   for (int i = 0; i < bh; i++)
+   {
+      for (int l = 0; l < ah; l++)
+      {
+         TZ t1 = 0.0;
+         for (int j = 0; j < bw; j++)
+         {
+            const TB t2 = bd[i + j * bh];
+            for (int k = 0; k < aw; k++)
+            {
+               t1 += t2 * r[j + k * bw] * ad[l + k * ah];
+            }
+         }
+         z[i + l * bh] = t1;
+      }
+   }
+}
+
+template<typename TA, typename TB, typename TC, typename TR, typename TZ>
+MFEM_HOST_DEVICE inline
+void KronMult(const int ah, const int aw, const TA *ad, const int bh, const int bw, const TB *bd, const int ch, const int cw, const TC *cd, TR *r, TZ *z)
+{
+   for (int i = 0; i < ch; i++)
+   {
+      for (int l = 0; l < ah * bh; l++)
+      {
+         TZ t1 = 0.0;
+         for (int j = 0; j < cw; j++)
+         {
+            const TB t2 = cd[i + j * ch];
+            for (int k = 0; k < aw * bw; k++)
+            {
+               const TA ta = ad[(l / bh) + (k / bw) * ah];
+               const TB tb = bd[(l % bh) + (k % bw) * bh];
+               t1 += t2 * r[j + k * cw] * ta * tb;
+            }
+         }
+         z[i + l * ch] = t1;
+      }
+   }
+}
+
 /// Symmetrize a square matrix with given @a size and @a data: A -> (A+A^T)/2.
 template<typename T>
 MFEM_HOST_DEVICE inline

linalg/linalg.hpp

@@ -31,6 +31,7 @@
 #include "invariants.hpp"
 #include "constraints.hpp"
 #include "auxiliary.hpp"
+#include "fdsolver.hpp"
 
 #ifdef MFEM_USE_AMGX
 #include "amgxsolver.hpp"

tests/unit/CMakeLists.txt

@@ -28,6 +28,7 @@ set(UNIT_TESTS_SRCS
   linalg/test_complex_operator.cpp
   linalg/test_constrainedsolver.cpp
   linalg/test_direct_solvers.cpp
+  linalg/test_fdsolver.cpp
   linalg/test_hypre_ilu.cpp
   linalg/test_ilu.cpp
   linalg/test_matrix_block.cpp

tests/unit/linalg/test_fdsolver.cpp

@@ -0,0 +1,126 @@
+// 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;
+
+#ifdef MFEM_USE_LAPACK
+
+TEST_CASE("FDSolver",
+          "[FDSolver]")
+{
+   double tol = 1e-10;
+
+   // SPD matrices
+   DenseMatrix A0(
+   {
+      {
+         1.29919, 0.61256, 0.82545
+      },
+      {0.61256, 0.57891, 0.39662},
+      {0.82545, 0.39662, 0.57541}
+   });
+
+   DenseMatrix A1(
+   {
+      {0.748236, 0.701663, 0.607517, 0.236740},
+      {0.701663, 0.809316, 0.713186, 0.256070},
+      {0.607517, 0.713186, 0.794221, 0.233943},
+      {0.236740, 0.256070, 0.233943, 0.083129}
+   });
+
+   DenseMatrix B0(
+   {
+      {0.13483, 0.51389, 0.43052},
+      {0.51389, 2.26750, 1.86331},
+      {0.43052, 1.86331, 1.59869}
+   });
+
+   DenseMatrix B1(
+   {
+      {0.94177, 1.02400, 1.14743, 0.35723},
+      {1.02400, 1.79087, 1.78708, 0.78304},
+      {1.14743, 1.78708, 2.06259, 0.80837},
+      {0.35723, 0.78304, 0.80837, 1.01798}
+   });
+
+
+   SECTION("2D")
+   {
+      Array<DenseMatrix *> A(2), B(2);
+      A[0] = &A0;  A[1] = &A1;
+      B[0] = &B0;  B[1] = &B1;
+
+      Vector y(12); y.Randomize(1);
+      Vector x(12), diff(12);
+
+      FDSolver S(A,B);
+      S.Mult(y,x);
+
+      DenseMatrix C1, C;
+      KronProd(A0, B1, C1);
+      KronProd(B0, A1, C);
+
+      C.Add(1., C1);
+
+      DenseMatrixInverse Cinv(C);
+      Cinv.Mult(y,diff);
+
+      diff-=x;
+      REQUIRE(diff.Norml2() < tol);
+   }
+
+
+   SECTION("3D")
+   {
+      DenseMatrix A2(
+      {
+         {1.14593, 0.76119},
+         {0.76119, 0.78993}
+      });
+      DenseMatrix B2(
+      {
+         {0.88088, 0.37899},
+         {0.37899, 0.45096}
+      });
+      Array<DenseMatrix *> A(3), B(3);
+      A[0] = &A0;  A[1] = &A1; A[2] = &A2;
+      B[0] = &B0;  B[1] = &B1; B[2] = &B2;
+
+      Vector y(24); y.Randomize(1);
+      Vector x(24), diff(24);
+
+      FDSolver S(A,B);
+      S.Mult(y,x);
+
+      DenseMatrix Temp, C0, C1, C;
+
+      KronProd(A0, B1, Temp);
+      KronProd(Temp, B2, C0);
+      KronProd(B0, A1, Temp);
+      KronProd(Temp, B2, C1);
+      KronProd(B0, B1, Temp);
+      KronProd(Temp, A2, C);
+
+      C.Add(1.,C0);
+      C.Add(1.,C1);
+
+      DenseMatrixInverse Cinv(C);
+      Cinv.Mult(y,diff);
+
+      diff-=x;
+      REQUIRE(diff.Norml2() < tol);
+   }
+}
+
+#endif // if MFEM_USE_LAPACK

tests/unit/linalg/test_matrix_dense.cpp

@@ -12,6 +12,8 @@
 #include "mfem.hpp"
 #include "unit_tests.hpp"
 #include "linalg/dtensor.hpp"
+#include "nvToolsExt.h"
+#include "nvToolsExtCudaRt.h"
 
 using namespace mfem;
 
@@ -239,6 +241,212 @@ TEST_CASE("DenseMatrix A*B^T methods",
    }
 }
 
+TEST_CASE("KronMult methods",
+          "[DenseMatrix], [CUDA]")
+{
+   double tol = 1e-12;
+   int nA = 3, mA = 4;
+   int nB = 5, mB = 6;
+   DenseMatrix A(nA,mA);
+   DenseMatrix B(nB,mB);
+
+   for (int i = 0; i<nA; i++)
+      for (int j = 0; j<mA; j++)
+      {
+         A(i,j) = ((double)rand()/(double)RAND_MAX);
+      }
+
+   for (int i = 0; i<nB; i++)
+      for (int j = 0; j<mB; j++)
+      {
+         B(i,j) = ((double)rand()/(double)RAND_MAX);
+      }
+
+   DenseMatrix AB;
+   KronProd(A,B,AB);
+   AB.HostRead();
+
+   // (A ⊗ B) r
+   SECTION("KronMultABr")
+   {
+      nvtxRangePush("KronMultABr");
+      Vector r(mA*mB);
+      MFEM_VERIFY(r.Size() == AB.Width(), "Check r size");
+
+      r.HostReadWrite();
+      r.Randomize();
+
+      Vector z0(AB.Height());
+      AB.Mult(r,z0);
+      //z0.HostRead();
+      Vector z1;
+      KronMult(A,B,r,z1);
+      MFEM_VERIFY(z0.Size() == z1.Size(), "Check z1 size");
+      z0-=z1;
+      REQUIRE(z0.Norml2() < tol);
+      nvtxRangePop();
+   }
+
+
+   // (A ⊗ B) R
+   SECTION("KronMultABR")
+   {
+      nvtxRangePush("KronMultABR");
+      int nR = mA*mB;
+      int mR = 7;
+      DenseMatrix R(nR, mR);
+      R.HostReadWrite();
+
+      for (int i = 0; i<nR; i++)
+         for (int j = 0; j<mR; j++)
+         {
+            R(i,j) = ((double)rand()/(double)RAND_MAX);
+         }
+
+      DenseMatrix Z0(nA*nB,mR);
+      Mult(AB,R,Z0);
+
+      DenseMatrix Z1;
+      KronMult(A,B,R,Z1);
+      MFEM_VERIFY(Z0.Height() == Z1.Height() &&
+                  Z0.Width() == Z1.Width(), "Check z1 size");
+      Z0-=Z1;
+
+      REQUIRE(Z0.MaxMaxNorm() < tol);
+
+      nvtxRangePop();
+   }
+
+   // (A ⊗ B ⊗ C) r
+   SECTION("KronMultABCr")
+   {
+      nvtxRangePush("KronMultABCr");
+
+      int nC = 7, mC = 2;
+      DenseMatrix C(nC, mC);
+
+      C.HostReadWrite();
+
+      for (int i = 0; i<nC; i++)
+         for (int j = 0; j<mC; j++)
+         {
+            C(i,j) = ((double)rand()/(double)RAND_MAX);
+         }
+
+      DenseMatrix ABC;
+      KronProd(AB,C,ABC);
+
+      Vector r(mA*mB*mC);
+      r.HostReadWrite();
+      r.Randomize();
+
+      MFEM_VERIFY(r.Size() == ABC.Width(), "Check r size");
+      Vector z0(nA*nB*nC);
+      ABC.Mult(r,z0);
+
+      Vector z1;
+      KronMult(A,B,C,r,z1);
+      MFEM_VERIFY(z0.Size() == z1.Size(), "Check z1 size");
+      z0-=z1;
+      REQUIRE(z0.Norml2() < tol);
+
+      nvtxRangePop();
+   }
+}
+
+TEST_CASE("KronMultInv methods",
+          "[DenseMatrixInverse]")
+{
+   double tol = 1e-12;
+   int nA = 3;
+   int nB = 2;
+   DenseMatrix A(
+   {
+      { 1.0,  0.2, 3.4},
+      {-2.0, -1.0, 3.1},
+      { 0.7,  1.4,-0.9}
+   });
+   DenseMatrix B(
+   {
+      {-10.1, 5.7},
+      {-3.0,  4.2}
+   });
+
+   DenseMatrixInverse Ainv(A);
+   DenseMatrixInverse Binv(B);
+
+   DenseMatrix AB;
+   KronProd(A,B,AB);
+
+   // (A^-1 ⊗ B^-1) r
+   SECTION("KronMultInvABr")
+   {
+
+      Vector r(nA*nB); r.Randomize();
+      MFEM_VERIFY(r.Size() == AB.Width(), "Check r size");
+      Vector z0(AB.Height());
+      DenseMatrixInverse ABinv(AB);
+      ABinv.Mult(r,z0);
+
+      Vector z1;
+      KronMult(Ainv,Binv,r,z1);
+      MFEM_VERIFY(z0.Size() == z1.Size(), "Check z1 size");
+      z0-=z1;
+      REQUIRE(z0.Norml2() < tol);
+   }
+
+   // (A^-1 ⊗ B^-1) R
+   SECTION("KronMultInvABR")
+   {
+      int nR = nA*nB;
+      int mR = 7;
+      DenseMatrix R(nR, mR);
+      for (int i = 0; i<nR; i++)
+         for (int j = 0; j<mR; j++)
+         {
+            R(i,j) = ((double)rand()/(double)RAND_MAX);
+         }
+
+      DenseMatrixInverse ABinv(AB);
+      DenseMatrix Z0(nA*nB,mR);
+      ABinv.Mult(R,Z0);
+
+      DenseMatrix Z1;
+      KronMult(Ainv,Binv,R,Z1);
+      MFEM_VERIFY(Z0.Height() == Z1.Height() &&
+                  Z0.Width() == Z1.Width(), "Check z1 size");
+      Z0-=Z1;
+
+      REQUIRE(Z0.MaxMaxNorm() < tol);
+   }
+
+   // (A^-1 ⊗ B^-1 ⊗ C^-1) r
+   SECTION("KronMultInvABCr")
+   {
+      int nC = 4;
+      DenseMatrix C(
+      {
+         {-2.1, 1.6, -3.4,  17.5},
+         {-7.1, 1.3, -7.5, -12.5},
+         { 0.5, 5.7, -6.0, -0.5},
+         { 9.2, 0.3, -1.4, -14.9}
+      });
+      DenseMatrix ABC;
+      KronProd(AB,C,ABC);
+      DenseMatrixInverse ABCInv(ABC);
+      Vector r(nA*nB*nC); r.Randomize();
+      MFEM_VERIFY(r.Size() == ABC.Width(), "Check r size");
+      Vector z0(nA*nB*nC);
+      ABCInv.Mult(r,z0);
+
+      DenseMatrixInverse Cinv(C);
+      Vector z1;
+      KronMult(Ainv,Binv,Cinv,r,z1);
+      MFEM_VERIFY(z0.Size() == z1.Size(), "Check z1 size");
+      z0-=z1;
+      REQUIRE(z0.Norml2() < tol);
+   }
+}
 
 TEST_CASE("LUFactors RightSolve", "[DenseMatrix]")
 {
@@ -314,6 +522,88 @@ TEST_CASE("DenseTensor LinearSolve methods",
    }
 }
 
+#ifdef MFEM_USE_LAPACK
+
+TEST_CASE("EigenSystem methods",
+          "[DenseMatrix]")
+{
+   double tol = 1e-12;
+   SECTION("SPD Matrix")
+   {
+      DenseMatrix A({{0.56806, 0.29211, 0.48315, 0.70024},
+         {0.29211, 0.85147, 0.68123, 0.70689},
+         {0.48315, 0.68123, 1.07229, 1.02681},
+         {0.70024, 0.70689, 1.02681, 1.15468}
+      });
+      DenseMatrix V, AV(4);
+      Vector Lambda;
+      for (bool sym: { false, true })
+      {
+         DenseMatrixEigensystem  EigA(A,sym);
+         EigA.Eval();
+         V = EigA.Eigenvectors();
+         Lambda = EigA.Eigenvalues();
+         Mult(A,V,AV);
+         V.RightScaling(Lambda);
+         AV -= V;
+         REQUIRE(AV.MaxMaxNorm() < tol);
+      }
+   }
+
+   SECTION("Indefinite Matrix")
+   {
+      DenseMatrix A({{0.486278, 0.041135, 0.480727, 0.616026},
+         {0.523599, 0.119827, 0.087808, 0.415241},
+         {0.214454, 0.661631, 0.909626, 0.744259},
+         {0.107007, 0.630604, 0.077862, 0.221006}
+      });
+      DenseMatrixEigensystem  EigA(A);
+      EigA.Eval();
+
+      Vector Lambda_r, Lambda_i;
+      // Real part of eigenvalues
+      Lambda_r = EigA.Eigenvalues();
+      // Imag part of eigenvalues
+      Lambda_i = EigA.Eigenvalues(true);
+
+      DenseMatrix V;
+      V = EigA.Eigenvectors();
+      // Real part of eigenvectors
+      DenseMatrix Vr(4), Vi(4);
+      Vr.SetCol(0,V.GetColumn(0));
+      Vr.SetCol(1,V.GetColumn(1));
+      Vr.SetCol(2,V.GetColumn(1));
+      Vr.SetCol(3,V.GetColumn(3));
+
+      // Imag part of eigenvectors
+      Vector vi(4); V.GetColumn(2,vi);
+      Vi.SetCol(0,0.);
+      Vi.SetCol(1,vi); vi *= -1.;
+      Vi.SetCol(2,vi);
+      Vi.SetCol(3,0.);
+
+      // Check that A*V = V * Lambda
+      // or A * (V_r + i V_i ) = (V_r + i V_i)*(Lamda_r + i Lambda_i)
+      // or  A * V_r = V_r * Lambda_r -  V_i * Lambda_i
+      // and A * V_i = V_r ( Lambda_i +  V_i * Lambda_r
+      DenseMatrix AVr(4), AVi(4);
+      Mult(A,Vr, AVr);
+      Mult(A,Vi, AVi);
+
+      DenseMatrix Vrlr = Vr; Vrlr.RightScaling(Lambda_r);
+      DenseMatrix Vrli = Vr; Vrli.RightScaling(Lambda_i);
+      DenseMatrix Vilr = Vi; Vilr.RightScaling(Lambda_r);
+      DenseMatrix Vili = Vi; Vili.RightScaling(Lambda_i);
+
+      AVr -= Vrlr; AVr+= Vili;
+      AVi -= Vrli; AVi-= Vilr;
+
+      REQUIRE(AVr.MaxMaxNorm() < tol);
+      REQUIRE(AVi.MaxMaxNorm() < tol);
+   }
+}
+
+#endif // if MFEM_USE_LAPACK
 TEST_CASE("DenseTensor copy", "[DenseMatrix][DenseTensor]")
 {
    DenseTensor t1(2,3,4);