Merge master in metal

.gitignore

@@ -28,6 +28,15 @@ CMakeFiles/
 # Typical build directory
 /build/
 
+# vscode settings
+/.vscode/
+
+# Metal objects and libraries
+general/metal.*
+*.metallib
+*.metal
+*.mlb
+
 # Generated files in main directory, config/ and docs/
 /deps.mk
 config/_config.hpp
@@ -122,6 +131,8 @@ examples/cond_mesh.*
 examples/port_mesh.*
 examples/port_mode.*
 
+# examples/metal-*
+
 examples/euler-*
 
 examples/amgx/ex1

.vscode/settings.json

@@ -0,0 +1,9 @@
+{
+    "editor.defaultFormatter": "chiehyu.vscode-astyle",
+    "cmake.environment": {
+        "MFEM_DEBUG": "1",
+        // "ASAN_OPTIONS": "detect_leaks=1",
+        "LSAN_OPTIONS": "suppressions=/Users/camierjs/.asan.sup"
+    },
+    "editor.formatOnSave": true
+}

CMakeLists.txt

@@ -711,6 +711,30 @@ endif()
 
 set(MFEM_CUSTOM_TARGET_PREFIX CACHE STRING "")
 
+#-------------------------------------------------------------------------------
+# Metal
+#-------------------------------------------------------------------------------
+if (MFEM_USE_METAL)
+  message(STATUS "[🤘METAL🤘] PROJECT_SOURCE_DIR: ${PROJECT_SOURCE_DIR}")
+  message(STATUS "[🤘METAL🤘] PROJECT_BINARY_DIR: ${PROJECT_BINARY_DIR}")
+  message(STATUS "[🤘METAL🤘] CMAKE_INSTALL_PREFIX: ${CMAKE_INSTALL_PREFIX}")
+  
+  set(mfem_metal_src "${PROJECT_SOURCE_DIR}/general/metal.metal")
+  set(mfem_metal_lib "${PROJECT_BINARY_DIR}/mfem.mlb")
+  message(STATUS "[🤘METAL🤘] mfem_metal_lib: 🎯🎯🎯 ${mfem_metal_lib} 🎯🎯🎯")
+
+  add_custom_command( 
+      OUTPUT ${mfem_metal_lib} 
+      DEPENDS ${mfem_metal_src}
+      COMMAND /usr/bin/xcrun -sdk macosx metal -O3 -o ${mfem_metal_lib} ${mfem_metal_src} 
+  )
+
+  install(FILES ${PROJECT_BINARY_DIR}/mfem.mlb DESTINATION ${CMAKE_INSTALL_PREFIX})
+
+  add_custom_target(MFEM_METAL_LIB ALL DEPENDS ${mfem_metal_lib})  
+  add_dependencies(mfem MFEM_METAL_LIB)
+endif(MFEM_USE_METAL)
+
 #-------------------------------------------------------------------------------
 # Examples, miniapps, benchmarks and testing
 #-------------------------------------------------------------------------------

config/cmake/config.hpp.in

@@ -173,6 +173,9 @@
 // Enable functionality based on the libCEED library.
 #cmakedefine MFEM_USE_CEED
 
+// Build the Metal-enabled version of the MFEM library.
+#cmakedefine MFEM_USE_METAL
+
 // Enable functionality based on the Caliper library.
 #cmakedefine MFEM_USE_CALIPER
 

config/config.hpp.in

@@ -19,16 +19,16 @@
 // #define MFEM_VERSION_STRING "@MFEM_VERSION_STRING@"
 
 // MFEM version type, see the MFEM_VERSION_TYPE_* constants below.
-#define MFEM_VERSION_TYPE ((MFEM_VERSION)%2)
+#define MFEM_VERSION_TYPE ((MFEM_VERSION) % 2)
 
 // MFEM version type constants.
 #define MFEM_VERSION_TYPE_RELEASE 0
 #define MFEM_VERSION_TYPE_DEVELOPMENT 1
 
 // Separate MFEM version numbers for major, minor, and patch.
-#define MFEM_VERSION_MAJOR ((MFEM_VERSION)/10000)
-#define MFEM_VERSION_MINOR (((MFEM_VERSION)/100)%100)
-#define MFEM_VERSION_PATCH ((MFEM_VERSION)%100)
+#define MFEM_VERSION_MAJOR ((MFEM_VERSION) / 10000)
+#define MFEM_VERSION_MINOR (((MFEM_VERSION) / 100) % 100)
+#define MFEM_VERSION_PATCH ((MFEM_VERSION) % 100)
 
 // The absolute path of the MFEM source prefix.
 // #define MFEM_SOURCE_DIR "@MFEM_SOURCE_DIR@"
@@ -156,6 +156,9 @@
 // Enable MFEM functionality based on the GSLIB library.
 // #define MFEM_USE_GSLIB
 
+// Build the Darwin GPU/Metal-enabled version of the MFEM library.
+// #define MFEM_USE_METAL
+
 // Build the NVIDIA GPU/CUDA-enabled version of the MFEM library.
 // Requires a CUDA compiler (nvcc).
 // #define MFEM_USE_CUDA
@@ -210,4 +213,4 @@
 // Enable the Enzyme LLVM plugin
 // #define MFEM_USE_ENZYME
 
-#endif // MFEM_CONFIG_HEADER
+#endif  // MFEM_CONFIG_HEADER

config/config.mk.in

@@ -50,6 +50,7 @@ MFEM_USE_CONDUIT       = @MFEM_USE_CONDUIT@
 MFEM_USE_PUMI          = @MFEM_USE_PUMI@
 MFEM_USE_HIOP          = @MFEM_USE_HIOP@
 MFEM_USE_GSLIB         = @MFEM_USE_GSLIB@
+MFEM_USE_METAL         = @MFEM_USE_METAL@
 MFEM_USE_CUDA          = @MFEM_USE_CUDA@
 MFEM_USE_HIP           = @MFEM_USE_HIP@
 MFEM_USE_RAJA          = @MFEM_USE_RAJA@

config/defaults.cmake

@@ -57,6 +57,7 @@ option(MFEM_USE_HIP "Enable HIP" OFF)
 option(MFEM_USE_OCCA "Enable OCCA" OFF)
 option(MFEM_USE_RAJA "Enable RAJA" OFF)
 option(MFEM_USE_CEED "Enable CEED" OFF)
+option(MFEM_USE_METAL "Enable METAL" OFF)
 option(MFEM_USE_UMPIRE "Enable Umpire" OFF)
 option(MFEM_USE_SIMD "Enable use of SIMD intrinsics" OFF)
 option(MFEM_USE_ADIOS2 "Enable ADIOS2" OFF)

config/defaults.mk

@@ -41,6 +41,19 @@ INSTALL = /usr/bin/install
 STATIC = YES
 SHARED = NO
 
+# METAL configuration options
+METAL_CXX = clang++
+ifeq (YES,$(MFEM_USE_METAL))
+   BASE_FLAGS = -std=c++17
+   # specific options for development, should be reverted to -O3 for release
+   OPTIM_FLAGS = -Wall -O2 -g $(BASE_FLAGS) 
+   # the header file Metal.hpp should be in the root directory
+   # MFEM_TPLFLAGS += -Imetal-cpp
+endif
+METAL_FLAGS = -fno-objc-arc
+METAL_LIBS = -framework Metal -framework MetalKit -framework Cocoa\
+ -framework Foundation -framework CoreGraphics -L/opt/homebrew/lib -lfmt
+
 # CUDA configuration options
 #
 # If you set MFEM_USE_ENZYME=YES, CUDA_CXX has to be configured to use cuda with

examples/CMakeLists.txt

@@ -249,3 +249,8 @@ endif()
 if(MFEM_USE_MOONOLITH)
   add_subdirectory(moonolith)
 endif()
+
+# Include the examples/metal directory if Metal is enabled.
+if (MFEM_USE_METAL)
+  add_subdirectory(metal)
+endif()

examples/makefile

@@ -22,7 +22,7 @@ MFEM_LIB_FILE = mfem_is_not_built
 
 SEQ_EXAMPLES = ex0 ex1 ex2 ex3 ex4 ex5 ex6 ex7 ex8 ex9 ex10 ex14 ex15 ex16 \
    ex17 ex18 ex19 ex20 ex21 ex22 ex23 ex24 ex25 ex26 ex27 ex28 ex29 ex30 \
-   ex31 ex33 ex34 ex36 ex37 ex38 ex39 ex40
+   ex31 ex33 ex34 ex36 ex37 ex38 ex39 ex40 metal
 PAR_EXAMPLES = ex0p ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex8p ex9p ex10p ex11p \
    ex12p ex13p ex14p ex15p ex16p ex17p ex18p ex19p ex20p ex21p ex22p ex24p \
    ex25p ex26p ex27p ex28p ex29p ex30p ex31p ex32p ex33p ex34p ex35p ex36p \

examples/metal/CMakeLists.txt

