diff --git a/flang/lib/Lower/ConvertCall.cpp b/flang/lib/Lower/ConvertCall.cpp
--- a/flang/lib/Lower/ConvertCall.cpp
+++ b/flang/lib/Lower/ConvertCall.cpp
@@ -1281,6 +1281,13 @@
                                                       *callContext.resultType);
     return {hlfir::EntityWithAttributes{sumOp.getResult()}};
   }
+  if (intrinsic.name == "matmul") {
+    llvm::SmallVector<mlir::Value> operands = getOperandVector(loweredActuals);
+    hlfir::MatmulOp matmulOp = builder.create<hlfir::MatmulOp>(
+        loc, operands[0], operands[1], *callContext.resultType);
+
+    return {hlfir::EntityWithAttributes{matmulOp.getResult()}};
+  }
 
   // TODO add hlfir operations for other transformational intrinsics here
 
diff --git a/flang/lib/Optimizer/HLFIR/IR/HLFIROps.cpp b/flang/lib/Optimizer/HLFIR/IR/HLFIROps.cpp
--- a/flang/lib/Optimizer/HLFIR/IR/HLFIROps.cpp
+++ b/flang/lib/Optimizer/HLFIR/IR/HLFIROps.cpp
@@ -580,6 +580,24 @@
   return mlir::success();
 }
 
+void hlfir::MatmulOp::build(mlir::OpBuilder &builder,
+                            mlir::OperationState &result, mlir::Value lhs,
+                            mlir::Value rhs, mlir::Type stmtResultType) {
+  assert(lhs && rhs && stmtResultType && "No arguments are optional");
+
+  // figure out result type
+  fir::SequenceType firResultTy =
+      hlfir::getFortranElementOrSequenceType(stmtResultType)
+          .dyn_cast<fir::SequenceType>();
+  assert(firResultTy && "result must be of array type");
+  mlir::Type scalarType = firResultTy.getEleTy();
+  llvm::ArrayRef<int64_t> shape = firResultTy.getShape();
+  mlir::Type resultType = hlfir::ExprType::get(
+      builder.getContext(), shape, scalarType, /*polymorphic=*/false);
+
+  build(builder, result, resultType, lhs, rhs);
+}
+
 //===----------------------------------------------------------------------===//
 // AssociateOp
 //===----------------------------------------------------------------------===//
diff --git a/flang/test/Lower/HLFIR/matmul.f90 b/flang/test/Lower/HLFIR/matmul.f90
new file mode 100644
--- /dev/null
+++ b/flang/test/Lower/HLFIR/matmul.f90
@@ -0,0 +1,19 @@
+! Test lowering of MATMUL intrinsic to HLFIR
+! RUN: bbc -emit-fir -hlfir -o - %s 2>&1 | FileCheck %s
+
+subroutine matmul1(lhs, rhs, res)
+  integer :: lhs(:,:), rhs(:,:), res(:,:)
+  res = MATMUL(lhs, rhs)
+endsubroutine
+! CHECK-LABEL: func.func @_QPmatmul1
+! CHECK:           %[[LHS:.*]]: !fir.box<!fir.array<?x?xi32>> {fir.bindc_name = "lhs"}
+! CHECK:           %[[RHS:.*]]: !fir.box<!fir.array<?x?xi32>> {fir.bindc_name = "rhs"}
+! CHECK:           %[[RES:.*]]: !fir.box<!fir.array<?x?xi32>> {fir.bindc_name = "res"}
+! CHECK-DAG:     %[[LHS_VAR:.*]]:2 = hlfir.declare %[[LHS]]
+! CHECK-DAG:     %[[RHS_VAR:.*]]:2 = hlfir.declare %[[RHS]]
+! CHECK-DAG:     %[[RES_VAR:.*]]:2 = hlfir.declare %[[RES]]
+! CHECK-NEXT:    %[[EXPR:.*]] = hlfir.matmul %[[LHS_VAR]]#0 %[[RHS_VAR]]#0 {fastmath = #arith.fastmath<contract>} : (!fir.box<!fir.array<?x?xi32>>, !fir.box<!fir.array<?x?xi32>>) -> !hlfir.expr<?x?xi32>
+! CHECK-NEXT:    hlfir.assign %[[EXPR]] to %[[RES_VAR]]#0 : !hlfir.expr<?x?xi32>, !fir.box<!fir.array<?x?xi32>>
+! CHECK-NEXT:    hlfir.destroy %[[EXPR]]
+! CHECK-NEXT:    return
+! CHECK-NEXT:   }