Merge branch 'main' into users/meinersbur/flang_runtime_split-headersusers/meinersbur/flang_runtime_split-headers

8 files changed, 836 insertions, 166 deletions

diff --git a/flang/lib/Optimizer/HLFIR/IR/HLFIRDialect.cpp b/flang/lib/Optimizer/HLFIR/IR/HLFIRDialect.cpp
index d67b5fa65980..cb77aef74acd 100644
--- a/flang/lib/Optimizer/HLFIR/IR/HLFIRDialect.cpp
+++ b/flang/lib/Optimizer/HLFIR/IR/HLFIRDialect.cpp
@@ -228,3 +228,12 @@ mlir::Type hlfir::getExprType(mlir::Type variableType) {
   return hlfir::ExprType::get(variableType.getContext(), typeShape, type,
                               isPolymorphic);
 }
+
+bool hlfir::isFortranIntegerScalarOrArrayObject(mlir::Type type) {
+  if (isBoxAddressType(type))
+    return false;
+
+  mlir::Type unwrappedType = fir::unwrapPassByRefType(fir::unwrapRefType(type));
+  mlir::Type elementType = getFortranElementType(unwrappedType);
+  return mlir::isa<mlir::IntegerType>(elementType);
+}
diff --git a/flang/lib/Optimizer/HLFIR/IR/HLFIROps.cpp b/flang/lib/Optimizer/HLFIR/IR/HLFIROps.cpp
index 3a172d1b8b54..82aac7cafa1d 100644
--- a/flang/lib/Optimizer/HLFIR/IR/HLFIROps.cpp
+++ b/flang/lib/Optimizer/HLFIR/IR/HLFIROps.cpp
@@ -1342,6 +1342,108 @@ void hlfir::MatmulTransposeOp::getEffects(
 }
 
 //===----------------------------------------------------------------------===//
+// CShiftOp
+//===----------------------------------------------------------------------===//
+
+llvm::LogicalResult hlfir::CShiftOp::verify() {
+  mlir::Value array = getArray();
+  fir::SequenceType arrayTy = mlir::cast<fir::SequenceType>(
+      hlfir::getFortranElementOrSequenceType(array.getType()));
+  llvm::ArrayRef<int64_t> inShape = arrayTy.getShape();
+  std::size_t arrayRank = inShape.size();
+  mlir::Type eleTy = arrayTy.getEleTy();
+  hlfir::ExprType resultTy = mlir::cast<hlfir::ExprType>(getResult().getType());
+  llvm::ArrayRef<int64_t> resultShape = resultTy.getShape();
+  std::size_t resultRank = resultShape.size();
+  mlir::Type resultEleTy = resultTy.getEleTy();
+  mlir::Value shift = getShift();
+  mlir::Type shiftTy = hlfir::getFortranElementOrSequenceType(shift.getType());
+
+  if (eleTy != resultEleTy) {
+    if (mlir::isa<fir::CharacterType>(eleTy) &&
+        mlir::isa<fir::CharacterType>(resultEleTy)) {
+      auto eleCharTy = mlir::cast<fir::CharacterType>(eleTy);
+      auto resultCharTy = mlir::cast<fir::CharacterType>(resultEleTy);
+      if (eleCharTy.getFKind() != resultCharTy.getFKind())
+        return emitOpError("kind mismatch between input and output arrays");
+      if (eleCharTy.getLen() != fir::CharacterType::unknownLen() &&
+          resultCharTy.getLen() != fir::CharacterType::unknownLen() &&
+          eleCharTy.getLen() != resultCharTy.getLen())
+        return emitOpError(
+            "character LEN mismatch between input and output arrays");
+    } else {
+      return emitOpError(
+          "input and output arrays should have the same element type");
+    }
+  }
+
+  if (arrayRank != resultRank)
+    return emitOpError("input and output arrays should have the same rank");
+
+  constexpr int64_t unknownExtent = fir::SequenceType::getUnknownExtent();
+  for (auto [inDim, resultDim] : llvm::zip(inShape, resultShape))
+    if (inDim != unknownExtent && resultDim != unknownExtent &&
+        inDim != resultDim)
+      return emitOpError(
+          "output array's shape conflicts with the input array's shape");
+
+  int64_t dimVal = -1;
+  if (!getDim())
+    dimVal = 1;
+  else if (auto dim = fir::getIntIfConstant(getDim()))
+    dimVal = *dim;
+
+  // The DIM argument may be statically invalid (e.g. exceed the
+  // input array rank) in dead code after constant propagation,
+  // so avoid some checks unless useStrictIntrinsicVerifier is true.
+  if (useStrictIntrinsicVerifier && dimVal != -1) {
+    if (dimVal < 1)
+      return emitOpError("DIM must be >= 1");
+    if (dimVal > static_cast<int64_t>(arrayRank))
+      return emitOpError("DIM must be <= input array's rank");
+  }
+
+  if (auto shiftSeqTy = mlir::dyn_cast<fir::SequenceType>(shiftTy)) {
+    // SHIFT is an array. Verify the rank and the shape (if DIM is constant).
+    llvm::ArrayRef<int64_t> shiftShape = shiftSeqTy.getShape();
+    std::size_t shiftRank = shiftShape.size();
+    if (shiftRank != arrayRank - 1)
+      return emitOpError(
+          "SHIFT's rank must be 1 less than the input array's rank");
+
+    if (useStrictIntrinsicVerifier && dimVal != -1) {
+      // SHIFT's shape must be [d(1), d(2), ..., d(DIM-1), d(DIM+1), ..., d(n)],
+      // where [d(1), d(2), ..., d(n)] is the shape of the ARRAY.
+      int64_t arrayDimIdx = 0;
+      int64_t shiftDimIdx = 0;
+      for (auto shiftDim : shiftShape) {
+        if (arrayDimIdx == dimVal - 1)
+          ++arrayDimIdx;
+
+        if (inShape[arrayDimIdx] != unknownExtent &&
+            shiftDim != unknownExtent && inShape[arrayDimIdx] != shiftDim)
+          return emitOpError("SHAPE(ARRAY)(" + llvm::Twine(arrayDimIdx + 1) +
+                             ") must be equal to SHAPE(SHIFT)(" +
+                             llvm::Twine(shiftDimIdx + 1) +
+                             "): " + llvm::Twine(inShape[arrayDimIdx]) +
+                             " != " + llvm::Twine(shiftDim));
+        ++arrayDimIdx;
+        ++shiftDimIdx;
+      }
+    }
+  }
+
+  return mlir::success();
+}
+
+void hlfir::CShiftOp::getEffects(
+    llvm::SmallVectorImpl<
+        mlir::SideEffects::EffectInstance<mlir::MemoryEffects::Effect>>
+        &effects) {
+  getIntrinsicEffects(getOperation(), effects);
+}
+
+//===----------------------------------------------------------------------===//
 // AssociateOp
 //===----------------------------------------------------------------------===//
 
@@ -1602,6 +1704,15 @@ hlfir::ShapeOfOp::canonicalize(ShapeOfOp shapeOf,
   return llvm::LogicalResult::success();
 }
 
+mlir::OpFoldResult hlfir::ShapeOfOp::fold(FoldAdaptor adaptor) {
+  if (matchPattern(getExpr(), mlir::m_Op<hlfir::ElementalOp>())) {
+    auto elementalOp =
+        mlir::cast<hlfir::ElementalOp>(getExpr().getDefiningOp());
+    return elementalOp.getShape();
+  }
+  return {};
+}
+
 //===----------------------------------------------------------------------===//
 // GetExtent
 //===----------------------------------------------------------------------===//