@@ -0,0 +1,37 @@
+# Copyright (c) 2010-2024, 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.
+
+set(METAL_EXAMPLES_SRCS)
+list(APPEND METAL_EXAMPLES_SRCS
+    axpy.cpp
+    )
+
+message(STATUS "[🤘METAL🤘] Adding target: metal")
+add_custom_target(metal ${CMAKE_CTEST_COMMAND} -R metal USES_TERMINAL)
+
+# Include the source directory where mfem.hpp and mfem-performance.hpp are.
+include_directories(BEFORE ${PROJECT_BINARY_DIR})
+
+# Add one executable per cpp file
+set(PREFIX metal_)
+add_mfem_examples(METAL_EXAMPLES_SRCS ${PREFIX} "" metal)
+
+# Add a test for each example 
+if (MFEM_ENABLE_TESTING)
+  foreach(SRC_FILE ${METAL_EXAMPLES_SRCS})
+    get_filename_component(SRC_FILENAME ${SRC_FILE} NAME)
+    string(REPLACE ".cpp" "" TEST_NAME ${PREFIX}${SRC_FILENAME})
+
+    # add_dependencies(${TEST_NAME} MFEM_METAL_LIB)
+    add_test(NAME ${TEST_NAME} COMMAND ${TEST_NAME})
+    message(STATUS "[🤘METAL🤘] Adding test: ${TEST_NAME}")
+  endforeach()
+endif()

examples/metal/axpy.cpp

@@ -0,0 +1,146 @@
+// Copyright (c) 2010-2024, 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"
+
+#define DBG_COLOR ::debug::kPaleGreen
+#include "../../general/debug.hpp"
+
+using mfem::real_t;
+using mfem::Vector;
+
+//// ///////////////////////////////////////////////////////////////////////////
+template <class T>
+std::enable_if_t<!std::numeric_limits<T>::is_integer, bool>
+AlmostEq(T x, T y, T tolerance = 10.0*std::numeric_limits<T>::epsilon())
+{
+   const T neg = std::abs(x - y);
+   constexpr T min = std::numeric_limits<T>::min();
+   constexpr T eps = std::numeric_limits<T>::epsilon();
+   const T min_abs = std::min(std::abs(x), std::abs(y));
+   if (std::abs(min_abs) == 0.0)
+   {
+      return neg < eps;
+   }
+   return (neg / (1.0 + std::max(min, min_abs))) < tolerance;
+}
+
+//// ///////////////////////////////////////////////////////////////////////////
+bool equalArray(const real_t *x, const real_t *y, size_t N)
+{
+   for (unsigned long index = 0; index < N; index++)
+   {
+      // dbg("x:{}, y:{}", x[index], y[index]);
+      if (!AlmostEq(x[index], y[index]))
+      {
+         printf("Compute ERROR: index=%lu x=%e vs y=%e\n", index, x[index],
+                y[index]);
+         return false;
+      };
+   }
+   return true;
+}
+
+/// ////////////////////////////////////////////////////////////////////////////
+int main(int argc, char *argv[])
+{
+   dbg();
+   mfem::OptionsParser args(argc, argv);
+   const char *device_config = "metal";
+
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(std::cout);
+      return 1;
+   }
+   args.PrintOptions(std::cout);
+
+   mfem::Device device(device_config);
+   device.Print();
+
+   srand(time(nullptr));
+
+   constexpr auto N = 1 * 1024 * 1024;
+   constexpr auto K = 1;
+
+   mfem::Vector h_a(N), h_b(N), h_c(N);
+   h_a.UseDevice(false), h_b.UseDevice(false), h_c.UseDevice(false);
+   h_a.Randomize(1), h_b.Randomize(2), h_c.Randomize(3);
+
+   mfem::Vector a(N), b(N), c(N);
+   a.UseDevice(true), b.UseDevice(true), c.UseDevice(true);
+   a.Randomize(1), b.Randomize(2), c.Randomize(4);
+
+   constexpr mfem::real_t alpha = M_PI, beta = M_SQRT1_2, gamma = M_LN10;
+
+   mfem::tic_toc.Clear(), mfem::tic_toc.Start();
+   for (int i = 0; i < K; i++)
+   {
+      add(h_a, alpha, h_b, h_c), h_c *= h_c, h_c *= alpha;
+      h_a /= beta, h_a /= h_b;
+      h_a -= beta, h_a -= h_b;
+      h_a += alpha, h_a += h_c;
+      h_a.Add(beta, h_b);
+      h_a.Set(gamma, h_a), h_c = h_a;
+      h_a.Neg();
+      h_a.Reciprocal(), h_c = h_a;
+      add(h_a, h_b, h_c), h_a = h_c;
+      add(h_a, alpha, h_b, h_c), h_a = h_c;
+      add(gamma, h_a, h_b, h_c), h_a = h_c;
+      add(gamma, h_a, beta, h_b, h_c), h_a = h_c;
+      h_a /= alpha, h_b /= alpha, h_c /= alpha;
+      subtract(h_b, h_a, h_c), h_a = h_c;
+      subtract(alpha, h_a, h_b, h_c), h_a = h_c;
+   }
+   mfem::tic_toc.Stop();
+   const auto cpu_time = mfem::tic_toc.RealTime();
+   dbg("\033[36mCPU time: {}", cpu_time);
+
+   mfem::tic_toc.Clear();
+   mfem::tic_toc.Start();
+   for (int i = 0; i < K; i++)
+   {
+      add(a, alpha, b, c), c *= c, c *= alpha;
+      a /= beta, a /= b;
+      a -= beta, a -= b;
+      a += alpha, a += c;
+      a.Add(beta, b);
+      a.Set(gamma, a), c = a;
+      a.Neg();
+      a.Reciprocal(), c = a;
+      add(a, b, c), a = c;
+      add(a, alpha, b, c), a = c;
+      add(gamma, a, b, c), a = c;
+      add(gamma, a, beta, b, c), a = c;
+      a /= alpha, b /= alpha, c /= alpha;
+      subtract(b, a, c), a = c;
+      subtract(alpha, a, b, c), a = c;
+   }
+   mfem::tic_toc.Stop();
+   const auto gpu_time = mfem::tic_toc.RealTime();
+   dbg("\033[32mGPU time: {}", gpu_time);
+
+   dbg("Speedup: {}", (int) (cpu_time / gpu_time));
+
+   if (!equalArray(c.HostRead(), h_c.HostRead(), N))
+   {
+      dbg("\033[31mKernel result error vs. CPU code!!");
+      return EXIT_FAILURE;
+   }
+   dbg("\033[32mKernel result is equal to CPU code");
+
+   // const real_t dot = a * b; dbg("dot: {}", dot);
+
+   return EXIT_SUCCESS;
+}

examples/metal/form.cpp

@@ -0,0 +1,62 @@
+#include "mfem.hpp"
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   const char *mesh_file = "../data/star.mesh";
+   int order = 2;
+   const char *device_config = "cpu";
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(std::cout);
+      return 1;
+   }
+   args.PrintOptions(std::cout);
+
+   Device device(device_config);
+   device.Print();
+
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+
+   H1_FECollection fec(order, dim);
+
+   FiniteElementSpace fespace(&mesh, &fec);
+
+   Array<int> ess_tdof_list;
+   if (mesh.bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(mesh.bdr_attributes.Max());
+      ess_bdr = 1;
+      fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   GridFunction x(&fespace), y(&fespace);
+   x = 0.0;
+
+   BilinearForm a(&fespace);
+   a.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+
+   ConstantCoefficient one(1.0);
+   a.AddDomainIntegrator(new DiffusionIntegrator(one));
+   //    a.AddDomainIntegrator(new BilinearFormIntegrator("u·v"));
+
+   a.Assemble();
+   a.Mult(x, y);
+
+   // ir
+   // fe, mesh = fes
+   // a = inner(grad(u), grad(v))*dx = ∇u·∇v *dx = ∇u·∇v
+   // a = inner(u, v)*dx = u·v *dx u·v
+
+   return 0;
+}

general/CMakeLists.txt

