fix year

Algoim compilation warnings

fem/CMakeLists.txt

@@ -165,6 +165,7 @@ set(SRCS
   hyperbolic.cpp
   integrator.cpp
   bounds.cpp
+  blockform.cpp
   )
 
 set(HDRS
@@ -277,6 +278,7 @@ set(HDRS
   hyperbolic.hpp
   integrator.hpp
   bounds.hpp
+  blockform.hpp
   )
 
 if (MFEM_USE_SIDRE)

fem/blockform.cpp

@@ -0,0 +1,353 @@
+// 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.
+
+// Implementation of class BlockForm
+
+#include "../config/config.hpp"
+
+#include "fem.hpp"
+
+namespace mfem
+{
+
+void BlockForm::BuildProlongation()
+{
+   P = new BlockMatrix(dof_offsets, tdof_offsets);
+   R = new BlockMatrix(tdof_offsets, dof_offsets);
+   P->owns_blocks = 0;
+   R->owns_blocks = 0;
+
+   for (int i = 0; i<nblocks; i++)
+   {
+      const SparseMatrix * P_ = fes[i]->GetConformingProlongation();
+      P->SetBlock(i,i,const_cast<SparseMatrix*>(P_));
+      const SparseMatrix * R_ = fes[i]->GetRestrictionMatrix();
+      R->SetBlock(i,i,const_cast<SparseMatrix*>(R_));
+   }
+}
+
+void BlockForm::Finalize(int skip_zeros)
+{
+   if (mat) { mat->Finalize(skip_zeros); }
+   if (mat_e) { mat_e->Finalize(skip_zeros); }
+}
+
+void BlockForm::ConformingAssemble()
+{
+   Finalize(0);
+   if (!P) { BuildProlongation(); }
+
+   BlockMatrix * Pt = Transpose(*P);
+   BlockMatrix * PtA = mfem::Mult(*Pt, *mat);
+   // mat->owns_blocks = 0;
+   for (int i = 0; i<nblocks; i++)
+   {
+      for (int j = 0; j<nblocks; j++)
+      {
+         if (mat->IsZeroBlock(i,j)) { continue; }
+         if (Pt->IsZeroBlock(i,i))
+         {
+            PtA->SetBlock(i,j,&mat->GetBlock(i,j));
+         }
+      }
+   }
+   delete mat;
+   if (mat_e)
+   {
+      BlockMatrix *PtAe = mfem::Mult(*Pt, *mat_e);
+      mat_e->owns_blocks = 0;
+      for (int i = 0; i<nblocks; i++)
+      {
+         for (int j = 0; j<nblocks; j++)
+         {
+            if (mat_e->IsZeroBlock(i,j)) { continue; }
+            SparseMatrix * tmp = &mat_e->GetBlock(i,j);
+            if (Pt->IsZeroBlock(i,i))
+            {
+               PtAe->SetBlock(i,j,tmp);
+            }
+            else
+            {
+               delete tmp;
+            }
+         }
+      }
+      delete mat_e;
+      mat_e = PtAe;
+      mat_e->owns_blocks = 1;
+   }
+   delete Pt;
+
+   mat = mfem::Mult(*PtA, *P);
+   PtA->owns_blocks = 0;
+   for (int i = 0; i<nblocks; i++)
+   {
+      for (int j = 0; j<nblocks; j++)
+      {
+         if (PtA->IsZeroBlock(j,i)) { continue; }
+         SparseMatrix * tmp = &PtA->GetBlock(j,i);
+         if (P->IsZeroBlock(i,i))
+         {
+            mat->SetBlock(j,i,tmp);
+         }
+         else
+         {
+            delete tmp;
+         }
+      }
+   }
+   delete PtA;
+
+   if (mat_e)
+   {
+      BlockMatrix *PtAeP = mfem::Mult(*mat_e, *P);
+      mat_e->owns_blocks = 0;
+      for (int i = 0; i<nblocks; i++)
+      {
+         for (int j = 0; j<nblocks; j++)
+         {
+            if (mat_e->IsZeroBlock(j,i)) { continue; }
+            SparseMatrix * tmp = &mat_e->GetBlock(j,i);
+            if (P->IsZeroBlock(i,i))
+            {
+               PtAeP->SetBlock(j,i,tmp);
+            }
+            else
+            {
+               delete tmp;
+            }
+         }
+      }
+
+      delete mat_e;
+      mat_e = PtAeP;
+   }
+   height = mat->Height();
+   width = mat->Width();
+}
+
+
+BlockForm::BlockForm(const Array<FiniteElementSpace*> fes_ ): fes(
+      fes_)
+{
+   nblocks = fes.Size();
+   bforms.SetSize(nblocks,nblocks);
+   mforms.SetSize(nblocks,nblocks);
+   dof_offsets.Append(0);
+   tdof_offsets.Append(0);
+   for (int i = 0; i<nblocks; i++)
+   {
+      dof_offsets.Append(fes[i]->GetVSize());
+      tdof_offsets.Append(fes[i]->GetTrueVSize());
+      for (int j = 0; j<nblocks; j++)
+      {
+         bforms(i,j) = nullptr;
+         mforms(i,j) = nullptr;
+      }
+   }
+   dof_offsets.PartialSum();
+   tdof_offsets.PartialSum();
+   diag_policy = mfem::Operator::DIAG_ONE;
+
+}
+
+void BlockForm::SetBlock(BilinearForm * bform, int row_idx, int col_idx)
+{
+   MFEM_VERIFY((row_idx >=0 && row_idx < nblocks), "row index out of bounds");
+   MFEM_VERIFY((col_idx >=0 && col_idx < nblocks), "col index out of bounds");
+   MFEM_VERIFY(!mforms(row_idx,col_idx), "Entry has already been set");
+   MFEM_VERIFY(!bforms(row_idx,col_idx), "Entry has already been set");
+   bforms(row_idx,col_idx) = bform;
+}
+void BlockForm::SetBlock(MixedBilinearForm * mform, int row_idx,
+                         int col_idx)
+{
+   MFEM_VERIFY((row_idx >=0 && row_idx < nblocks), "row index out of bounds");
+   MFEM_VERIFY((col_idx >=0 && col_idx < nblocks), "col index out of bounds");
+   MFEM_VERIFY(!mforms(row_idx,col_idx), "Entry has already been set");
+   MFEM_VERIFY(!bforms(row_idx,col_idx), "Entry has already been set");
+   mforms(row_idx,col_idx) = mform;
+}
+
+/// Assemble the local matrix
+void BlockForm::Assemble(int skip_zeros)
+{
+   mat = new BlockMatrix(dof_offsets);
+   for (int i = 0; i<nblocks; i++)
+   {
+      int h = dof_offsets[i+1]-dof_offsets[i];
+      for (int j = 0; j<nblocks; j++)
+      {
+         int w = dof_offsets[j+1]-dof_offsets[j];
+         if (bforms(i,j))
+         {
+            bforms(i,j)->Assemble(skip_zeros);
+            MFEM_VERIFY(h = bforms(i,j)->Height(), "inconsistent height of bilinear form");
+            MFEM_VERIFY(w = bforms(i,j)->Width(), "inconsistent width of bilinear form");
+            mat->SetBlock(i,j,&bforms(i,j)->SpMat());
+         }
+         else if (mforms(i,j))
+         {
+            mforms(i,j)->Assemble(skip_zeros);
+            MFEM_VERIFY(h = mforms(i,j)->Height(),
+                        "inconsistent height of MixedBilinear form");
+            MFEM_VERIFY(w = mforms(i,j)->Width(),
+                        "inconsistent width of Mixedbilinear form");
+            mat->SetBlock(i,j,&mforms(i,j)->SpMat());
+         }
+         else
+         {
+            mat->SetBlock(i,j,nullptr);
+         }
+      }
+   }
+}
+
+void BlockForm::FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x,
+                                 Vector & b,
+                                 OperatorHandle &A, Vector &X,
+                                 Vector &B, int copy_interior)
+{
+   FormSystemMatrix(ess_tdof_list, A);
+
+   if (!P)
+   {
+      EliminateVDofsInRHS(ess_tdof_list, x, b);
+      X.MakeRef(x, 0, x.Size());
+      B.MakeRef(b, 0, b.Size());
+      if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
+   }
+   else // non conforming space
+   {
+      B.SetSize(P->Width());
+
+      P->MultTranspose(b, B);
+      real_t *data = b.GetData();
+      Vector tmp;
+      for (int i = 0; i<nblocks; i++)
+      {
+         if (P->IsZeroBlock(i,i))
+         {
+            int offset = tdof_offsets[i];
+            tmp.SetDataAndSize(&data[offset],tdof_offsets[i+1]-tdof_offsets[i]);
+            B.SetVector(tmp,offset);
+         }
+      }
+
+      X.SetSize(R->Height());
+
+      R->Mult(x, X);
+      data = x.GetData();
+      for (int i = 0; i<nblocks; i++)
+      {
+         if (R->IsZeroBlock(i,i))
+         {
+            int offset = tdof_offsets[i];
+            tmp.SetDataAndSize(&data[offset],tdof_offsets[i+1]-tdof_offsets[i]);
+            X.SetVector(tmp,offset);
+         }
+      }
+
+      EliminateVDofsInRHS(ess_tdof_list, X, B);
+      if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
+   }
+}
+
+void BlockForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
+                                 OperatorHandle &A)
+{
+   if (!mat_e)
+   {
+      bool conforming = true;
+      for (int i = 0; i<nblocks; i++)
+      {
+         const SparseMatrix *P_ = fes[i]->GetConformingProlongation();
+         if (P_)
+         {
+            conforming = false;
+            break;
+         }
+      }
+      if (!conforming) { ConformingAssemble(); }
+      const int remove_zeros = 0;
+      EliminateVDofs(ess_tdof_list, diag_policy);
+      Finalize(remove_zeros);
+   }
+   A.Reset(mat, false);
+}
+
+void BlockForm::RecoverFEMSolution(const Vector &X, Vector &x)
+{
+   if (!P)
+   {
+      x.SyncMemory(X);
+   }
+   else
+   {
+      x.SetSize(P->Height());
+      P->Mult(X, x);
+      real_t *data = X.GetData();
+      Vector tmp;
+      for (int i = 0; i<nblocks; i++)
+      {
+         if (P->IsZeroBlock(i,i))
+         {
+            int offset = tdof_offsets[i];
+            tmp.SetDataAndSize(&data[offset],tdof_offsets[i+1]-tdof_offsets[i]);
+            x.SetVector(tmp,offset);
+         }
+      }
+   }
+}
+
+void BlockForm::EliminateVDofs(const Array<int> &vdofs,
+                               Operator::DiagonalPolicy dpolicy)
+{
+   if (mat_e == NULL)
+   {
+      Array<int> offsets;
+
+      offsets.MakeRef( (P) ? tdof_offsets : dof_offsets);
+
+      mat_e = new BlockMatrix(offsets);
+      mat_e->owns_blocks = 1;
+      for (int i = 0; i<mat_e->NumRowBlocks(); i++)
+      {
+         int h = offsets[i+1] - offsets[i];
+         for (int j = 0; j<mat_e->NumColBlocks(); j++)
+         {
+            int w = offsets[j+1] - offsets[j];
+            mat_e->SetBlock(i,j,new SparseMatrix(h, w));
+         }
+      }
+   }
+
+   mat->EliminateRowCols(vdofs,mat_e,diag_policy);
+}
+
+void BlockForm::EliminateVDofsInRHS(
+   const Array<int> &vdofs, const Vector &x, Vector &b)
+{
+   mat_e->AddMult(x,b,-1.);
+   mat->PartMult(vdofs,x,b);
+}
+
+BlockForm::~BlockForm()
+{
+   delete mat_e;
+   mat_e = nullptr;
+   delete mat;
+   mat = nullptr;
+   delete P;
+   delete R;
+}
+
+};