diff --git a/flang/lib/Optimizer/HLFIR/Transforms/BufferizeHLFIR.cpp b/flang/lib/Optimizer/HLFIR/Transforms/BufferizeHLFIR.cpp
index 347f0a563077..48e3db0eb39b 100644
--- a/flang/lib/Optimizer/HLFIR/Transforms/BufferizeHLFIR.cpp
+++ b/flang/lib/Optimizer/HLFIR/Transforms/BufferizeHLFIR.cpp
@@ -806,8 +806,12 @@ struct ElementalOpConversion
     // elemental is a "view" over a variable (e.g parentheses or transpose).
     if (auto asExpr = elementValue.getDefiningOp<hlfir::AsExprOp>()) {
       if (asExpr->hasOneUse() && !asExpr.isMove()) {
-        elementValue = hlfir::Entity{asExpr.getVar()};
-        rewriter.eraseOp(asExpr);
+        // Check that the asExpr is the final operation before the yield,
+        // otherwise, clean-ups could impact the memory being re-used.
+        if (asExpr->getNextNode() == yield.getOperation()) {
+          elementValue = hlfir::Entity{asExpr.getVar()};
+          rewriter.eraseOp(asExpr);
+        }
       }
     }
     rewriter.eraseOp(yield);
diff --git a/flang/lib/Optimizer/HLFIR/Transforms/InlineElementals.cpp b/flang/lib/Optimizer/HLFIR/Transforms/InlineElementals.cpp
index 769e14b1316d..b68fe6ee0c74 100644
--- a/flang/lib/Optimizer/HLFIR/Transforms/InlineElementals.cpp
+++ b/flang/lib/Optimizer/HLFIR/Transforms/InlineElementals.cpp
@@ -125,7 +125,7 @@ public:
     mlir::RewritePatternSet patterns(context);
     patterns.insert<InlineElementalConversion>(context);
 