@@ -20,6 +20,7 @@ list(APPEND SRCS
   hash.cpp
   isockstream.cpp
   mem_manager.cpp
+  metal.cpp
   occa.cpp
   optparser.cpp
   osockstream.cpp
@@ -50,6 +51,8 @@ list(APPEND HDRS
   kdtree.hpp
   mem_alloc.hpp
   mem_manager.hpp
+  metal.hpp
+  metal.h
   occa.hpp
   forall.hpp
   optparser.hpp

general/backends.hpp

@@ -31,6 +31,10 @@
 #include "occa.hpp"
 #endif
 
+#ifdef MFEM_USE_METAL
+#include "metal.hpp"
+#endif
+
 #ifdef MFEM_USE_RAJA
 // The following two definitions suppress CUB and THRUST deprecation warnings
 // about requiring c++14 with c++11 deprecated but still supported (to be

general/debug.hpp

@@ -0,0 +1,178 @@
+#pragma once
+#include <cstring>
+// https://fmt.dev/11.0/api/
+#include </opt/homebrew/include/fmt/format.h>
+
+#include <iomanip>
+#include <iostream>
+#include <utility>
+
+namespace debug
+{
+
+///////////////////////////////////////////////////////////////////////////////
+// https://en.wikipedia.org/wiki/Web_colors#Extended_colors
+// http://www.calmar.ws/vim/256-xterm-24bit-rgb-color-chart.html
+// clang-format off
+enum color_names
+{
+   kBlack = 0, kNavyBlue, kDarkBlue, kMediumBlue, kBlue, kDarkGreen, kWebGreen, kTeal,
+   kDarkCyan, kDeepSkyBlue, kDarkTurquoise, kMediumSpringGreen, kGreen, kLime,
+   kSpringGreen, kAqua, kCyan, kMidnightBlue, kDodgerBlue, kLightSeaGreen, kForestGreen,
+   kSeaGreen, kDarkSlateGray, kLimeGreen, kMediumSeaGreen, kTurquoise, kRoyalBlue,
+   kSteelBlue, kDarkSlateBlue, kMediumTurquoise, kIndigo, kDarkOliveGreen, kCadetBlue,
+   kCornflower, kRebeccaPurple, kMediumAquamarine, kDimGray, kSlateBlue, kOliveDrab,
+   kSlateGray, kLightSlateGray, kMediumSlateBlue, kLawnGreen, kWebMaroon, kWebPurple,
+   kChartreuse, kAquamarine, kOlive, kWebGray, kSkyBlue, kLightSkyBlue, kBlueViolet,
+   kDarkRed, kDarkMagenta, kSaddleBrown, kDarkSeaGreen, kLightGreen, kMediumPurple,
+   kDarkViolet, kPaleGreen, kDarkOrchid, kYellowGreen, kPurple, kSienna, kBrown,
+   kDarkGray, kLightBlue, kGreenYellow, kPaleTurquoise, kMaroon, kLightSteelBlue,
+   kPowderBlue, kFirebrick, kDarkGoldenrod, kMediumOrchid, kRosyBrown, kDarkKhaki,
+   kGray, kSilver, kMediumVioletRed, kIndianRed, kPeru, kChocolate, kTan, kLightGray,
+   kThistle, kOrchid, kGoldenrod, kPaleVioletRed, kCrimson, kGainsboro, kPlum, kBurlywood,
+   kLightCyan, kLavender, kDarkSalmon, kViolet, kPaleGoldenrod, kLightCoral, kKhaki,
+   kAliceBlue, kHoneydew, kAzure, kSandyBrown, kWheat, kBeige, kWhiteSmoke, kMintCream,
+   kGhostWhite, kSalmon, kAntiqueWhite, kLinen, kLightGoldenrod, kOldLace, kRed,
+   kFuchsia, kMagenta, kDeepPink, kOrangeRed, kTomato, kHotPink, kCoral, kDarkOrange,
+   kLightSalmon, kOrange, kLightPink, kPink, kGold, kPeachPuff, kNavajoWhite, kMoccasin,
+   kBisque, kMistyRose, kBlanchedAlmond, kPapayaWhip, kLavenderBlush, kSeashell,
+   kCornsilk, kLemonChiffon, kFloralWhite, kSnow, kYellow, kLightYellow, kIvory, kWhite,
+   kNvidia
+};
+// clang-format on
+
+static constexpr int kNumHexColors = 146;
+static constexpr std::array<uint32_t, kNumHexColors> kHexColors =
+{
+   {
+      0x000000, 0x000080, 0x00008B, 0x0000CD, 0x0000FF, 0x006400, 0x008000,
+      0x008080, 0x008B8B, 0x00BFFF, 0x00CED1, 0x00FA9A, 0x00FF00, 0x00FF00,
+      0x00FF7F, 0x00FFFF, 0x00FFFF, 0x191970, 0x1E90FF, 0x20B2AA, 0x228B22,
+      0x2E8B57, 0x2F4F4F, 0x32CD32, 0x3CB371, 0x40E0D0, 0x4169E1, 0x4682B4,
+      0x483D8B, 0x48D1CC, 0x4B0082, 0x556B2F, 0x5F9EA0, 0x6495ED, 0x663399,
+      0x66CDAA, 0x696969, 0x6A5ACD, 0x6B8E23, 0x708090, 0x778899, 0x7B68EE,
+      0x7CFC00, 0x7F0000, 0x7F007F, 0x7FFF00, 0x7FFFD4, 0x808000, 0x808080,
+      0x87CEEB, 0x87CEFA, 0x8A2BE2, 0x8B0000, 0x8B008B, 0x8B4513, 0x8FBC8F,
+      0x90EE90, 0x9370DB, 0x9400D3, 0x98FB98, 0x9932CC, 0x9ACD32, 0xA020F0,
+      0xA0522D, 0xA52A2A, 0xA9A9A9, 0xADD8E6, 0xADFF2F, 0xAFEEEE, 0xB03060,
+      0xB0C4DE, 0xB0E0E6, 0xB22222, 0xB8860B, 0xBA55D3, 0xBC8F8F, 0xBDB76B,
+      0xBEBEBE, 0xC0C0C0, 0xC71585, 0xCD5C5C, 0xCD853F, 0xD2691E, 0xD2B48C,
+      0xD3D3D3, 0xD8BFD8, 0xDA70D6, 0xDAA520, 0xDB7093, 0xDC143C, 0xDCDCDC,
+      0xDDA0DD, 0xDEB887, 0xE0FFFF, 0xE6E6FA, 0xE9967A, 0xEE82EE, 0xEEE8AA,
+      0xF08080, 0xF0E68C, 0xF0F8FF, 0xF0FFF0, 0xF0FFFF, 0xF4A460, 0xF5DEB3,
+      0xF5F5DC, 0xF5F5F5, 0xF5FFFA, 0xF8F8FF, 0xFA8072, 0xFAEBD7, 0xFAF0E6,
+      0xFAFAD2, 0xFDF5E6, 0xFF0000, 0xFF00FF, 0xFF00FF, 0xFF1493, 0xFF4500,
+      0xFF6347, 0xFF69B4, 0xFF7F50, 0xFF8C00, 0xFFA07A, 0xFFA500, 0xFFB6C1,
+      0xFFC0CB, 0xFFD700, 0xFFDAB9, 0xFFDEAD, 0xFFE4B5, 0xFFE4C4, 0xFFE4E1,
+      0xFFEBCD, 0xFFEFD5, 0xFFF0F5, 0xFFF5EE, 0xFFF8DC, 0xFFFACD, 0xFFFAF0,
+      0xFFFAFA, 0xFFFF00, 0xFFFFE0, 0xFFFFF0, 0xFFFFFF, 0x76B900
+   }
+};
+
+///////////////////////////////////////////////////////////////////////////////
+constexpr size_t static_strlen(const char *str)
+{
+   return *str == '\0' ? 0 : static_strlen(str + 1) + 1;
+}
+
+constexpr uint8_t static_checksum8(const char *bfr)
+{
+   unsigned int chk = 0;
+   size_t len = static_strlen(bfr);
+   for (; len; len--, bfr++) { chk += static_cast<unsigned int>(*bfr); }
+   return static_cast<uint8_t>(chk);
+}
+
+constexpr char *static_strrnchr(const char *str, const char c, int n)
+{
+   size_t len = static_strlen(str);
+   char *p = const_cast<char *>(str) + len - 1;
+   for (; n; n--, p--, len--)
+   {
+      for (; len; p--, len--)
+         if (*p == c) { break; }
+      if (!len) { return nullptr; }
+      if (n == 1) { return p; }
+   }
+   return nullptr;
+}
+
+inline uint32_t static_color(const uint8_t COLOR) { return kHexColors[COLOR]; }
+
+///////////////////////////////////////////////////////////////////////////////
+struct Debug
+{
+   const bool debug = false;
+
+   inline Debug() = default;
+
+   inline Debug(const char *FILE, const int LINE, const char *FUNC,
+                uint8_t COLOR)
+      : debug(true)
+   {
+      const char *base = static_strrnchr(FILE, '/', 2);
+      const char *file = base ? base + 1 : FILE;
+      const uint32_t rgb = static_color(COLOR);
+      const uint8_t r = (rgb >> 16) & 0xFF, g = (rgb >> 8) & 0xFF, b = rgb & 0xFF;
+      std::cout << "\033[38;2;";
+      std::cout << std::to_string(r) << ";";
+      std::cout << std::to_string(g) << ";";
+      std::cout << std::to_string(b) << "m";
+      std::cout << 0 << std::setw(64) << file << ":";
+      std::cout << "\033[2m" << std::setw(4) << std::left << LINE << "\033[22m: ";
+      if (FUNC) { std::cout << "[" << FUNC << "] "; }
+      std::cout << std::right << "\033[1m";
+   }
+
+   inline ~Debug()
+   {
+      if (debug) { std::cout << "\033[m\n" << std::flush; }
+   }
+
+   template <typename T>
+   inline void operator<<(const T &arg) const noexcept
+   {
+      if (debug) { std::cout << arg; }
+   }
+
+   template <typename T>
+   inline void operator()(const T &arg) const noexcept
+   {
+      if (debug) { this->operator<<(arg); }
+   }
+
+   template <typename... Args>
+   inline void operator()(const char *fmt, Args &&...args) const noexcept
+   {
+      if (debug) { std::cout << fmt::format(fmt, std::forward<Args>(args)...); }
+   }
+
+   inline void operator()() const noexcept {}
+
+   static Debug Set(const char *FILE, const int LINE, const char *FUNC,
+                    uint8_t COLOR)
+   {
+      static bool env_dbg = false;
+      if (static bool ini = false; !std::exchange(ini, true))
+      {
+         env_dbg = (getenv("MFEM_DEBUG") != nullptr);
+      }
+      const bool debug = env_dbg;
+      return debug ? Debug(FILE, LINE, FUNC, COLOR) : Debug();
+   }
+};
+
+// Helpers to generate unique variable names
+#define DBG_FLF __FILE__, __LINE__, __FUNCTION__
+#define DBG_PRIVATE_NAME(prefix) DBG_PRIVATE_CONCAT(prefix, __LINE__)
+#define DBG_PRIVATE_CONCAT(a, b) DBG_PRIVATE_CONCAT2(a, b)
+#define DBG_PRIVATE_CONCAT2(a, b) a##b
+
+#ifndef DBG_COLOR
+#define DBG_COLOR ::debug::kBlack
+#endif
+
+// Debug console traces, unnamed
+#define dbg(...) ::debug::Debug::Set(DBG_FLF, DBG_COLOR).operator()(__VA_ARGS__)
+
+}  // namespace debug

general/device.cpp

@@ -22,6 +22,9 @@
 #include <string>
 #include <map>
 
+#define DBG_COLOR ::debug::kHotPink
+#include "general/debug.hpp"
+
 namespace mfem
 {
 
@@ -48,7 +51,8 @@ static const Backend::Id backend_list[Backend::NUM_BACKENDS] =
    Backend::CEED_CUDA, Backend::OCCA_CUDA, Backend::RAJA_CUDA, Backend::CUDA,
    Backend::CEED_HIP, Backend::RAJA_HIP, Backend::HIP, Backend::DEBUG_DEVICE,
    Backend::OCCA_OMP, Backend::RAJA_OMP, Backend::OMP,
-   Backend::CEED_CPU, Backend::OCCA_CPU, Backend::RAJA_CPU, Backend::CPU
+   Backend::CEED_CPU, Backend::OCCA_CPU, Backend::RAJA_CPU,
+   Backend::METAL, Backend::CPU
 };
 
 // Backend names listed by priority, high to low:
@@ -57,7 +61,8 @@ static const char *backend_name[Backend::NUM_BACKENDS] =
    "ceed-cuda", "occa-cuda", "raja-cuda", "cuda",
    "ceed-hip", "raja-hip", "hip", "debug",
    "occa-omp", "raja-omp", "omp",
-   "ceed-cpu", "occa-cpu", "raja-cpu", "cpu"
+   "ceed-cpu", "occa-cpu", "raja-cpu", "metal",
+   "cpu"
 };
 
 } // namespace mfem::internal
