diff --git a/flang/lib/Evaluate/fold-real.cpp b/flang/lib/Evaluate/fold-real.cpp
--- a/flang/lib/Evaluate/fold-real.cpp
+++ b/flang/lib/Evaluate/fold-real.cpp
@@ -294,6 +294,8 @@
     return FoldSum<T>(context, std::move(funcRef));
   } else if (name == "tiny") {
     return Expr<T>{Scalar<T>::TINY()};
+  } else if (name == "__builtin_fma") {
+    CHECK(args.size() == 3);
   } else if (name == "__builtin_ieee_next_after") {
     if (const auto *yExpr{UnwrapExpr<Expr<SomeReal>>(args[1])}) {
       return common::visit(
diff --git a/flang/lib/Evaluate/intrinsics.cpp b/flang/lib/Evaluate/intrinsics.cpp
--- a/flang/lib/Evaluate/intrinsics.cpp
+++ b/flang/lib/Evaluate/intrinsics.cpp
@@ -880,6 +880,8 @@
             {"back", AnyLogical, Rank::elemental, Optionality::optional},
             DefaultingKIND},
         KINDInt},
+    {"__builtin_fma", {{"f1", SameReal}, {"f2", SameReal}, {"f3", SameReal}},
+        SameReal},
     {"__builtin_ieee_is_nan", {{"a", AnyFloating}}, DefaultLogical},
     {"__builtin_ieee_is_negative", {{"a", AnyFloating}}, DefaultLogical},
     {"__builtin_ieee_is_normal", {{"a", AnyFloating}}, DefaultLogical},
diff --git a/flang/lib/Optimizer/Builder/IntrinsicCall.cpp b/flang/lib/Optimizer/Builder/IntrinsicCall.cpp
--- a/flang/lib/Optimizer/Builder/IntrinsicCall.cpp
+++ b/flang/lib/Optimizer/Builder/IntrinsicCall.cpp
@@ -1023,6 +1023,10 @@
     // math::FloorOp returns a real, while Fortran FLOOR returns integer.
     {"floor", "floorf", genF32F32FuncType, genMathOp<mlir::math::FloorOp>},
     {"floor", "floor", genF64F64FuncType, genMathOp<mlir::math::FloorOp>},
+    {"fma", "llvm.fma.f32", genF32F32F32F32FuncType,
+     genMathOp<mlir::math::FmaOp>},
+    {"fma", "llvm.fma.f64", genF64F64F64F64FuncType,
+     genMathOp<mlir::math::FmaOp>},
     {"gamma", "tgammaf", genF32F32FuncType, genLibCall},
     {"gamma", "tgamma", genF64F64FuncType, genLibCall},
     {"hypot", "hypotf", genF32F32F32FuncType, genLibCall},
diff --git a/flang/module/__fortran_builtins.f90 b/flang/module/__fortran_builtins.f90
--- a/flang/module/__fortran_builtins.f90
+++ b/flang/module/__fortran_builtins.f90
@@ -52,6 +52,7 @@
     __builtin_threadIdx, __builtin_blockDim, __builtin_blockIdx, __builtin_gridDim
   integer, parameter :: __builtin_warpsize = 32
 
+  intrinsic :: __builtin_fma
   intrinsic :: __builtin_ieee_is_nan, __builtin_ieee_is_negative, &
     __builtin_ieee_is_normal
   intrinsic :: __builtin_ieee_next_after, __builtin_ieee_next_down, &
diff --git a/flang/module/ieee_arithmetic.f90 b/flang/module/ieee_arithmetic.f90
--- a/flang/module/ieee_arithmetic.f90
+++ b/flang/module/ieee_arithmetic.f90
@@ -15,6 +15,7 @@
   use __Fortran_ieee_exceptions
 
   use __Fortran_builtins, only: &
+    ieee_fma => __builtin_fma, &
     ieee_is_nan => __builtin_ieee_is_nan, &
     ieee_is_negative => __builtin_ieee_is_negative, &
     ieee_is_normal => __builtin_ieee_is_normal, &
@@ -220,16 +221,6 @@
   PRIVATE_RR(IEEE_COPY_SIGN)
 #undef IEEE_COPY_SIGN_RR
 
-#define IEEE_FMA_R(AKIND) \
-  elemental real(AKIND) function ieee_fma_a##AKIND(a, b, c); \
-    real(AKIND), intent(in) :: a, b, c; \
-  end function ieee_fma_a##AKIND;
-  interface ieee_fma
-    SPECIFICS_R(IEEE_FMA_R)
-  end interface ieee_fma
-  PRIVATE_R(IEEE_FMA)
-#undef IEEE_FMA_R
-
 #define IEEE_GET_ROUNDING_MODE_I(RKIND) \
   subroutine ieee_get_rounding_mode_i##RKIND(round_value, radix); \
     import ieee_round_type; \
diff --git a/flang/test/Lower/Intrinsics/fma.f90 b/flang/test/Lower/Intrinsics/fma.f90
new file mode 100644
--- /dev/null
+++ b/flang/test/Lower/Intrinsics/fma.f90
@@ -0,0 +1,22 @@
+! RUN: bbc -emit-fir %s -o - | FileCheck --check-prefix=CHECK-FIR %s
+! RUN: %flang_fc1 -emit-llvm %s -o - | FileCheck --check-prefix=CHECK-LLVMIR %s
+
+function test_real4(a, x, y)
+  use ieee_arithmetic, only: ieee_fma
+  real :: a, x, y
+  test_real4 = ieee_fma(a, x, y)
+end function
+
+! CHECK-LABEL: @_QPtest_real4
+! CHECK-FIR: {{%[A-Za-z0-9._]+}} = math.fma {{%[0-9]+}}, {{%[0-9]+}}, {{%[0-9]+}} {{.*}} : f32
+! CHECK-LLVMIR: {{%[A-Za-z0-9._]+}} = call {{.*}} float @llvm.fma.f32(float {{%[0-9]+}}, float {{%[0-9]+}}, float {{%[0-9]+}})
+
+function test_real8(a, x, y)
+  use ieee_arithmetic, only: ieee_fma
+  real(8) :: a, x, y
+  test_real8 = ieee_fma(a, x, y)
+end function
+
+! CHECK-LABEL: @_QPtest_real8
+! CHECK-FIR: {{%[A-Za-z0-9._]+}} = math.fma {{%[0-9]+}}, {{%[0-9]+}}, {{%[0-9]+}} {{.*}} : f64
+! CHECK-LLVMIR: {{%[A-Za-z0-9._]+}} = call {{.*}} double @llvm.fma.f64(double {{%[0-9]+}}, double {{%[0-9]+}}, double {{%[0-9]+}})