fem/blockform.hpp

@@ -0,0 +1,88 @@
+// 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_BLOCKFORM
+#define MFEM_BLOCKFORM
+
+#include "../config/config.hpp"
+
+#include "fespace.hpp"
+#include "gridfunc.hpp"
+#include "bilinearform.hpp"
+
+namespace mfem
+{
+// square block forms
+class BlockForm
+{
+
+private:
+   int height, width;
+   int nblocks;
+   Array<int> dof_offsets;
+   Array<int> tdof_offsets;
+   // BilinearForms
+   Array2D<BilinearForm * > bforms;
+   Array2D<MixedBilinearForm * > mforms;
+   Array<FiniteElementSpace *> fes;
+
+   // Block Prolongation
+   BlockMatrix * P = nullptr;
+   // Block Restriction
+   BlockMatrix * R = nullptr;
+
+   BlockMatrix * mat = nullptr;
+   BlockMatrix * mat_e = nullptr;
+
+   void BuildProlongation();
+   void ConformingAssemble();
+
+   mfem::Operator::DiagonalPolicy diag_policy;
+
+public:
+
+   BlockForm(const Array<FiniteElementSpace*> pfes_ );
+
+   void SetBlock(BilinearForm * bform, int row_idx, int col_idx);
+   void SetBlock(MixedBilinearForm * mform, int row_idx, int col_idx);
+
+   /// Assemble the local matrix
+   void Assemble(int skip_zeros = 1);
+
+   void FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x, Vector & b,
+                         OperatorHandle &A, Vector &X,
+                         Vector &B, int copy_interior = 0);
+
+   void FormSystemMatrix(const Array<int> &ess_tdof_list,
+                         OperatorHandle &A);
+
+   void RecoverFEMSolution(const Vector &X, Vector &x);
+
+   ///  Finalizes the matrix initialization.
+   void Finalize(int skip_zeros = 1);
+
+   void EliminateVDofs(const Array<int> &vdofs,
+                       Operator::DiagonalPolicy dpolicy = Operator::DIAG_ONE);
+
+   void EliminateVDofsInRHS(const Array<int> &vdofs, const Vector &x, Vector &b);
+
+   void SetDiagonalPolicy(Operator::DiagonalPolicy policy)
+   {
+      diag_policy = policy;
+   }
+
+   /// Destroys bilinear form.
+   ~BlockForm();
+};
+
+}
+
+#endif