@@ -182,6 +187,8 @@ Device::~Device()
 
 void Device::Configure(const std::string &device, const int device_id)
 {
+   dbg();
+
    // If a device was configured via the environment, skip the configuration,
    // and avoid the 'singleton_device' to destroy the mm.
    if (device_env)
@@ -197,7 +204,7 @@ void Device::Configure(const std::string &device, const int device_id)
       bmap[internal::backend_name[i]] = internal::backend_list[i];
    }
    std::string::size_type beg = 0, end, option;
-   while (1)
+   while (true)
    {
       end = device.find(',', beg);
       end = (end != std::string::npos) ? end : device.size();
@@ -206,7 +213,7 @@ void Device::Configure(const std::string &device, const int device_id)
       if (option==std::string::npos) // No option
       {
          const std::string backend = bname;
-         std::map<std::string, Backend::Id>::iterator it = bmap.find(backend);
+         auto it = bmap.find(backend);
          MFEM_VERIFY(it != bmap.end(), "invalid backend name: '" << backend << '\'');
          Get().MarkBackend(it->second);
       }
@@ -215,7 +222,7 @@ void Device::Configure(const std::string &device, const int device_id)
          const std::string backend = bname.substr(0, option);
          const std::string boption = bname.substr(option+1);
          Get().device_option = strdup(boption.c_str());
-         std::map<std::string, Backend::Id>::iterator it = bmap.find(backend);
+         auto it = bmap.find(backend);
          MFEM_VERIFY(it != bmap.end(), "invalid backend name: '" << backend << '\'');
          Get().MarkBackend(it->second);
       }
@@ -410,6 +417,15 @@ static void CudaDeviceSetup(const int dev, int &ngpu)
 #endif
 }
 