-    if (mlir::failed(mlir::applyPatternsAndFoldGreedily(
+    if (mlir::failed(mlir::applyPatternsGreedily(
             getOperation(), std::move(patterns), config))) {
       mlir::emitError(getOperation()->getLoc(),
                       "failure in HLFIR elemental inlining");
diff --git a/flang/lib/Optimizer/HLFIR/Transforms/LowerHLFIRIntrinsics.cpp b/flang/lib/Optimizer/HLFIR/Transforms/LowerHLFIRIntrinsics.cpp
index bfb4148d93bf..091ed7ed999d 100644
--- a/flang/lib/Optimizer/HLFIR/Transforms/LowerHLFIRIntrinsics.cpp
+++ b/flang/lib/Optimizer/HLFIR/Transforms/LowerHLFIRIntrinsics.cpp
@@ -463,6 +463,37 @@ struct MatmulTransposeOpConversion
   }
 };
 
+class CShiftOpConversion : public HlfirIntrinsicConversion<hlfir::CShiftOp> {
+  using HlfirIntrinsicConversion<hlfir::CShiftOp>::HlfirIntrinsicConversion;
+
+  llvm::LogicalResult
+  matchAndRewrite(hlfir::CShiftOp cshift,
+                  mlir::PatternRewriter &rewriter) const override {
+    fir::FirOpBuilder builder{rewriter, cshift.getOperation()};
+    const mlir::Location &loc = cshift->getLoc();
+
+    llvm::SmallVector<IntrinsicArgument, 3> inArgs;
+    mlir::Value array = cshift.getArray();
+    inArgs.push_back({array, array.getType()});
+    mlir::Value shift = cshift.getShift();
+    inArgs.push_back({shift, shift.getType()});
+    inArgs.push_back({cshift.getDim(), builder.getI32Type()});
+
+    auto *argLowering = fir::getIntrinsicArgumentLowering("cshift");
+    llvm::SmallVector<fir::ExtendedValue, 3> args =
+        lowerArguments(cshift, inArgs, rewriter, argLowering);
+
+    mlir::Type scalarResultType =
+        hlfir::getFortranElementType(cshift.getType());
+
+    auto [resultExv, mustBeFreed] =
+        fir::genIntrinsicCall(builder, loc, "cshift", scalarResultType, args);
+
+    processReturnValue(cshift, resultExv, mustBeFreed, builder, rewriter);
+    return mlir::success();
+  }
+};
+
 class LowerHLFIRIntrinsics
     : public hlfir::impl::LowerHLFIRIntrinsicsBase<LowerHLFIRIntrinsics> {
 public:
@@ -475,7 +506,8 @@ public:
                 AllOpConversion, AnyOpConversion, SumOpConversion,
                 ProductOpConversion, TransposeOpConversion, CountOpConversion,
                 DotProductOpConversion, MaxvalOpConversion, MinvalOpConversion,
-                MinlocOpConversion, MaxlocOpConversion>(context);
+                MinlocOpConversion, MaxlocOpConversion, CShiftOpConversion>(
+            context);
 
     // While conceptually this pass is performing dialect conversion, we use
     // pattern rewrites here instead of dialect conversion because this pass
@@ -488,8 +520,8 @@ public:
     config.enableRegionSimplification =
         mlir::GreedySimplifyRegionLevel::Disabled;
 
-    if (mlir::failed(mlir::applyPatternsAndFoldGreedily(
-            module, std::move(patterns), config))) {
+    if (mlir::failed(
+            mlir::applyPatternsGreedily(module, std::move(patterns), config))) {
       mlir::emitError(mlir::UnknownLoc::get(context),
                       "failure in HLFIR intrinsic lowering");
       signalPassFailure();
diff --git a/flang/lib/Optimizer/HLFIR/Transforms/LowerHLFIROrderedAssignments.cpp b/flang/lib/Optimizer/HLFIR/Transforms/LowerHLFIROrderedAssignments.cpp
index 424566462e8f..cba1bfc74e92 100644
--- a/flang/lib/Optimizer/HLFIR/Transforms/LowerHLFIROrderedAssignments.cpp
+++ b/flang/lib/Optimizer/HLFIR/Transforms/LowerHLFIROrderedAssignments.cpp
@@ -660,10 +660,7 @@ OrderedAssignmentRewriter::generateYieldedEntity(
     return castIfNeeded(loc, builder, {maskedValue, std::nullopt}, castToType);
   }
 
-  assert(region.hasOneBlock() && "region must contain one block");
   auto oldYield = getYield(region);
-  mlir::Block::OpListType &ops = region.back().getOperations();
-
   // Inside Forall, scalars that do not depend on forall indices can be hoisted
   // here because their evaluation is required to only call pure procedures, and
   // if they depend on a variable previously assigned to in a forall assignment,
@@ -674,24 +671,24 @@ OrderedAssignmentRewriter::generateYieldedEntity(
   bool hoistComputation = false;
   if (fir::isa_trivial(oldYield.getEntity().getType()) &&
       !constructStack.empty()) {
-    hoistComputation = true;
-    for (mlir::Operation &op : ops)
-      if (llvm::any_of(op.getOperands(), [](mlir::Value value) {
-            return isForallIndex(value);
-          })) {
-        hoistComputation = false;
-        break;
-      }
+    mlir::WalkResult walkResult =
+        region.walk([&](mlir::Operation *op) -> mlir::WalkResult {
+          if (llvm::any_of(op->getOperands(), [](mlir::Value value) {
+                return isForallIndex(value);
+              }))
+            return mlir::WalkResult::interrupt();
+          return mlir::WalkResult::advance();
+        });
+    hoistComputation = !walkResult.wasInterrupted();
   }
   auto insertionPoint = builder.saveInsertionPoint();
   if (hoistComputation)
     builder.setInsertionPoint(constructStack[0]);
 
   // Clone all operations except the final hlfir.yield.
-  assert(!ops.empty() && "yield block cannot be empty");
-  auto end = ops.end();
-  for (auto opIt = ops.begin(); std::next(opIt) != end; ++opIt)
-    (void)builder.clone(*opIt, mapper);
+  assert(region.hasOneBlock() && "region must contain one block");
+  for (auto &op : region.back().without_terminator())
+    (void)builder.clone(op, mapper);
   // Get the value for the yielded entity, it may be the result of an operation
   // that was cloned, or it may be the same as the previous value if the yield
   // operand was created before the ordered assignment tree.
diff --git a/flang/lib/Optimizer/HLFIR/Transforms/OptimizedBufferization.cpp b/flang/lib/Optimizer/HLFIR/Transforms/OptimizedBufferization.cpp
index 9327e7ad5875..bf3cf861e46f 100644
--- a/flang/lib/Optimizer/HLFIR/Transforms/OptimizedBufferization.cpp
+++ b/flang/lib/Optimizer/HLFIR/Transforms/OptimizedBufferization.cpp
@@ -159,28 +159,162 @@ containsReadOrWriteEffectOn(const mlir::MemoryEffects::EffectInstance &effect,
   return mlir::AliasResult::NoAlias;
 }
 
-// Returns true if the given array references represent identical
-// or completely disjoint array slices. The callers may use this
-// method when the alias analysis reports an alias of some kind,
-// so that we can run Fortran specific analysis on the array slices
-// to see if they are identical or disjoint. Note that the alias
-// analysis are not able to give such an answer about the references.
-static bool areIdenticalOrDisjointSlices(mlir::Value ref1, mlir::Value ref2) {
+// Helper class for analyzing two array slices represented
+// by two hlfir.designate operations.
+class ArraySectionAnalyzer {
+public:
+  // The result of the analyzis is one of the values below.
+  enum class SlicesOverlapKind {
+    // Slices overlap is unknown.
+    Unknown,
+    // Slices are definitely identical.
+    DefinitelyIdentical,
+    // Slices are definitely disjoint.
+    DefinitelyDisjoint,
+    // Slices may be either disjoint or identical,
+    // i.e. there is definitely no partial overlap.
+    EitherIdenticalOrDisjoint
+  };
+
+  // Analyzes two hlfir.designate results and returns the overlap kind.
+  // The callers may use this method when the alias analysis reports
+  // an alias of some kind, so that we can run Fortran specific analysis
+  // on the array slices to see if they are identical or disjoint.
+  // Note that the alias analysis are not able to give such an answer
+  // about the references.
+  static SlicesOverlapKind analyze(mlir::Value ref1, mlir::Value ref2);
+
+private:
+  struct SectionDesc {
+    // An array section is described by <lb, ub, stride> tuple.
+    // If the designator's subscript is not a triple, then
+    // the section descriptor is constructed as <lb, nullptr, nullptr>.
+    mlir::Value lb, ub, stride;
+
+    SectionDesc(mlir::Value lb, mlir::Value ub, mlir::Value stride)
+        : lb(lb), ub(ub), stride(stride) {
+      assert(lb && "lower bound or index must be specified");
+      normalize();
+    }
+
+    // Normalize the section descriptor:
+    //   1. If UB is nullptr, then it is set to LB.
+    //   2. If LB==UB, then stride does not matter,
+    //      so it is reset to nullptr.
+    //   3. If STRIDE==1, then it is reset to nullptr.
+    void normalize() {
+      if (!ub)
+        ub = lb;
+      if (lb == ub)
+        stride = nullptr;
+      if (stride)
+        if (auto val = fir::getIntIfConstant(stride))
+          if (*val == 1)
+            stride = nullptr;
+    }
+
+    bool operator==(const SectionDesc &other) const {
+      return lb == other.lb && ub == other.ub && stride == other.stride;
+    }
+  };
+
+  // Given an operand_iterator over the indices operands,
+  // read the subscript values and return them as SectionDesc
+  // updating the iterator. If isTriplet is true,
+  // the subscript is a triplet, and the result is <lb, ub, stride>.
+  // Otherwise, the subscript is a scalar index, and the result
+  // is <index, nullptr, nullptr>.
+  static SectionDesc readSectionDesc(mlir::Operation::operand_iterator &it,
+                                     bool isTriplet) {
+    if (isTriplet)
+      return {*it++, *it++, *it++};
+    return {*it++, nullptr, nullptr};
+  }
+
+  // Return the ordered lower and upper bounds of the section.
+  // If stride is known to be non-negative, then the ordered
+  // bounds match the <lb, ub> of the descriptor.
+  // If stride is known to be negative, then the ordered
+  // bounds are <ub, lb> of the descriptor.
+  // If stride is unknown, we cannot deduce any order,
+  // so the result is <nullptr, nullptr>
+  static std::pair<mlir::Value, mlir::Value>
+  getOrderedBounds(const SectionDesc &desc) {
+    mlir::Value stride = desc.stride;
+    // Null stride means stride=1.
+    if (!stride)
+      return {desc.lb, desc.ub};
+    // Reverse the bounds, if stride is negative.
+    if (auto val = fir::getIntIfConstant(stride)) {
+      if (*val >= 0)
+        return {desc.lb, desc.ub};
+      else
+        return {desc.ub, desc.lb};
+    }
+
+    return {nullptr, nullptr};
+  }
+
+  // Given two array sections <lb1, ub1, stride1> and
+  // <lb2, ub2, stride2>, return true only if the sections
+  // are known to be disjoint.
+  //
+  // For example, for any positive constant C:
+  //   X:Y does not overlap with (Y+C):Z
+  //   X:Y does not overlap with Z:(X-C)
+  static bool areDisjointSections(const SectionDesc &desc1,
+                                  const SectionDesc &desc2) {
+    auto [lb1, ub1] = getOrderedBounds(desc1);
+    auto [lb2, ub2] = getOrderedBounds(desc2);
+    if (!lb1 || !lb2)
+      return false;
+    // Note that this comparison must be made on the ordered bounds,
+    // otherwise 'a(x:y:1) = a(z:x-1:-1) + 1' may be incorrectly treated
+    // as not overlapping (x=2, y=10, z=9).
+    if (isLess(ub1, lb2) || isLess(ub2, lb1))
+      return true;
+    return false;
+  }
+
+  // Given two array sections <lb1, ub1, stride1> and
+  // <lb2, ub2, stride2>, return true only if the sections
+  // are known to be identical.
+  //
+  // For example:
+  //   <x, x, stride>
+  //   <x, nullptr, nullptr>
+  //
+  // These sections are identical, from the point of which array
+  // elements are being addresses, even though the shape
+  // of the array slices might be different.
+  static bool areIdenticalSections(const SectionDesc &desc1,
+                                   const SectionDesc &desc2) {
+    if (desc1 == desc2)
+      return true;
+    return false;
+  }
+
+  // Return true, if v1 is known to be less than v2.
+  static bool isLess(mlir::Value v1, mlir::Value v2);
+};
+
+ArraySectionAnalyzer::SlicesOverlapKind
+ArraySectionAnalyzer::analyze(mlir::Value ref1, mlir::Value ref2) {
   if (ref1 == ref2)
-    return true;
+    return SlicesOverlapKind::DefinitelyIdentical;
 
   auto des1 = ref1.getDefiningOp<hlfir::DesignateOp>();
   auto des2 = ref2.getDefiningOp<hlfir::DesignateOp>();
   // We only support a pair of designators right now.
   if (!des1 || !des2)
-    return false;
+    return SlicesOverlapKind::Unknown;
 
   if (des1.getMemref() != des2.getMemref()) {
     // If the bases are different, then there is unknown overlap.
     LLVM_DEBUG(llvm::dbgs() << "No identical base for:\n"
                             << des1 << "and:\n"
                             << des2 << "\n");
-    return false;
+    return SlicesOverlapKind::Unknown;
   }
 
   // Require all components of the designators to be the same.
@@ -194,104 +328,105 @@ static bool areIdenticalOrDisjointSlices(mlir::Value ref1, mlir::Value ref2) {
     LLVM_DEBUG(llvm::dbgs() << "Different designator specs for:\n"
                             << des1 << "and:\n"
                             << des2 << "\n");
-    return false;
-  }
-
-  if (des1.getIsTriplet() != des2.getIsTriplet()) {
-    LLVM_DEBUG(llvm::dbgs() << "Different sections for:\n"
-                            << des1 << "and:\n"
-                            << des2 << "\n");
-    return false;
+    return SlicesOverlapKind::Unknown;
   }
 
   // Analyze the subscripts.
-  // For example:
-  //   hlfir.designate %6#0 (%c2:%c7999:%c1, %c1:%c120:%c1, %0)  shape %9
-  //   hlfir.designate %6#0 (%c2:%c7999:%c1, %c1:%c120:%c1, %1)  shape %9
-  //
-  // If all the triplets (section speficiers) are the same, then
-  // we do not care if %0 is equal to %1 - the slices are either
-  // identical or completely disjoint.
   auto des1It = des1.getIndices().begin();
   auto des2It = des2.getIndices().begin();
   bool identicalTriplets = true;
-  for (bool isTriplet : des1.getIsTriplet()) {
-    if (isTriplet) {
-      for (int i = 0; i < 3; ++i)
-        if (*des1It++ != *des2It++) {
-          LLVM_DEBUG(llvm::dbgs() << "Triplet mismatch for:\n"
-                                  << des1 << "and:\n"
-                                  << des2 << "\n");
-          identicalTriplets = false;
-          break;
-        }
-    } else {
-      ++des1It;
-      ++des2It;
+  bool identicalIndices = true;
+  for (auto [isTriplet1, isTriplet2] :
+       llvm::zip(des1.getIsTriplet(), des2.getIsTriplet())) {
+    SectionDesc desc1 = readSectionDesc(des1It, isTriplet1);
+    SectionDesc desc2 = readSectionDesc(des2It, isTriplet2);
+
+    // See if we can prove that any of the sections do not overlap.
+    // This is mostly a Polyhedron/nf performance hack that looks for
+    // particular relations between the lower and upper bounds
+    // of the array sections, e.g. for any positive constant C:
+    //   X:Y does not overlap with (Y+C):Z
+    //   X:Y does not overlap with Z:(X-C)
+    if (areDisjointSections(desc1, desc2))
+      return SlicesOverlapKind::DefinitelyDisjoint;
+
+    if (!areIdenticalSections(desc1, desc2)) {
+      if (isTriplet1 || isTriplet2) {
+        // For example:
+        //   hlfir.designate %6#0 (%c2:%c7999:%c1, %c1:%c120:%c1, %0)
+        //   hlfir.designate %6#0 (%c2:%c7999:%c1, %c1:%c120:%c1, %1)
+        //
+        // If all the triplets (section speficiers) are the same, then
+        // we do not care if %0 is equal to %1 - the slices are either
+        // identical or completely disjoint.
+        //
+        // Also, treat these as identical sections:
+        //   hlfir.designate %6#0 (%c2:%c2:%c1)
+        //   hlfir.designate %6#0 (%c2)
+        identicalTriplets = false;
+        LLVM_DEBUG(llvm::dbgs() << "Triplet mismatch for:\n"
+                                << des1 << "and:\n"
+                                << des2 << "\n");
+      } else {
+        identicalIndices = false;
+        LLVM_DEBUG(llvm::dbgs() << "Indices mismatch for:\n"
+                                << des1 << "and:\n"
+                                << des2 << "\n");
+      }
     }
   }
-  if (identicalTriplets)
-    return true;
 
-  // See if we can prove that any of the triplets do not overlap.
-  // This is mostly a Polyhedron/nf performance hack that looks for
-  // particular relations between the lower and upper bounds
-  // of the array sections, e.g. for any positive constant C:
-  //   X:Y does not overlap with (Y+C):Z
-  //   X:Y does not overlap with Z:(X-C)
-  auto displacedByConstant = [](mlir::Value v1, mlir::Value v2) {
-    auto removeConvert = [](mlir::Value v) -> mlir::Operation * {
-      auto *op = v.getDefiningOp();
-      while (auto conv = mlir::dyn_cast_or_null<fir::ConvertOp>(op))
-        op = conv.getValue().getDefiningOp();
-      return op;
-    };
+  if (identicalTriplets) {
+    if (identicalIndices)
+      return SlicesOverlapKind::DefinitelyIdentical;
+    else
+      return SlicesOverlapKind::EitherIdenticalOrDisjoint;
+  }
 
-    auto isPositiveConstant = [](mlir::Value v) -> bool {
-      if (auto conOp =
-              mlir::dyn_cast<mlir::arith::ConstantOp>(v.getDefiningOp()))
-        if (auto iattr = mlir::dyn_cast<mlir::IntegerAttr>(conOp.getValue()))
-          return iattr.getInt() > 0;
-      return false;
-    };
+  LLVM_DEBUG(llvm::dbgs() << "Different sections for:\n"
+                          << des1 << "and:\n"
+                          << des2 << "\n");
+  return SlicesOverlapKind::Unknown;
+}
 
-    auto *op1 = removeConvert(v1);
-    auto *op2 = removeConvert(v2);
-    if (!op1 || !op2)
-      return false;
-    if (auto addi = mlir::dyn_cast<mlir::arith::AddIOp>(op2))
-      if ((addi.getLhs().getDefiningOp() == op1 &&
-           isPositiveConstant(addi.getRhs())) ||
-          (addi.getRhs().getDefiningOp() == op1 &&
-           isPositiveConstant(addi.getLhs())))
-        return true;
-    if (auto subi = mlir::dyn_cast<mlir::arith::SubIOp>(op1))
-      if (subi.getLhs().getDefiningOp() == op2 &&
-          isPositiveConstant(subi.getRhs()))
-        return true;
+bool ArraySectionAnalyzer::isLess(mlir::Value v1, mlir::Value v2) {
+  auto removeConvert = [](mlir::Value v) -> mlir::Operation * {
+    auto *op = v.getDefiningOp();
+    while (auto conv = mlir::dyn_cast_or_null<fir::ConvertOp>(op))
+      op = conv.getValue().getDefiningOp();
+    return op;
+  };
+
+  auto isPositiveConstant = [](mlir::Value v) -> bool {
+    if (auto val = fir::getIntIfConstant(v))
+      return *val > 0;
     return false;
   };
 
-  des1It = des1.getIndices().begin();
-  des2It = des2.getIndices().begin();
-  for (bool isTriplet : des1.getIsTriplet()) {
-    if (isTriplet) {
-      mlir::Value des1Lb = *des1It++;
-      mlir::Value des1Ub = *des1It++;
-      mlir::Value des2Lb = *des2It++;
-      mlir::Value des2Ub = *des2It++;
-      // Ignore strides.
-      ++des1It;
-      ++des2It;
-      if (displacedByConstant(des1Ub, des2Lb) ||
-          displacedByConstant(des2Ub, des1Lb))
-        return true;
-    } else {
-      ++des1It;
-      ++des2It;
-    }
-  }
+  auto *op1 = removeConvert(v1);
+  auto *op2 = removeConvert(v2);
+  if (!op1 || !op2)
+    return false;
 
+  // Check if they are both constants.
+  if (auto val1 = fir::getIntIfConstant(op1->getResult(0)))
+    if (auto val2 = fir::getIntIfConstant(op2->getResult(0)))
+      return *val1 < *val2;
+
+  // Handle some variable cases (C > 0):
+  //   v2 = v1 + C
+  //   v2 = C + v1
+  //   v1 = v2 - C
+  if (auto addi = mlir::dyn_cast<mlir::arith::AddIOp>(op2))
+    if ((addi.getLhs().getDefiningOp() == op1 &&
+         isPositiveConstant(addi.getRhs())) ||
+        (addi.getRhs().getDefiningOp() == op1 &&
+         isPositiveConstant(addi.getLhs())))
+      return true;
+  if (auto subi = mlir::dyn_cast<mlir::arith::SubIOp>(op1))
+    if (subi.getLhs().getDefiningOp() == op2 &&
+        isPositiveConstant(subi.getRhs()))
+      return true;
   return false;
 }
 
@@ -328,44 +463,36 @@ ElementalAssignBufferization::findMatch(hlfir::ElementalOp elemental) {
   match.array = match.assign.getLhs();
   mlir::Type arrayType = mlir::dyn_cast<fir::SequenceType>(
       fir::unwrapPassByRefType(match.array.getType()));
-  if (!arrayType)
+  if (!arrayType) {
+    LLVM_DEBUG(llvm::dbgs() << "AssignOp's result is not an array\n");
     return std::nullopt;
+  }
 
   // require that the array elements are trivial
   // TODO: this is just to make the pass easier to think about. Not an inherent
   // limitation
   mlir::Type eleTy = hlfir::getFortranElementType(arrayType);
-  if (!fir::isa_trivial(eleTy))
+  if (!fir::isa_trivial(eleTy)) {
+    LLVM_DEBUG(llvm::dbgs() << "AssignOp's data type is not trivial\n");
     return std::nullopt;
+  }
 
-  // the array must have the same shape as the elemental. CSE should have
-  // deduplicated the fir.shape operations where they are provably the same
-  // so we just have to check for the same ssa value
-  // TODO: add more ways of getting the shape of the array
-  mlir::Value arrayShape;
-  if (match.array.getDefiningOp())
-    arrayShape =
-        mlir::TypeSwitch<mlir::Operation *, mlir::Value>(
-            match.array.getDefiningOp())
-            .Case([](hlfir::DesignateOp designate) {
-              return designate.getShape();
-            })
-            .Case([](hlfir::DeclareOp declare) { return declare.getShape(); })
-            .Default([](mlir::Operation *) { return mlir::Value{}; });
-  if (!arrayShape) {
-    LLVM_DEBUG(llvm::dbgs() << "Can't get shape of " << match.array << " at "
-                            << elemental->getLoc() << "\n");
+  // The array must have the same shape as the elemental.
+  //
+  // f2018 10.2.1.2 (3) requires the lhs and rhs of an assignment to be
+  // conformable unless the lhs is an allocatable array. In HLFIR we can
+  // see this from the presence or absence of the realloc attribute on
+  // hlfir.assign. If it is not a realloc assignment, we can trust that
+  // the shapes do conform.
+  //
+  // TODO: the lhs's shape is dynamic, so it is hard to prove that
+  // there is no reallocation of the lhs due to the assignment.
+  // We can probably try generating multiple versions of the code
+  // with checking for the shape match, length parameters match, etc.
+  if (match.assign.isAllocatableAssignment()) {
+    LLVM_DEBUG(llvm::dbgs() << "AssignOp may involve (re)allocation of LHS\n");
     return std::nullopt;
   }
-  if (arrayShape != elemental.getShape()) {
-    // f2018 10.2.1.2 (3) requires the lhs and rhs of an assignment to be
-    // conformable unless the lhs is an allocatable array. In HLFIR we can
-    // see this from the presence or absence of the realloc attribute on
-    // hlfir.assign. If it is not a realloc assignment, we can trust that
-    // the shapes do conform
-    if (match.assign.getRealloc())
-      return std::nullopt;
-  }
 
   // the transformation wants to apply the elemental in a do-loop at the
   // hlfir.assign, check there are no effects which make this unsafe
@@ -419,21 +546,27 @@ ElementalAssignBufferization::findMatch(hlfir::ElementalOp elemental) {
     if (!res.isPartial()) {
       if (auto designate =
               effect.getValue().getDefiningOp<hlfir::DesignateOp>()) {
-        if (!areIdenticalOrDisjointSlices(match.array, designate.getMemref())) {
+        ArraySectionAnalyzer::SlicesOverlapKind overlap =
+            ArraySectionAnalyzer::analyze(match.array, designate.getMemref());
+        if (overlap ==
+            ArraySectionAnalyzer::SlicesOverlapKind::DefinitelyDisjoint)
+          continue;
+
+        if (overlap == ArraySectionAnalyzer::SlicesOverlapKind::Unknown) {
           LLVM_DEBUG(llvm::dbgs() << "possible read conflict: " << designate
                                   << " at " << elemental.getLoc() << "\n");
           return std::nullopt;
         }
         auto indices = designate.getIndices();
         auto elementalIndices = elemental.getIndices();
-        if (indices.size() != elementalIndices.size()) {
-          LLVM_DEBUG(llvm::dbgs() << "possible read conflict: " << designate
-                                  << " at " << elemental.getLoc() << "\n");
-          return std::nullopt;
-        }
-        if (std::equal(indices.begin(), indices.end(), elementalIndices.begin(),
+        if (indices.size() == elementalIndices.size() &&
+            std::equal(indices.begin(), indices.end(), elementalIndices.begin(),
                        elementalIndices.end()))
           continue;
+
+        LLVM_DEBUG(llvm::dbgs() << "possible read conflict: " << designate
+                                << " at " << elemental.getLoc() << "\n");
+        return std::nullopt;
       }
     }
     LLVM_DEBUG(llvm::dbgs() << "disallowed side-effect: " << effect.getValue()
@@ -1239,7 +1372,7 @@ public:
     // patterns.insert<ReductionMaskConversion<hlfir::MaxvalOp>>(context);
     // patterns.insert<ReductionMaskConversion<hlfir::MinvalOp>>(context);
 
-    if (mlir::failed(mlir::applyPatternsAndFoldGreedily(
+    if (mlir::failed(mlir::applyPatternsGreedily(
             getOperation(), std::move(patterns), config))) {
       mlir::emitError(getOperation()->getLoc(),
                       "failure in HLFIR optimized bufferization");
diff --git a/flang/lib/Optimizer/HLFIR/Transforms/SimplifyHLFIRIntrinsics.cpp b/flang/lib/Optimizer/HLFIR/Transforms/SimplifyHLFIRIntrinsics.cpp
index 60b06437e6a9..bf3d261e7e88 100644
--- a/flang/lib/Optimizer/HLFIR/Transforms/SimplifyHLFIRIntrinsics.cpp
+++ b/flang/lib/Optimizer/HLFIR/Transforms/SimplifyHLFIRIntrinsics.cpp
@@ -10,24 +10,29 @@
 // into the calling function.
 //===----------------------------------------------------------------------===//
 
+#include "flang/Optimizer/Builder/Complex.h"
 #include "flang/Optimizer/Builder/FIRBuilder.h"
 #include "flang/Optimizer/Builder/HLFIRTools.h"
+#include "flang/Optimizer/Builder/IntrinsicCall.h"
 #include "flang/Optimizer/Dialect/FIRDialect.h"
 #include "flang/Optimizer/HLFIR/HLFIRDialect.h"
 #include "flang/Optimizer/HLFIR/HLFIROps.h"
 #include "flang/Optimizer/HLFIR/Passes.h"
 #include "mlir/Dialect/Arith/IR/Arith.h"
-#include "mlir/Dialect/Func/IR/FuncOps.h"
-#include "mlir/IR/BuiltinDialect.h"
 #include "mlir/IR/Location.h"
 #include "mlir/Pass/Pass.h"
-#include "mlir/Transforms/DialectConversion.h"
+#include "mlir/Transforms/GreedyPatternRewriteDriver.h"
 
 namespace hlfir {
 #define GEN_PASS_DEF_SIMPLIFYHLFIRINTRINSICS
 #include "flang/Optimizer/HLFIR/Passes.h.inc"
 } // namespace hlfir
 
+static llvm::cl::opt<bool>
+    simplifySum("flang-simplify-hlfir-sum",
+                llvm::cl::desc("Expand hlfir.sum into an inline sequence"),
+                llvm::cl::init(true));
+
 namespace {
 
 class TransposeAsElementalConversion
@@ -38,9 +43,15 @@ public:
   llvm::LogicalResult
   matchAndRewrite(hlfir::TransposeOp transpose,
                   mlir::PatternRewriter &rewriter) const override {
+    hlfir::ExprType expr = transpose.getType();
+    // TODO: hlfir.elemental supports polymorphic data types now,
+    // so this can be supported.
+    if (expr.isPolymorphic())
+      return rewriter.notifyMatchFailure(transpose,
+                                         "TRANSPOSE of polymorphic type");
+
     mlir::Location loc = transpose.getLoc();
     fir::FirOpBuilder builder{rewriter, transpose.getOperation()};
-    hlfir::ExprType expr = transpose.getType();
     mlir::Type elementType = expr.getElementType();
     hlfir::Entity array = hlfir::Entity{transpose.getArray()};
     mlir::Value resultShape = genResultShape(loc, builder, array);
@@ -90,25 +101,398 @@ private:
   }
 };
 
+// Expand the SUM(DIM=CONSTANT) operation into .
+class SumAsElementalConversion : public mlir::OpRewritePattern<hlfir::SumOp> {
+public:
+  using mlir::OpRewritePattern<hlfir::SumOp>::OpRewritePattern;
+
+  llvm::LogicalResult
+  matchAndRewrite(hlfir::SumOp sum,
+                  mlir::PatternRewriter &rewriter) const override {
+    if (!simplifySum)
+      return rewriter.notifyMatchFailure(sum, "SUM simplification is disabled");
+
+    hlfir::Entity array = hlfir::Entity{sum.getArray()};
+    bool isTotalReduction = hlfir::Entity{sum}.getRank() == 0;
+    mlir::Value dim = sum.getDim();
+    int64_t dimVal = 0;
+    if (!isTotalReduction) {
+      // In case of partial reduction we should ignore the operations
+      // with invalid DIM values. They may appear in dead code
+      // after constant propagation.
+      auto constDim = fir::getIntIfConstant(dim);
+      if (!constDim)
+        return rewriter.notifyMatchFailure(sum, "Nonconstant DIM for SUM");
+      dimVal = *constDim;
+
+      if ((dimVal <= 0 || dimVal > array.getRank()))
+        return rewriter.notifyMatchFailure(
+            sum, "Invalid DIM for partial SUM reduction");
+    }
+
+    mlir::Location loc = sum.getLoc();
+    fir::FirOpBuilder builder{rewriter, sum.getOperation()};
+    mlir::Type elementType = hlfir::getFortranElementType(sum.getType());
+    mlir::Value mask = sum.getMask();
+
+    mlir::Value resultShape, dimExtent;
+    llvm::SmallVector<mlir::Value> arrayExtents;
+    if (isTotalReduction)
+      arrayExtents = genArrayExtents(loc, builder, array);
+    else
+      std::tie(resultShape, dimExtent) =
+          genResultShapeForPartialReduction(loc, builder, array, dimVal);
+
+    // If the mask is present and is a scalar, then we'd better load its value
+    // outside of the reduction loop making the loop unswitching easier.
+    mlir::Value isPresentPred, maskValue;
+    if (mask) {
+      if (mlir::isa<fir::BaseBoxType>(mask.getType())) {
+        // MASK represented by a box might be dynamically optional,
+        // so we have to check for its presence before accessing it.
+        isPresentPred =
+            builder.create<fir::IsPresentOp>(loc, builder.getI1Type(), mask);
+      }
+
+      if (hlfir::Entity{mask}.isScalar())
+        maskValue = genMaskValue(loc, builder, mask, isPresentPred, {});
+    }
+
+    auto genKernel = [&](mlir::Location loc, fir::FirOpBuilder &builder,
+                         mlir::ValueRange inputIndices) -> hlfir::Entity {
+      // Loop over all indices in the DIM dimension, and reduce all values.
+      // If DIM is not present, do total reduction.
+
+      // Initial value for the reduction.
+      mlir::Value reductionInitValue = genInitValue(loc, builder, elementType);
+
+      // The reduction loop may be unordered if FastMathFlags::reassoc
+      // transformations are allowed. The integer reduction is always
+      // unordered.
+      bool isUnordered = mlir::isa<mlir::IntegerType>(elementType) ||
+                         static_cast<bool>(sum.getFastmath() &
+                                           mlir::arith::FastMathFlags::reassoc);
+
+      llvm::SmallVector<mlir::Value> extents;
+      if (isTotalReduction)
+        extents = arrayExtents;
+      else
+        extents.push_back(
+            builder.createConvert(loc, builder.getIndexType(), dimExtent));
+
+      auto genBody = [&](mlir::Location loc, fir::FirOpBuilder &builder,
+                         mlir::ValueRange oneBasedIndices,
+                         mlir::ValueRange reductionArgs)
+          -> llvm::SmallVector<mlir::Value, 1> {
+        // Generate the reduction loop-nest body.
+        // The initial reduction value in the innermost loop
+        // is passed via reductionArgs[0].
+        llvm::SmallVector<mlir::Value> indices;
+        if (isTotalReduction) {
+          indices = oneBasedIndices;
+        } else {
+          indices = inputIndices;
+          indices.insert(indices.begin() + dimVal - 1, oneBasedIndices[0]);
+        }
+
+        mlir::Value reductionValue = reductionArgs[0];
+        fir::IfOp ifOp;
+        if (mask) {
+          // Make the reduction value update conditional on the value
+          // of the mask.
+          if (!maskValue) {
+            // If the mask is an array, use the elemental and the loop indices
+            // to address the proper mask element.
+            maskValue =
+                genMaskValue(loc, builder, mask, isPresentPred, indices);
+          }
+          mlir::Value isUnmasked = builder.create<fir::ConvertOp>(
+              loc, builder.getI1Type(), maskValue);
+          ifOp = builder.create<fir::IfOp>(loc, elementType, isUnmasked,
+                                           /*withElseRegion=*/true);
+          // In the 'else' block return the current reduction value.
+          builder.setInsertionPointToStart(&ifOp.getElseRegion().front());
+          builder.create<fir::ResultOp>(loc, reductionValue);
+
+          // In the 'then' block do the actual addition.
+          builder.setInsertionPointToStart(&ifOp.getThenRegion().front());
+        }
+
+        hlfir::Entity element =
+            hlfir::getElementAt(loc, builder, array, indices);
+        hlfir::Entity elementValue =
+            hlfir::loadTrivialScalar(loc, builder, element);
+        // NOTE: we can use "Kahan summation" same way as the runtime
+        // (e.g. when fast-math is not allowed), but let's start with
+        // the simple version.
+        reductionValue =
+            genScalarAdd(loc, builder, reductionValue, elementValue);
+
+        if (ifOp) {
+          builder.create<fir::ResultOp>(loc, reductionValue);
+          builder.setInsertionPointAfter(ifOp);
+          reductionValue = ifOp.getResult(0);
+        }
+
+        return {reductionValue};
+      };
+
+      llvm::SmallVector<mlir::Value, 1> reductionFinalValues =
+          hlfir::genLoopNestWithReductions(loc, builder, extents,
+                                           {reductionInitValue}, genBody,
+                                           isUnordered);
+      return hlfir::Entity{reductionFinalValues[0]};
+    };
+
+    if (isTotalReduction) {
+      hlfir::Entity result = genKernel(loc, builder, mlir::ValueRange{});
+      rewriter.replaceOp(sum, result);
+      return mlir::success();
+    }
+
+    hlfir::ElementalOp elementalOp = hlfir::genElementalOp(
+        loc, builder, elementType, resultShape, {}, genKernel,
+        /*isUnordered=*/true, /*polymorphicMold=*/nullptr,
+        sum.getResult().getType());
+
+    // it wouldn't be safe to replace block arguments with a different
+    // hlfir.expr type. Types can differ due to differing amounts of shape
+    // information
+    assert(elementalOp.getResult().getType() == sum.getResult().getType());
+
+    rewriter.replaceOp(sum, elementalOp);
+    return mlir::success();
+  }
+
+private:
+  static llvm::SmallVector<mlir::Value>
+  genArrayExtents(mlir::Location loc, fir::FirOpBuilder &builder,
+                  hlfir::Entity array) {
+    mlir::Value inShape = hlfir::genShape(loc, builder, array);
+    llvm::SmallVector<mlir::Value> inExtents =
+        hlfir::getExplicitExtentsFromShape(inShape, builder);
+    if (inShape.getUses().empty())
+      inShape.getDefiningOp()->erase();
+    return inExtents;
+  }
+
+  // Return fir.shape specifying the shape of the result
+  // of a SUM reduction with DIM=dimVal. The second return value
+  // is the extent of the DIM dimension.
+  static std::tuple<mlir::Value, mlir::Value>
+  genResultShapeForPartialReduction(mlir::Location loc,
+                                    fir::FirOpBuilder &builder,
+                                    hlfir::Entity array, int64_t dimVal) {
+    llvm::SmallVector<mlir::Value> inExtents =
+        genArrayExtents(loc, builder, array);
+    assert(dimVal > 0 && dimVal <= static_cast<int64_t>(inExtents.size()) &&
+           "DIM must be present and a positive constant not exceeding "
+           "the array's rank");
+
+    mlir::Value dimExtent = inExtents[dimVal - 1];
+    inExtents.erase(inExtents.begin() + dimVal - 1);
+    return {builder.create<fir::ShapeOp>(loc, inExtents), dimExtent};
+  }
+
+  // Generate the initial value for a SUM reduction with the given
+  // data type.
+  static mlir::Value genInitValue(mlir::Location loc,
+                                  fir::FirOpBuilder &builder,
+                                  mlir::Type elementType) {
+    if (auto ty = mlir::dyn_cast<mlir::FloatType>(elementType)) {
+      const llvm::fltSemantics &sem = ty.getFloatSemantics();
+      return builder.createRealConstant(loc, elementType,
+                                        llvm::APFloat::getZero(sem));
+    } else if (auto ty = mlir::dyn_cast<mlir::ComplexType>(elementType)) {
+      mlir::Value initValue = genInitValue(loc, builder, ty.getElementType());
+      return fir::factory::Complex{builder, loc}.createComplex(ty, initValue,
+                                                               initValue);
+    } else if (mlir::isa<mlir::IntegerType>(elementType)) {
+      return builder.createIntegerConstant(loc, elementType, 0);
+    }
+
+    llvm_unreachable("unsupported SUM reduction type");
+  }
+
+  // Generate scalar addition of the two values (of the same data type).
+  static mlir::Value genScalarAdd(mlir::Location loc,
+                                  fir::FirOpBuilder &builder,
+                                  mlir::Value value1, mlir::Value value2) {
+    mlir::Type ty = value1.getType();
+    assert(ty == value2.getType() && "reduction values' types do not match");
+    if (mlir::isa<mlir::FloatType>(ty))
+      return builder.create<mlir::arith::AddFOp>(loc, value1, value2);
+    else if (mlir::isa<mlir::ComplexType>(ty))
+      return builder.create<fir::AddcOp>(loc, value1, value2);
+    else if (mlir::isa<mlir::IntegerType>(ty))
+      return builder.create<mlir::arith::AddIOp>(loc, value1, value2);
+
+    llvm_unreachable("unsupported SUM reduction type");
+  }
+
+  static mlir::Value genMaskValue(mlir::Location loc,
+                                  fir::FirOpBuilder &builder, mlir::Value mask,
+                                  mlir::Value isPresentPred,
+                                  mlir::ValueRange indices) {
+    mlir::OpBuilder::InsertionGuard guard(builder);
+    fir::IfOp ifOp;
+    mlir::Type maskType =
+        hlfir::getFortranElementType(fir::unwrapPassByRefType(mask.getType()));
+    if (isPresentPred) {
+      ifOp = builder.create<fir::IfOp>(loc, maskType, isPresentPred,
+                                       /*withElseRegion=*/true);
+
+      // Use 'true', if the mask is not present.
+      builder.setInsertionPointToStart(&ifOp.getElseRegion().front());
+      mlir::Value trueValue = builder.createBool(loc, true);
+      trueValue = builder.createConvert(loc, maskType, trueValue);
+      builder.create<fir::ResultOp>(loc, trueValue);
+
+      // Load the mask value, if the mask is present.
+      builder.setInsertionPointToStart(&ifOp.getThenRegion().front());
+    }
+
+    hlfir::Entity maskVar{mask};
+    if (maskVar.isScalar()) {
+      if (mlir::isa<fir::BaseBoxType>(mask.getType())) {
+        // MASK may be a boxed scalar.
+        mlir::Value addr = hlfir::genVariableRawAddress(loc, builder, maskVar);
+        mask = builder.create<fir::LoadOp>(loc, hlfir::Entity{addr});
+      } else {
+        mask = hlfir::loadTrivialScalar(loc, builder, maskVar);
+      }
+    } else {
+      // Load from the mask array.
+      assert(!indices.empty() && "no indices for addressing the mask array");
+      maskVar = hlfir::getElementAt(loc, builder, maskVar, indices);
+      mask = hlfir::loadTrivialScalar(loc, builder, maskVar);
+    }
+
+    if (!isPresentPred)
+      return mask;
+
+    builder.create<fir::ResultOp>(loc, mask);
+    return ifOp.getResult(0);
+  }
+};
+
+class CShiftAsElementalConversion
+    : public mlir::OpRewritePattern<hlfir::CShiftOp> {
+public:
+  using mlir::OpRewritePattern<hlfir::CShiftOp>::OpRewritePattern;
+
+  llvm::LogicalResult
+  matchAndRewrite(hlfir::CShiftOp cshift,
+                  mlir::PatternRewriter &rewriter) const override {
+    using Fortran::common::maxRank;
+
+    hlfir::ExprType expr = mlir::dyn_cast<hlfir::ExprType>(cshift.getType());
+    assert(expr &&
+           "expected an expression type for the result of hlfir.cshift");
+    unsigned arrayRank = expr.getRank();
+    // When it is a 1D CSHIFT, we may assume that the DIM argument
+    // (whether it is present or absent) is equal to 1, otherwise,
+    // the program is illegal.
+    int64_t dimVal = 1;
+    if (arrayRank != 1)
+      if (mlir::Value dim = cshift.getDim()) {
+        auto constDim = fir::getIntIfConstant(dim);
+        if (!constDim)
+          return rewriter.notifyMatchFailure(cshift,
+                                             "Nonconstant DIM for CSHIFT");
+        dimVal = *constDim;
+      }
+
+    if (dimVal <= 0 || dimVal > arrayRank)
+      return rewriter.notifyMatchFailure(cshift, "Invalid DIM for CSHIFT");
+
+    mlir::Location loc = cshift.getLoc();
+    fir::FirOpBuilder builder{rewriter, cshift.getOperation()};
+    mlir::Type elementType = expr.getElementType();
+    hlfir::Entity array = hlfir::Entity{cshift.getArray()};
+    mlir::Value arrayShape = hlfir::genShape(loc, builder, array);
+    llvm::SmallVector<mlir::Value> arrayExtents =
+        hlfir::getExplicitExtentsFromShape(arrayShape, builder);
+    llvm::SmallVector<mlir::Value, 1> typeParams;
+    hlfir::genLengthParameters(loc, builder, array, typeParams);
+    hlfir::Entity shift = hlfir::Entity{cshift.getShift()};
+    // The new index computation involves MODULO, which is not implemented
+    // for IndexType, so use I64 instead.
+    mlir::Type calcType = builder.getI64Type();
+
+    mlir::Value one = builder.createIntegerConstant(loc, calcType, 1);
+    mlir::Value shiftVal;
+    if (shift.isScalar()) {
+      shiftVal = hlfir::loadTrivialScalar(loc, builder, shift);
+      shiftVal = builder.createConvert(loc, calcType, shiftVal);
+    }
+
+    auto genKernel = [&](mlir::Location loc, fir::FirOpBuilder &builder,
+                         mlir::ValueRange inputIndices) -> hlfir::Entity {
+      llvm::SmallVector<mlir::Value, maxRank> indices{inputIndices};
+      if (!shift.isScalar()) {
+        // When the array is not a vector, section
+        // (s(1), s(2), ..., s(dim-1), :, s(dim+1), ..., s(n)
+        // of the result has a value equal to:
+        // CSHIFT(ARRAY(s(1), s(2), ..., s(dim-1), :, s(dim+1), ..., s(n)),
+        //        SH, 1),
+        // where SH is either SHIFT (if scalar) or
+        // SHIFT(s(1), s(2), ..., s(dim-1), s(dim+1), ..., s(n)).
+        llvm::SmallVector<mlir::Value, maxRank> shiftIndices{indices};
+        shiftIndices.erase(shiftIndices.begin() + dimVal - 1);
+        hlfir::Entity shiftElement =
+            hlfir::getElementAt(loc, builder, shift, shiftIndices);
+        shiftVal = hlfir::loadTrivialScalar(loc, builder, shiftElement);
+        shiftVal = builder.createConvert(loc, calcType, shiftVal);
+      }
+
+      // Element i of the result (1-based) is element
+      // 'MODULO(i + SH - 1, SIZE(ARRAY)) + 1' (1-based) of the original
+      // ARRAY (or its section, when ARRAY is not a vector).
+      mlir::Value index =
+          builder.createConvert(loc, calcType, inputIndices[dimVal - 1]);
+      mlir::Value extent = arrayExtents[dimVal - 1];
+      mlir::Value newIndex =
+          builder.create<mlir::arith::AddIOp>(loc, index, shiftVal);
+      newIndex = builder.create<mlir::arith::SubIOp>(loc, newIndex, one);
+      newIndex = fir::IntrinsicLibrary{builder, loc}.genModulo(
+          calcType, {newIndex, builder.createConvert(loc, calcType, extent)});
+      newIndex = builder.create<mlir::arith::AddIOp>(loc, newIndex, one);
+      newIndex = builder.createConvert(loc, builder.getIndexType(), newIndex);
+
+      indices[dimVal - 1] = newIndex;
+      hlfir::Entity element = hlfir::getElementAt(loc, builder, array, indices);
+      return hlfir::loadTrivialScalar(loc, builder, element);
+    };
+
+    hlfir::ElementalOp elementalOp = hlfir::genElementalOp(
+        loc, builder, elementType, arrayShape, typeParams, genKernel,
+        /*isUnordered=*/true,
+        array.isPolymorphic() ? static_cast<mlir::Value>(array) : nullptr,
+        cshift.getResult().getType());
+    rewriter.replaceOp(cshift, elementalOp);
+    return mlir::success();
+  }
+};
+
 class SimplifyHLFIRIntrinsics
     : public hlfir::impl::SimplifyHLFIRIntrinsicsBase<SimplifyHLFIRIntrinsics> {
 public:
   void runOnOperation() override {
     mlir::MLIRContext *context = &getContext();
+
+    mlir::GreedyRewriteConfig config;
+    // Prevent the pattern driver from merging blocks
+    config.enableRegionSimplification =
+        mlir::GreedySimplifyRegionLevel::Disabled;
+
     mlir::RewritePatternSet patterns(context);
     patterns.insert<TransposeAsElementalConversion>(context);
-    mlir::ConversionTarget target(*context);
-    // don't transform transpose of polymorphic arrays (not currently supported
-    // by hlfir.elemental)
-    target.addDynamicallyLegalOp<hlfir::TransposeOp>(
-        [](hlfir::TransposeOp transpose) {
-          return mlir::cast<hlfir::ExprType>(transpose.getType())
-              .isPolymorphic();
-        });
-    target.markUnknownOpDynamicallyLegal(
-        [](mlir::Operation *) { return true; });
-    if (mlir::failed(mlir::applyFullConversion(getOperation(), target,
-                                               std::move(patterns)))) {
+    patterns.insert<SumAsElementalConversion>(context);
+    patterns.insert<CShiftAsElementalConversion>(context);
+
+    if (mlir::failed(mlir::applyPatternsGreedily(
+            getOperation(), std::move(patterns), config))) {
       mlir::emitError(getOperation()->getLoc(),
                       "failure in HLFIR intrinsic simplification");
       signalPassFailure();