fem/conduitdatacollection.cpp

@@ -912,7 +912,7 @@ ConduitDataCollection::GridFunctionToBlueprintField(mfem::GridFunction *gf,
 
    if (vdim == 1) // scalar case
    {
-      n_field["values"].set_external(gf->GetData(),
+      n_field["values"].set_external(const_cast<real_t *>(gf->HostRead()),
                                      ndofs);
    }
    else // vector case
@@ -925,18 +925,18 @@ ConduitDataCollection::GridFunctionToBlueprintField(mfem::GridFunction *gf,
       int vdim_stride  = (ordering == Ordering::byNODES ? ndofs : 1);
 
       index_t offset = 0;
-      index_t stride = sizeof(double) * entry_stride;
+      index_t stride = sizeof(real_t) * entry_stride;
 
       for (int d = 0;  d < vdim; d++)
       {
          std::ostringstream oss;
          oss << "v" << d;
          std::string comp_name = oss.str();
-         n_field["values"][comp_name].set_external(gf->GetData(),
+         n_field["values"][comp_name].set_external(const_cast<real_t *>(gf->HostRead()),
                                                    ndofs,
                                                    offset,
                                                    stride);
-         offset +=  sizeof(double) * vdim_stride;
+         offset +=  sizeof(real_t) * vdim_stride;
       }
    }
 

fem/fem.hpp

@@ -50,8 +50,8 @@
 #include "dgmassinv.hpp"
 #include "hyperbolic.hpp"
 #include "bounds.hpp"
-
 #include "dfem/doperator.hpp"
+#include "blockform.hpp"
 
 #ifdef MFEM_USE_MPI
 #include "pfespace.hpp"

fem/intrules_cut.hpp

@@ -346,13 +346,13 @@ private:
       template<typename T>
       T operator() (const blitz::TinyVector<T,3>& x) const
       {
-         int el_order=el->GetOrder();
-         T u1[el_order+1];
-         T u2[el_order+1];
-         T u3[el_order+1];
-         TmplPoly_1D::CalcBernstein(el_order, x[0], u1);
-         TmplPoly_1D::CalcBernstein(el_order, x[1], u2);
-         TmplPoly_1D::CalcBernstein(el_order, x[2], u3);
+         const int el_order = el->GetOrder();
+         std::vector<T> u1(el_order+1);
+         std::vector<T> u2(el_order+1);
+         std::vector<T> u3(el_order+1);
+         TmplPoly_1D::CalcBernstein(el_order, x[0], u1.data());
+         TmplPoly_1D::CalcBernstein(el_order, x[1], u2.data());
+         TmplPoly_1D::CalcBernstein(el_order, x[2], u3.data());
 
          const Array<int>& dof_map=el->GetDofMap();
 
@@ -370,17 +370,17 @@ private:
       template<typename T>
       blitz::TinyVector<T,3> grad(const blitz::TinyVector<T,3>& x) const
       {
-         int el_order=el->GetOrder();
-         T u1[el_order+1];
-         T u2[el_order+1];
-         T u3[el_order+1];
-         T d1[el_order+1];
-         T d2[el_order+1];
-         T d3[el_order+1];
+         const int el_order = el->GetOrder();
+         std::vector<T> u1(el_order+1);
+         std::vector<T> u2(el_order+1);
+         std::vector<T> u3(el_order+1);
+         std::vector<T> d1(el_order+1);
+         std::vector<T> d2(el_order+1);
+         std::vector<T> d3(el_order+1);
 
-         TmplPoly_1D::CalcBernstein(el_order,x[0], u1, d1);
-         TmplPoly_1D::CalcBernstein(el_order,x[1], u2, d2);
-         TmplPoly_1D::CalcBernstein(el_order,x[2], u3, d3);
+         TmplPoly_1D::CalcBernstein(el_order,x[0], u1.data(), d1.data());
+         TmplPoly_1D::CalcBernstein(el_order,x[1], u2.data(), d2.data());
+         TmplPoly_1D::CalcBernstein(el_order,x[2], u3.data(), d3.data());
 
          blitz::TinyVector<T,3> res(T(0.0),T(0.0),T(0.0));
 
@@ -415,11 +415,11 @@ private:
       template<typename T>
       T operator() (const blitz::TinyVector<T,2>& x) const
       {
-         int el_order=el->GetOrder();
-         T u1[el_order+1];
-         T u2[el_order+1];
-         TmplPoly_1D::CalcBernstein(el_order, x[0], u1);
-         TmplPoly_1D::CalcBernstein(el_order, x[1], u2);
+         const int el_order = el->GetOrder();
+         std::vector<T> u1(el_order+1);
+         std::vector<T> u2(el_order+1);
+         TmplPoly_1D::CalcBernstein(el_order, x[0], u1.data());
+         TmplPoly_1D::CalcBernstein(el_order, x[1], u2.data());
 
          const Array<int>& dof_map=el->GetDofMap();
 
@@ -437,14 +437,14 @@ private:
       template<typename T>
       blitz::TinyVector<T,2> grad(const blitz::TinyVector<T,2>& x) const
       {
-         int el_order=el->GetOrder();
-         T u1[el_order+1];
-         T u2[el_order+1];
-         T d1[el_order+1];
-         T d2[el_order+1];
+         const int el_order = el->GetOrder();
+         std::vector<T> u1(el_order+1);
+         std::vector<T> u2(el_order+1);
+         std::vector<T> d1(el_order+1);
+         std::vector<T> d2(el_order+1);
 
-         TmplPoly_1D::CalcBernstein(el_order,x[0], u1, d1);
-         TmplPoly_1D::CalcBernstein(el_order,x[1], u2, d2);
+         TmplPoly_1D::CalcBernstein(el_order,x[0], u1.data(), d1.data());
+         TmplPoly_1D::CalcBernstein(el_order,x[1], u2.data(), d2.data());
 
          blitz::TinyVector<T,2> res(T(0.0),T(0.0));
 

miniapps/performance/makefile

@@ -105,8 +105,11 @@ MFEM_PERF_CXXFLAGS_xlc = -mcpu=native
 # - Clang extra options:
 ifeq ($(MFEM_MACHINE),riscv64)
   MFEM_PERF_CXXFLAGS_clang += -march=rv64gc
-else ifneq ($(MFEM_MACHINE),arm64)
-  # -march=native is unavailable on clang/ARM64 as of 05/2021: support could be added later.
+else ifneq (,$(findstring ppc,$(MFEM_MACHINE)))
+  MFEM_PERF_CXXFLAGS_clang += -mcpu=native -mtune=native
+else ifeq ($(MFEM_MACHINE),arm64)
+  MFEM_PERF_CXXFLAGS_clang += -mcpu=native -mtune=native
+else
   MFEM_PERF_CXXFLAGS_clang += -march=native
 endif
 MFEM_PERF_CXXFLAGS_clang += $(PEDANTIC_FLAG) -Wall