+static void MetalDeviceSetup(const int dev, int &ngpu)
+{
+   dbg();
+   MFEM_CONTRACT_VAR(dev);
+   MFEM_CONTRACT_VAR(ngpu);
+   auto metal = MTL::CreateSystemDefaultDevice();
+   dbg("Running on {}", metal->name()->utf8String());
+}
+
 static void HipDeviceSetup(const int dev, int &ngpu)
 {
 #ifdef MFEM_USE_HIP
@@ -511,6 +527,7 @@ static void CeedDeviceSetup(const char* ceed_spec)
 
 void Device::Setup(const int device_id)
 {
+   dbg();
    MFEM_VERIFY(ngpu == -1, "the mfem::Device is already configured!");
 
    ngpu = 0;
@@ -543,6 +560,7 @@ void Device::Setup(const int device_id)
                "Only one CEED backend can be enabled at a time!");
 #endif
    if (Allows(Backend::CUDA)) { CudaDeviceSetup(dev, ngpu); }
+   if (Allows(Backend::METAL)) { MetalDeviceSetup(dev, ngpu); }
    if (Allows(Backend::HIP)) { HipDeviceSetup(dev, ngpu); }
    if (Allows(Backend::RAJA_CUDA) || Allows(Backend::RAJA_HIP))
    { RajaDeviceSetup(dev, ngpu); }

general/device.hpp

@@ -73,7 +73,8 @@ struct Backend
           (using separate host/device memory pools and host <-> device
           transfers) without any GPU hardware. As 'DEBUG' is sometimes used
           as a macro, `_DEVICE` has been added to avoid conflicts. */
-      DEBUG_DEVICE = 1 << 14
+      METAL = 1 << 14,
+      DEBUG_DEVICE = 1 << 15
    };
 
    /** @brief Additional useful constants. For example, the *_MASK constants can
@@ -81,7 +82,7 @@ struct Backend
    enum
    {
       /// Number of backends: from (1 << 0) to (1 << (NUM_BACKENDS-1)).
-      NUM_BACKENDS = 15,
+      NUM_BACKENDS = 16,
 
       /// Biwise-OR of all CPU backends
       CPU_MASK = CPU | RAJA_CPU | OCCA_CPU | CEED_CPU,
@@ -94,7 +95,7 @@ struct Backend
       /// Bitwise-OR of all CEED backends
       CEED_MASK = CEED_CPU | CEED_CUDA | CEED_HIP,
       /// Biwise-OR of all device backends
-      DEVICE_MASK = CUDA_MASK | HIP_MASK | DEBUG_DEVICE,
+      DEVICE_MASK = CUDA_MASK | HIP_MASK | METAL | DEBUG_DEVICE,
 
       /// Biwise-OR of all RAJA backends
       RAJA_MASK = RAJA_CPU | RAJA_OMP | RAJA_CUDA | RAJA_HIP,
@@ -145,7 +146,7 @@ private:
    /// Current Device MemoryClass
    MemoryClass device_mem_class = MemoryClass::HOST;
 
-   char *device_option = NULL;
+   char *device_option = nullptr;
    Device(Device const&);
    void operator=(Device const&);
    static Device& Get() { return device_singleton; }

general/forall.hpp

@@ -753,6 +753,35 @@ inline void ForallWrap(const bool use_dev, const int N, lambda &&body,
 template<typename lambda>
 inline void forall(int N, lambda &&body) { ForallWrap<1>(true, N, body); }
 
+// forall with METAL backend
+template<typename lambda, typename ...Args>
+inline void forall(int N, lambda &&body, const char* kernel_name,
+                   [[maybe_unused]] const char*kernel_ops,
+                   Args... args)
+{
+   if (Device::Allows(mfem::Backend::METAL))
+   {
+      return metal::Kernel_1D(N, kernel_name, kernel_ops, args...);
+   }
+
+   ForallWrap<1>(true, N, body, body);
+}
+
+// forall_switch with METAL backend
+template<typename lambda, typename ...Args>
+inline void forall_switch(bool use_dev, int N, lambda &&body,
+                          const char* kernel_name,
+                          [[maybe_unused]] const char*kernel_ops,
+                          Args... args)
+{
+   if (use_dev && Device::Allows(mfem::Backend::METAL))
+   {
+      return metal::Kernel_1D(N, kernel_name, kernel_ops, args...);
+   }
+
+   ForallWrap<1>(use_dev, N, body);
+}
+
 template<typename lambda>
 inline void forall_switch(bool use_dev, int N, lambda &&body)
 {

general/mem_manager.cpp

@@ -51,6 +51,14 @@
 #define MAP_ANONYMOUS MAP_ANON
 #endif
 
+#ifdef MFEM_USE_METAL
+#include "metal.hpp"
+#endif
+
+#define DBG_COLOR ::debug::kYellow
+#include "general/debug.hpp"
+
+
 // Internal debug option, useful for tracking some memory manager operations.
 // #define MFEM_TRACK_MEM_MANAGER
 
@@ -186,13 +194,22 @@ struct Alias
 };
 
 /// Maps for the Memory and the Alias classes
-typedef std::unordered_map<const void*, Memory> MemoryMap;
-typedef std::unordered_map<const void*, Alias> AliasMap;
+using MemoryMap = std::unordered_map<const void *, Memory>;
+using AliasMap = std::unordered_map<const void *, Alias>;
+
+/// Buffer class that holds the Metal buffer
+#ifdef MFEM_USE_METAL
+using Buffer = MTL::Buffer*;
+using BufferMap = std::unordered_map<const void *, Buffer>;
+#endif
 
 struct Maps
 {
    MemoryMap memories;
    AliasMap aliases;
+#ifdef MFEM_USE_METAL
+   BufferMap buffers;
+#endif
 };
 
 } // namespace mfem::internal
@@ -468,18 +485,21 @@ public:
    {
 #ifdef MFEM_USE_CUDA
       CuMemAllocHostPinned(ptr, bytes);
-#endif
-#ifdef MFEM_USE_HIP
+#elif MFEM_USE_HIP
       HipMemAllocHostPinned(ptr, bytes);
+#else
+      MFEM_CONTRACT_VAR(ptr);
+      MFEM_CONTRACT_VAR(bytes);
 #endif
    }
    void Dealloc(void *ptr) override
    {
 #ifdef MFEM_USE_CUDA
       CuMemFreeHostPinned(ptr);
-#endif
-#ifdef MFEM_USE_HIP
+#elif MFEM_USE_HIP
       HipMemFreeHostPinned(ptr);
+#else
+      MFEM_CONTRACT_VAR(ptr);
 #endif
    }
 };
@@ -566,6 +586,42 @@ public:
    { return std::memcpy(dst, src, bytes); }
 };
 
+#ifdef MFEM_USE_METAL
+/// The METAL device memory space
+class MetalDeviceMemorySpace: public DeviceMemorySpace
+{
+   MTL::Device *device;
+public:
+   MetalDeviceMemorySpace(): DeviceMemorySpace(),
+      device(MTL::CreateSystemDefaultDevice()) { }
+   ~MetalDeviceMemorySpace() { device->autorelease(); }
+   void Alloc(Memory &base) override
+   {
+      constexpr auto options = MTL::ResourceStorageModeManaged;
+      auto mtl_buffer = device->newBuffer(base.bytes, options);
+      base.d_ptr = mtl_buffer->contents();
+      maps->buffers[base.d_ptr] = mtl_buffer;
+   }
+   void Dealloc(Memory &) override { }
+   void *HtoD(void *dst, const void *src, size_t bytes) override
+   {
+      // dbg("src: {} dst: {} bytes: {}", src, dst, bytes);
+      // could inform the GPU that the data is ready
+      return std::memcpy(dst, src, bytes);
+   }
+   void *DtoD(void* dst, const void* src, size_t bytes) override
+   {
+      // dbg("src: {} dst: {} bytes: {}", src, dst, bytes);
+      return std::memcpy(dst, src, bytes);
+   }
+   void *DtoH(void *dst, const void *src, size_t bytes) override
+   {
+      // dbg("src: {} dst: {} bytes: {}", src, dst, bytes);
+      return std::memcpy(dst, src, bytes);
+   }
+};
+#endif // MFEM_USE_METAL
+
 #ifdef MFEM_USE_UMPIRE
 class UmpireMemorySpace
 {
@@ -782,6 +838,8 @@ private:
             return new CudaDeviceMemorySpace();
 #elif defined(MFEM_USE_HIP)
             return new HipDeviceMemorySpace();
+#elif defined(MFEM_USE_METAL)
+            return new MetalDeviceMemorySpace();
 #else
             MFEM_ABORT("No device memory controller!");
             break;
@@ -1502,6 +1560,21 @@ void *MemoryManager::GetDevicePtr(const void *h_ptr, size_t bytes,
    return mem.d_ptr;
 }
 
+#ifdef MFEM_USE_METAL
+const MTL::Buffer *MemoryManager::GetDeviceBfr(void *d_ptr)
+{
+   MFEM_VERIFY(d_ptr, "cannot get device buffer for NULL device pointer");
+   MFEM_VERIFY(maps->buffers.find(d_ptr) != maps->buffers.end(),
+               "device buffer not found");
+   return maps->buffers.at(d_ptr);
+}
+
+const MTL::Buffer *MemoryManager::GetDeviceBfr(const void *d_ptr)
+{
+   return GetDeviceBfr(const_cast<void*>(d_ptr));
+}
+#endif
+
 void *MemoryManager::GetAliasDevicePtr(const void *alias_ptr, size_t bytes,
                                        bool copy)
 {

general/mem_manager.hpp

@@ -19,6 +19,10 @@
 #include <type_traits> // std::is_const
 #include <cstddef> // std::max_align_t
 
+#ifdef MFEM_USE_METAL
+#include "metal.h"
+#endif
+
 #ifdef MFEM_USE_MPI
 // Enable internal hypre timing routines
 #define HYPRE_TIMING
@@ -855,6 +859,12 @@ public:
    /// returning the number of printed pointers
    int PrintAliases(std::ostream &out = mfem::out);
 
+#ifdef MFEM_USE_METAL
+   /// Return the corresponding device buffer of d_ptr.
+   static const MTL::Buffer *GetDeviceBfr(void *d_ptr);
+   static const MTL::Buffer *GetDeviceBfr(const void *d_ptr);
+#endif
+
    static MemoryType GetHostMemoryType() { return host_mem_type; }
    static MemoryType GetDeviceMemoryType() { return device_mem_type; }
 

general/metal.cpp

@@ -0,0 +1,418 @@
+// Copyright (c) 2010-2024, 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 <string>
+#include <cctype>
+#include <stack>
+#include <string>
+#include <memory>
+#include <thread>
+
+#define DBG_COLOR ::debug::kMagenta
+#include "debug.hpp"
+
+#include "error.hpp"
+
+#define NS_PRIVATE_IMPLEMENTATION
+#define MTL_PRIVATE_IMPLEMENTATION
+// #define MTK_PRIVATE_IMPLEMENTATION
+// #define CA_PRIVATE_IMPLEMENTATION
+#include "metal.h"
+
+#include "metal.hpp"
+
+//// ///////////////////////////////////////////////////////////////////////////
+namespace NS
+{
+static NS::Error *error = nullptr;
+
+inline void Check(const bool test)
+{
+   if (test) { return; }
+   __builtin_printf("\033[31m%s\033[m\n",
+                    error->localizedDescription()->utf8String());
+   fflush(nullptr);
+   assert(false);
+}
+
+} // namespace NS
+
+//// ///////////////////////////////////////////////////////////////////////////
+namespace mfem
+{
+
+namespace metal
+{
+
+/// ///////////////////////////////////////////////////////////////////////////
+template <typename T, typename... Types>
+struct Base : public Base<T>, public Base<Types...>
+{
+   using Base<T>::Accept;
+   using Base<Types...>::Accept;
+};
+template <typename T> struct Base<T>
+{
+   virtual void Accept(T &) = 0;
+   virtual bool IsOps() { return false; };
+};
+
+template <typename T, typename... Types>
+struct Visitors : public Visitors<T>, public Visitors<Types...>
+{
+   using Visitors<T>::Visit;
+   using Visitors<Types...>::Visit;
+};
+template <typename T> struct Visitors<T>
+{
+   virtual void Visit(T &) = 0;
+};
+
+struct Node;
+
+struct Rule;
+struct Token;
+struct Const;
+struct Scalar;
+
+struct Visitor : public Visitors<Rule, Token, Const, Scalar> { };
+
+/// ///////////////////////////////////////////////////////////////////////////
+struct Node : public Base<Visitor>
+{
+   const char c;
+   std::shared_ptr<Node> left, right;
+   Node(const char c) : c(c), left(nullptr), right(nullptr) { }
+
+   void Accept(Visitor &) override = 0;
+   virtual ~Node() {}
+};
+
+/// ///////////////////////////////////////////////////////////////////////////
+struct Rule : public Node
+{
+   Rule(const char c) : Node(c) {}
+   void Accept(Visitor &me) override { me.Visit(*this); }
+   bool IsOps() override { return true;}
+};
+
+struct Token : public Node
+{
+   Token(const char c) : Node(c) {}
+   void Accept(Visitor &me) override { me.Visit(*this); }
+};
+
+struct Const : public Node
+{
+   Const(const char c) : Node(c) {}
+   void Accept(Visitor &me) override { me.Visit(*this); }
+};
+
+struct Scalar : public Node
+{
+   Scalar(const char c) : Node(c) {}
+   void Accept(Visitor &me) override { me.Visit(*this); }
+};
+
+using node_t = std::shared_ptr<Node>;
+
+/// ///////////////////////////////////////////////////////////////////////////
+[[maybe_unused]] static void Doze()
+{
+   std::this_thread::sleep_for(std::chrono::milliseconds(100));
+}
+
+/// ///////////////////////////////////////////////////////////////////////////
+// https://en.wikipedia.org/wiki/Shunting_yard_algorithm
+static auto ShuntingYard(const std::string& input)
+{
+   std::string infix(input);
+
+   auto priority = [](const char c) -> int
+   {
+      if (c == '=') { return 1; }
+      if (c == '*' || c == '/') { return 3; }
+      if (c == '+' || c == '-') { return 2; }
+      return 0;
+   };
+
+   auto isOperator = [](const char c) -> bool
+   {
+      return (c == '+' || c == '-' || c == '*' || c == '/' || c == '^' || c == '=');
+   };
+
+   char op_eq = 0;
+   std::stack<char> ops;
+   std::stack<node_t> ast;
+
+   auto NewRule = [&ast](char c) { ast.push(std::make_shared<Rule>(c)); };
+   auto NewToken = [&ast](char c) { ast.push(std::make_shared<Token>(c)); };
+   auto NewConst = [&ast](char c) { ast.push(std::make_shared<Const>(c)); };
+   auto NewScalar = [&ast](char c) { ast.push(std::make_shared<Scalar>(c)); };
+
+   auto NewOp = [&ast, &op_eq, NewRule,NewToken](char c)
+   {
+      // dbg("\033[33m[NewOp] c: {}", c), Doze();
+      if (c == '=')
+      {
+         const auto top = ast.top(); ast.pop();
+         NewToken('z'), ast.push(top);
+         if (op_eq == 0) { op_eq = '='; }
+      }
+      if (ast.size() == 0) { op_eq = c; return; }
+      assert(ast.size() > 1);
+      const auto right = ast.top(); ast.pop();
+      const auto left = ast.top(); ast.pop();
+      NewRule(c);
+      ast.top()->right = right, ast.top()->left = left;
+   };
+
+   for (char c : infix)
+   {
+      // dbg("new '{}'", c), Doze();
+      if (isdigit(c)) { NewConst(c); }
+      else if (std::isspace(c)) { continue; }
+      else if (isalpha(c) && c >= 'x') { NewToken(c); }
+      else if (isalpha(c) && c < 'x') { NewScalar(c); }
+      else if (isOperator(c))
+      {
+         while (!ops.empty() && priority(ops.top()) >= priority(c))
+         {
+            NewOp(ops.top()), ops.pop();
+         }
+         ops.push(c);
+      }
+      else if (c == '(') { ops.push(c); }
+      else if (c == ')')
+      {
+         while (!ops.empty() && ops.top() != '(') { NewOp(ops.top()), ops.pop(); }
+         if (!ops.empty() && ops.top() == '(') {  ops.pop(); }
+      }
+      else
+      {
+         MFEM_ABORT("unsupported character");
+      }
+   }
+
+   while (!ops.empty()) { NewOp(ops.top()), ops.pop(); }
+
+   return std::make_tuple(ast.top(), op_eq);
+}
+
+/// ///////////////////////////////////////////////////////////////////////////
+[[maybe_unused]] static void DfsPreOrder(const node_t &n, Visitor& m)
+{
+   if (!n) { return; }
+   n->Accept(m);
+   if (n->left) { DfsPreOrder(n->left, m); }
+   if (n->right) { DfsPreOrder(n->right, m); }
+}
+
+/// ///////////////////////////////////////////////////////////////////////////
+static void DfsInOrder(const node_t &n, Visitor& m)
+{
+   if (!n) { return; }
+   if (n->left) { DfsInOrder(n->left, m); }
+   n->Accept(m);
+   if (n->right) { DfsInOrder(n->right, m); }
+}
+
+/// ///////////////////////////////////////////////////////////////////////////
+[[maybe_unused]] static void DfsPostOrder(const node_t &n, Visitor& m)
+{
+   if (!n) { return; }
+   if (n->left) { DfsPostOrder(n->left, m); }
+   if (n->right) { DfsPostOrder(n->right, m); }
+   n->Accept(m);
+}
+
+/// ///////////////////////////////////////////////////////////////////////////
+void printAST(const node_t& node, int depth = 0)
+{
+   if (!node) { return; }
+
+   for (int i = 0; i < depth; ++i) { mfem::out << "  "; }
+   mfem::out << node->c << std::endl;
+
+   printAST(node->left, depth + 1);
+   printAST(node->right, depth + 1);
+}
+
+/// ///////////////////////////////////////////////////////////////////////////
+struct KernelDump : public Visitor
+{
+   void Visit(Rule &n) override { dbg("\033[33mrule {}",n.c); }
+   void Visit(Token &n) override { dbg("\033[33mtoken {}",n.c); }
+   void Visit(Const &n) override { dbg("\033[33mconst {}",n.c); }
+   void Visit(Scalar &n) override { dbg("\033[33mscalar {}",n.c); }
+};
+
+/// ///////////////////////////////////////////////////////////////////////////
+struct KernelSignature : public Visitor
+{
+   const std::string name;
+   std::ostringstream oss;
+
+   KernelSignature(std::string name): name(name)
+   {
+      oss << "\n[[kernel]] void " << name << "(\n";
+   }
+
+   void Visit(Rule &) override { /* nothing to do */}
+
+   void Visit(Token &n) override
+   {
+      if (n.c == 'z') { return; }
+      dbg("token {}", n.c);
+      oss << "\tdevice const float* " << n.c << ",\n";
+   }
+
+   void Visit(Scalar &n) override
+   {
+      dbg("const {}",n.c);
+      oss << "\tconstant float& " << n.c << ",\n";
+   }
+
+   void Visit(Const &) override { /* nothing to do */ }
+
+   std::string operator()()
+   {
+      oss << "\tdevice float* z,\n";
+      oss << "\tconst uint i [[thread_position_in_grid]]){\n";
+      return oss.str();
+   }
+};
+
+/// ///////////////////////////////////////////////////////////////////////////
+struct KernelBody : public Visitor
+{
+   const char op_eq;
+   std::stack<std::string> stack;
+
+   KernelBody(const char op_eq): op_eq(op_eq) { dbg("op_eq: {}", op_eq); }
+
+   void Visit(Rule &n) override
+   {
+      dbg("Rule {}", n.c);
+      assert(stack.size() > 1);
+      std::string op1 = stack.top();
+      stack.pop();
+      std::string op2 = stack.top();
+      stack.pop();
+
+      std::string op;
+      op = op2 + " ";
+      if (n.c == '=' && op_eq != '=') { op += op_eq; }
+      std::string eq(1,n.c);
+      op += eq + " " + op1;
+
+      if (n.c != '=') { stack.push("(" + op + ")"); }
+      else { stack.push(op); }
+   }
+
+   void Visit(Token &n) override { stack.push(std::string(1, n.c) + "[i]"); }
+
+   void Visit(Scalar &n) override { stack.emplace(1, n.c); }
+
+   void Visit(Const &n) override { stack.emplace(1, n.c); }
+
+   std::string operator()()
+   {
+
+      return std::string(op_eq == 0 ? "z[i] = ":"") + stack.top() + ";\n}";
+   }
+};
+
+/// ///////////////////////////////////////////////////////////////////////////
+std::string KernelOps(const char *name, const char *ops)
+{
+   dbg("Kernel '{}': \033[33m{}", name, ops);
+
+   std::string infix {ops};
+   auto [ast, op_eq] = ShuntingYard(infix);
+   printAST(ast);
+
+   KernelDump kd;
+   // dbg("\033[32mPreOrder");
+   // DfsPreOrder(ast, kd);
+
+   // dbg("\033[32mInOrder");
+   // DfsInOrder(ast, kd);
+
+   dbg("\033[32mPostOrder");
+   DfsPostOrder(ast, kd);
+
+
+   KernelSignature ker_signature(name);
+   KernelBody ker_body(op_eq);
+
+   DfsInOrder(ast, ker_signature);
+   DfsPostOrder(ast, ker_body);
+
+   std::ostringstream oss;
+   oss << ker_signature() << ker_body();
+
+   std::string kernel = oss.str();
+   dbg("Kernel: {}", kernel);
+   // dbg("\033[31mEXIT"); std::exit(0);
+
+   return kernel;
+}
+
+/// ////////////////////////////////////////////////////////////////////////////
+setup_t KernelSetup(const char* name, const char *src)
+{
+   static auto *device = MTL::CreateSystemDefaultDevice();
+#ifdef MFEM_USE_METAL_JIT
+   const MTL::CompileOptions *options = nullptr;
+   auto kernel_str = NS::String::string(src, NS::UTF8StringEncoding);
+   auto library = device->newLibrary(kernel_str, options, &NS::error);
+   NS::Check(library);
+#else // MFEM_USE_METAL_JIT
+   constexpr auto path = MFEM_SOURCE_DIR "/build/mfem.mlb";
+   // constexpr auto path = MFEM_INSTALL_DIR "/mfem.mlb";
+   const static auto filepath = NS::String::string(path, NS::ASCIIStringEncoding);
+   static auto library = device->newLibrary(filepath, &NS::error);
+   NS::Check(library);
+#endif // MFEM_USE_METAL_JIT
+
+   // function
+   auto function = library->newFunction(
+                      NS::String::string(name, NS::UTF8StringEncoding));
+   assert(function);
+
+   // kernel
+   auto kernel = device->newComputePipelineState(function, &NS::error);
+   assert(kernel);
+
+   // queue
+   auto Q = device->newCommandQueue();
+   assert(Q);
+
+   // commands
+   auto commands = Q->commandBuffer();
+   assert(commands);
+
+   // encoder
+   auto encoder = commands->computeCommandEncoder();
+   assert(encoder);
+
+   // enqueue the kernel
+   encoder->setComputePipelineState(kernel);
+
+   const auto MaxThreadsPerGroup = kernel->maxTotalThreadsPerThreadgroup();
+
+   return std::make_tuple(device, commands, encoder, MaxThreadsPerGroup);
+}
+
+} // namespace metal
+
+} // namespace mfem

general/metal.hpp

@@ -0,0 +1,96 @@
+// Copyright (c) 2010-2024, 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.
+
+#pragma once
+
+#include <tuple>
+
+#include "metal.h"
+
+#include "../config/config.hpp"
+#include "mem_manager.hpp"
+
+using mfem::real_t;
+
+namespace mfem
+{
+
+namespace metal
+{
+
+/// ////////////////////////////////////////////////////////////////////////////
+using setup_t =
+   std::tuple<
+   MTL::Device*,                 // device
+   MTL::CommandBuffer*,          // commands
+   MTL::ComputeCommandEncoder*,  // encoder
+   NS::UInteger>;                // maxTotalThreadsPerThreadgroup
+
+/// ////////////////////////////////////////////////////////////////////////////
+std::string KernelOps(const char *name, const char *ops);
+setup_t KernelSetup(const char* name, const char *src);
+
+/// ////////////////////////////////////////////////////////////////////////////
+template <typename F, typename D, typename E, typename... Args>
+void KernelApply(F func, D &dev, E &enc, Args&&... args)
+{
+   (func(dev, enc, std::forward<Args>(args)), ...);
+}
+
+/// ////////////////////////////////////////////////////////////////////////////
+template <typename... Args>
+void Kernel_1D(const size_t N, const char* name, const char *ops,
+               Args &&...args)
+{
+   const std::string src = KernelOps(name, ops);
+   auto [device, commands, encoder, KerMaxTPG] = KernelSetup(name, src.c_str());
+
+   // enqueue the argument buffers
+   int k = 0;
+   auto enqueue = [&k](auto &device, auto&encoder, auto &arg)
+   {
+      // dbg("\t#{} {}", k, typeid(arg).name());
+      using arg_t = std::decay_t<decltype(arg)>;
+      if constexpr (std::is_same_v<arg_t, real_t*> ||
+                    std::is_same_v<arg_t, const real_t*>)
+      {
+         encoder->setBuffer(mfem::MemoryManager::GetDeviceBfr(arg), 0, k++);
+      }
+      else if constexpr (std::is_same_v<arg_t, float>)
+      {
+         // turn real_t into a buffer
+         auto N_MTL = device->newBuffer(sizeof(float), MTL::ResourceStorageModeManaged);
+         *static_cast<float *>(N_MTL->contents()) = arg;
+         encoder->setBuffer(N_MTL, 0, k++);
+      }
+      else
+      {
+         MFEM_ABORT("unsupported type");
+      }
+   };
+
+   KernelApply(enqueue, device, encoder, args...);
+
+   const auto threadGroupSize = KerMaxTPG > N ? N : KerMaxTPG;
+   const auto threadsPerThreadgroup = MTL::Size::Make(threadGroupSize, 1, 1);
+   const auto threadsPerGrid = MTL::Size::Make(N, 1, 1);
+
+   // Encode the compute command
+   encoder->dispatchThreads(threadsPerGrid, threadsPerThreadgroup);
+   encoder->endEncoding();
+
+   commands->commit();
+   commands->waitUntilCompleted();
+}
+
+} // namespace metal
+
+} // namespace mfem

linalg/vector.cpp

@@ -152,7 +152,8 @@ Vector &Vector::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; });
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] = value; },
+   "linalg_vector_eq_c", "c", value, y);
    return *this;
 }
 
@@ -161,7 +162,8 @@ Vector &Vector::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; });
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] *= c; },
+   "linalg_vector_mul_eq_c", "*=c", c, y);
    return *this;
 }
 
@@ -171,9 +173,10 @@ Vector &Vector::operator*=(const Vector &v)
 
    const bool use_dev = UseDevice() || v.UseDevice();
    const int N = size;
+   const auto x = v.Read(use_dev);
    auto y = ReadWrite(use_dev);
-   auto x = v.Read(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] *= x[i]; });
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] *= x[i]; },
+   "linalg_vector_mul_eq_x", "*=x", x, y);
    return *this;
 }
 
@@ -183,7 +186,8 @@ Vector &Vector::operator/=(real_t c)
    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; });
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] *= m; },
+   "linalg_vector_div_eq_c", "*=c", m, y);
    return *this;
 }
 
@@ -193,9 +197,10 @@ Vector &Vector::operator/=(const Vector &v)
 
    const bool use_dev = UseDevice() || v.UseDevice();
    const int N = size;
+   const auto x = v.Read(use_dev);
    auto y = ReadWrite(use_dev);
-   auto x = v.Read(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] /= x[i]; });
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] /= x[i]; },
+   "linalg_vector_div_eq_x", "/=x", x, y);
    return *this;
 }
 
@@ -204,7 +209,8 @@ Vector &Vector::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; });
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] -= c; },
+   "linalg_vector_sub_eq_c", "-=c", c, y);
    return *this;
 }
 
@@ -214,9 +220,10 @@ Vector &Vector::operator-=(const Vector &v)
 
    const bool use_dev = UseDevice() || v.UseDevice();
    const int N = size;
+   const auto x = v.Read(use_dev);
    auto y = ReadWrite(use_dev);
-   auto x = v.Read(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] -= x[i]; });
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] -= x[i]; },
+   "linalg_vector_seq_x", "-=x", x, y);
    return *this;
 }
 
@@ -225,7 +232,8 @@ Vector &Vector::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; });
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] += c; },
+   "linalg_vector_add_eq_c", "+=c", c, y);
    return *this;
 }
 
@@ -235,9 +243,10 @@ Vector &Vector::operator+=(const Vector &v)
 
    const bool use_dev = UseDevice() || v.UseDevice();
    const int N = size;
+   const auto x = v.Read(use_dev);
    auto y = ReadWrite(use_dev);
-   auto x = v.Read(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] += x[i]; });
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] += x[i]; },
+   "linalg_vector_add_eq_x", "+=x", x, y);
    return *this;
 }
 
@@ -249,9 +258,10 @@ Vector &Vector::Add(const real_t a, const Vector &Va)
    {
       const int N = size;
       const bool use_dev = UseDevice() || Va.UseDevice();
+      const auto x = Va.Read(use_dev);
       auto y = ReadWrite(use_dev);
-      auto x = Va.Read(use_dev);
-      mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] += a * x[i]; });
+      mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] += a * x[i]; },
+      "linalg_vector_add_eq_ax", "+=a*x", a, x, y);
    }
    return *this;
 }
@@ -264,7 +274,8 @@ Vector &Vector::Set(const real_t a, const Vector &Va)
    const int N = size;
    auto x = Va.Read(use_dev);
    auto y = Write(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] = a * x[i]; });
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] = a * x[i]; },
+   "linalg_vector_set_ax", "a*x", a, x, y);
    return *this;
 }
 
@@ -299,7 +310,8 @@ void Vector::Neg()
    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] = -y[i]; });
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] = -y[i]; },
+   "linalg_vector_neg", "(0-y)", y, y);
 }
 
 void Vector::Reciprocal()
@@ -307,7 +319,8 @@ void Vector::Reciprocal()
    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] = 1.0/y[i]; });
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] = 1.0/y[i]; },
+   "linalg_vector_reciprocal", "1/y", y, y);
 }
 
 void add(const Vector &v1, const Vector &v2, Vector &v)
@@ -319,10 +332,11 @@ void add(const Vector &v1, const Vector &v2, Vector &v)
    const bool use_dev = v1.UseDevice() || v2.UseDevice() || v.UseDevice();
    const int N = v.size;
    // Note: get read access first, in case v is the same as v1/v2.
-   auto x1 = v1.Read(use_dev);
-   auto x2 = v2.Read(use_dev);
+   const auto x1 = v1.Read(use_dev);
+   const auto x2 = v2.Read(use_dev);
    auto y = v.Write(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] = x1[i] + x2[i]; });
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] = x1[i] + x2[i]; },
+   "linalg_vector_add_x_y", "x+y", x1, x2, y);
 #else
    #pragma omp parallel for
    for (int i = 0; i < v.size; i++)
@@ -351,13 +365,14 @@ void add(const Vector &v1, real_t alpha, const Vector &v2, Vector &v)
       const bool use_dev = v1.UseDevice() || v2.UseDevice() || v.UseDevice();
       const int N = v.size;
       // Note: get read access first, in case v is the same as v1/v2.
-      auto d_x = v1.Read(use_dev);
-      auto d_y = v2.Read(use_dev);
+      const auto d_x = v1.Read(use_dev);
+      const auto d_y = v2.Read(use_dev);
       auto d_z = v.Write(use_dev);
       mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
       {
          d_z[i] = d_x[i] + alpha * d_y[i];
-      });
+      },
+      "linalg_vector_axpy", "x+a*y", d_x, alpha, d_y, d_z);
 #else
       const real_t *v1p = v1.data, *v2p = v2.data;
       real_t *vp = v.data;
@@ -390,13 +405,14 @@ void add(const real_t a, const Vector &x, const Vector &y, Vector &z)
       const bool use_dev = x.UseDevice() || y.UseDevice() || z.UseDevice();
       const int N = x.size;
       // Note: get read access first, in case z is the same as x/y.
-      auto xd = x.Read(use_dev);
-      auto yd = y.Read(use_dev);
+      const auto xd = x.Read(use_dev);
+      const auto yd = y.Read(use_dev);
       auto zd = z.Write(use_dev);
       mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
       {
          zd[i] = a * (xd[i] + yd[i]);
-      });
+      },
+      "linalg_vector_a_xpy", "a*(x+y)", a, xd, yd, zd);
 #else
       const real_t *xp = x.data;
       const real_t *yp = y.data;
@@ -445,13 +461,14 @@ void add(const real_t a, const Vector &x,
       const bool use_dev = x.UseDevice() || y.UseDevice() || z.UseDevice();
       const int N = x.size;
       // Note: get read access first, in case z is the same as x/y.
-      auto xd = x.Read(use_dev);
-      auto yd = y.Read(use_dev);
+      const auto xd = x.Read(use_dev);
+      const auto yd = y.Read(use_dev);
       auto zd = z.Write(use_dev);
       mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
       {
          zd[i] = a * xd[i] + b * yd[i];
-      });
+      },
+      "linalg_vector_ax_plus_by", "a*x + b*y", a, xd, b, yd, zd);
 #else
       const real_t *xp = x.data;
       const real_t *yp = y.data;
@@ -475,13 +492,14 @@ void subtract(const Vector &x, const Vector &y, Vector &z)
    const bool use_dev = x.UseDevice() || y.UseDevice() || z.UseDevice();
    const int N = x.size;
    // Note: get read access first, in case z is the same as x/y.
-   auto xd = x.Read(use_dev);
-   auto yd = y.Read(use_dev);
+   const auto xd = x.Read(use_dev);
+   const auto yd = y.Read(use_dev);
    auto zd = z.Write(use_dev);
    mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
    {
       zd[i] = xd[i] - yd[i];
-   });
+   },
+   "linalg_vector_subtract", "x-y", xd, yd, zd);
 #else
    const real_t *xp = x.data;
    const real_t *yp = y.data;
@@ -514,13 +532,14 @@ void subtract(const real_t a, const Vector &x, const Vector &y, Vector &z)
       const bool use_dev = x.UseDevice() || y.UseDevice() || z.UseDevice();
       const int N = x.size;
       // Note: get read access first, in case z is the same as x/y.
-      auto xd = x.Read(use_dev);
-      auto yd = y.Read(use_dev);
+      const auto xd = x.Read(use_dev);
+      const auto yd = y.Read(use_dev);
       auto zd = z.Write(use_dev);
       mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
       {
          zd[i] = a * (xd[i] - yd[i]);
-      });
+      },
+      "linalg_vector_subtract_a", "=a*(x-y)", a, xd, yd, zd);
 #else
       const real_t *xp = x.data;
       const real_t *yp = y.data;

makefile

@@ -26,6 +26,7 @@ MFEM makefile targets:
    make parallel
    make debug
    make pdebug
+   make metal
    make cuda
    make hip
    make pcuda
@@ -62,6 +63,8 @@ make debug
    A shortcut to configure and build the serial debug version of the library.
 make pdebug
    A shortcut to configure and build the parallel debug version of the library.
+make metal
+   A shortcut to configure and build the serial GPU/Metal optimized version of the library.
 make cuda
    A shortcut to configure and build the serial GPU/CUDA optimized version of the library.
 make pcuda
@@ -177,7 +180,7 @@ $(call mfem-info, BLD       = $(BLD))
 
 # Include $(CONFIG_MK) unless some of the $(SKIP_INCLUDE_TARGETS) are given
 SKIP_INCLUDE_TARGETS = help config clean distclean serial parallel debug pdebug\
- cuda hip pcuda phip cudebug hipdebug pcudebug phipdebug hpc style
+ metal cuda hip pcuda phip cudebug hipdebug pcudebug phipdebug hpc style
 HAVE_SKIP_INCLUDE_TARGET = $(filter $(SKIP_INCLUDE_TARGETS),$(MAKECMDGOALS))
 ifeq (,$(HAVE_SKIP_INCLUDE_TARGET))
    $(call mfem-info, Including $(CONFIG_MK))
@@ -234,13 +237,19 @@ else
 endif
 
 # Default configuration
-ifeq ($(MFEM_USE_CUDA)$(MFEM_USE_HIP),NONO)
+ifeq ($(MFEM_USE_CUDA)$(MFEM_USE_HIP)$(MFEM_USE_METAL),NONONO)
    MFEM_CXX ?= $(HOST_CXX)
    MFEM_HOST_CXX ?= $(MFEM_CXX)
    XCOMPILER = $(CXX_XCOMPILER)
    XLINKER   = $(CXX_XLINKER)
 endif
 
+ifeq ($(MFEM_USE_METAL),YES)
+   MFEM_CXX ?= $(METAL_CXX)
+   CXXFLAGS += $(METAL_FLAGS)
+   ALL_LIBS += $(METAL_LIBS)
+endif
+
 ifeq ($(MFEM_USE_CUDA),YES)
    MFEM_CXX ?= $(CUDA_CXX)
    MFEM_HOST_CXX ?= $(HOST_CXX)
@@ -307,7 +316,7 @@ ifeq ($(MAKECMDGOALS),config)
 endif
 
 # List of MFEM dependencies, processed below
-MFEM_DEPENDENCIES = $(MFEM_REQ_LIB_DEPS) LIBUNWIND OPENMP CUDA HIP
+MFEM_DEPENDENCIES = $(MFEM_REQ_LIB_DEPS) LIBUNWIND OPENMP CUDA HIP METAL
 
 # List of deprecated MFEM dependencies, processed below
 MFEM_LEGACY_DEPENDENCIES = OPENMP
@@ -352,7 +361,7 @@ MFEM_DEFINES = MFEM_VERSION MFEM_VERSION_STRING MFEM_GIT_STRING MFEM_USE_MPI\
  MFEM_USE_SUITESPARSE MFEM_USE_GINKGO MFEM_USE_SUPERLU MFEM_USE_SUPERLU5\
  MFEM_USE_STRUMPACK MFEM_USE_GNUTLS MFEM_USE_NETCDF MFEM_USE_PETSC\
  MFEM_USE_SLEPC MFEM_USE_MPFR MFEM_USE_SIDRE MFEM_USE_FMS MFEM_USE_CONDUIT\
- MFEM_USE_PUMI MFEM_USE_HIOP MFEM_USE_GSLIB MFEM_USE_CUDA MFEM_USE_HIP\
+ MFEM_USE_PUMI MFEM_USE_HIOP MFEM_USE_GSLIB MFEM_USE_METAL MFEM_USE_CUDA MFEM_USE_HIP\
  MFEM_USE_OCCA MFEM_USE_MOONOLITH MFEM_USE_CEED MFEM_USE_RAJA MFEM_USE_UMPIRE\
  MFEM_USE_SIMD MFEM_USE_ADIOS2 MFEM_USE_MKL_CPARDISO MFEM_USE_MKL_PARDISO MFEM_USE_AMGX\
  MFEM_USE_MAGMA MFEM_USE_MUMPS MFEM_USE_ADFORWARD MFEM_USE_CODIPACK MFEM_USE_CALIPER\
@@ -445,7 +454,7 @@ OKL_DIRS = fem
 
 .PHONY: lib all clean distclean install config status info deps serial parallel	\
 	debug pdebug cuda hip pcuda cudebug pcudebug hpc style check test unittest \
-	deprecation-warnings
+	deprecation-warnings metal
 
 .SUFFIXES:
 .SUFFIXES: .cpp .o
@@ -496,9 +505,9 @@ $(BLD)libmfem.$(SO_VER): $(OBJECT_FILES)
 	   $(EXT_LIBS) -o $(@)
 
 # Shortcut targets options
-serial debug cuda hip cudebug hipdebug:           M_MPI=NO
+serial debug cuda hip cudebug hipdebug metal:     M_MPI=NO
 parallel pdebug pcuda pcudebug phip phipdebug:    M_MPI=YES
-serial parallel cuda pcuda hip phip:              M_DBG=NO
+serial parallel cuda pcuda hip phip metal:        M_DBG=NO
 debug pdebug cudebug pcudebug hipdebug phipdebug: M_DBG=YES
 cuda pcuda cudebug pcudebug:                      M_CUDA=YES
 hip phip hipdebug phipdebug:                      M_HIP=YES
@@ -508,6 +517,11 @@ serial parallel debug pdebug:
 	   $(MAKEOVERRIDES_SAVE)
 	$(MAKE) $(MAKEOVERRIDES_SAVE)
 
+metal metaldebug:
+	$(MAKE) -f $(THIS_MK) config MFEM_USE_MPI=$(M_MPI) MFEM_DEBUG=$(M_DBG) \
+	   MFEM_USE_METAL=YES $(MAKEOVERRIDES_SAVE)
+	$(MAKE) $(MAKEOVERRIDES_SAVE)
+
 cuda pcuda cudebug pcudebug:
 	$(MAKE) -f $(THIS_MK) config MFEM_USE_MPI=$(M_MPI) MFEM_DEBUG=$(M_DBG) \
 	   MFEM_USE_CUDA=$(M_CUDA) $(MAKEOVERRIDES_SAVE)
@@ -733,6 +747,7 @@ status info:
 	$(info MFEM_USE_PUMI          = $(MFEM_USE_PUMI))
 	$(info MFEM_USE_HIOP          = $(MFEM_USE_HIOP))
 	$(info MFEM_USE_GSLIB         = $(MFEM_USE_GSLIB))
+	$(info MFEM_USE_METAL         = $(MFEM_USE_METAL))
 	$(info MFEM_USE_CUDA          = $(MFEM_USE_CUDA))
 	$(info MFEM_USE_HIP           = $(MFEM_USE_HIP))
 	$(info MFEM_USE_RAJA          = $(MFEM_USE_RAJA))
@@ -788,7 +803,7 @@ FORMAT_EXCLUDE = general/tinyxml2.cpp tests/unit/catch.hpp
 FORMAT_LIST = $(filter-out $(FORMAT_EXCLUDE),$(wildcard $(FORMAT_FILES)))
 
 COUT_CERR_FILES = $(foreach dir,$(DIRS),$(dir)/*.[ch]pp)
-COUT_CERR_EXCLUDE = '^general/error\.cpp' '^general/globals\.[ch]pp'
+COUT_CERR_EXCLUDE = '^general/error\.cpp' '^general/globals\.[ch]pp' '^general/debug\.hpp'
 
 DEPRECATION_WARNING := \
 "This feature is planned for removal in the next release."\