Fortran: Implement the F2018 reduce intrinsic [PR85836]

2025-03-21  Paul Thomas  <pault@gcc.gnu.org>

gcc/fortran
	PR fortran/85836
	* check.cc (get_ul_from_cst_cl): New function used in
	check_operation.
	(check_operation): New function used in check_reduce and
	check_co_reduce.
	(gfc_check_co_reduce): Use it.
	(gfc_check_reduce): New function.
	(gfc_check_rename): Add prototype for intrinsic with 6 arguments.
	* gfortran.h : Add isym id for reduce and prototype for f6.
	* intrinsic.cc (do_check): Add another argument expression and use
	it in the call to the six argument specific check.
	(add_sym_6): New function.
	(add_functions): Add the discription of the reduce intrinsic and
	add it to the intrinsic list.
	* intrinsic.h : Add prototypes for gfc_check_reduce and
	gfc_resolve_reduce.
	* iresolve.cc (generate_reduce_op_wrapper): Generate a wrapper
	subroutine for the 'operation' function to enable the library
	implementation to be type agnostic and use pointer arithmetic
	throughout.
	(gfc_resolve_reduce): New function.
	* trans-expr.cc (gfc_conv_procedure_call): Add flag for scalar
	reduce. Generate a return variable 'sr' for scalar reduce, pass its
	address to the library function and return it as the scalar result.
	* trans-intrinsic.cc (gfc_conv_intrinsic_function): Array valued
	reduce is called in same way as reshape. Fall through for call to
	the scalar version.

gcc/testsuite/
	PR fortran/85836
	* gfortran.dg/reduce_1.f90: New test
	* gfortran.dg/reduce_2.f90: New test

libgfortran/
	PR libfortran/85836
	* Makefile.am : Add reduce.c
	* Makefile.in : Regenerated
	* gfortran.map : Add _gfortran_reduce, _gfortran_reduce_scalar,
	_gfortran_reduce_c and _gfortran_reduce_scalar_c to the list.
	* intrinsics/reduce.c (reduce, reduce_scalar, reduce_c,
	reduce_scalar_c): New functions and prototypes

9 files changed, 774 insertions, 51 deletions

diff --git a/gcc/fortran/check.cc b/gcc/fortran/check.cc
index 3545864..d2c8816 100644
--- a/gcc/fortran/check.cc
+++ b/gcc/fortran/check.cc
@@ -2442,31 +2442,24 @@ gfc_check_co_broadcast (gfc_expr *a, gfc_expr *source_image, gfc_expr *stat,
 }
 
 
-bool
-gfc_check_co_reduce (gfc_expr *a, gfc_expr *op, gfc_expr *result_image,
-		     gfc_expr *stat, gfc_expr *errmsg)
+/* Helper function for character arguments in gfc_check_[co_]reduce.  */
+
+static unsigned long
+get_ul_from_cst_cl (const gfc_charlen *cl)
+{
+  return cl && cl->length && cl->length->expr_type == EXPR_CONSTANT
+	 ? mpz_get_ui (cl->length->value.integer) : 0;
+};
+
+
+/* Checks shared between co_reduce and reduce.  */
+static bool
+check_operation (gfc_expr *op, gfc_expr *a, bool is_co_reduce)
 {
   symbol_attribute attr;
   gfc_formal_arglist *formal;
   gfc_symbol *sym;
 
-  if (a->ts.type == BT_CLASS)
-    {
-      gfc_error ("The A argument at %L of CO_REDUCE shall not be polymorphic",
-		 &a->where);
-      return false;
-    }
-
-  if (gfc_expr_attr (a).alloc_comp)
-    {
-      gfc_error ("Support for the A argument at %L with allocatable components"
-                 " is not yet implemented", &a->where);
-      return false;
-    }
-
-  if (!check_co_collective (a, result_image, stat, errmsg, true))
-    return false;
-
   if (!gfc_resolve_expr (op))
     return false;
 
@@ -2483,8 +2476,9 @@ gfc_check_co_reduce (gfc_expr *a, gfc_expr *op, gfc_expr *result_image,
       /* None of the intrinsics fulfills the criteria of taking two arguments,
 	 returning the same type and kind as the arguments and being permitted
 	 as actual argument.  */
-      gfc_error ("Intrinsic function %s at %L is not permitted for CO_REDUCE",
-		 op->symtree->n.sym->name, &op->where);
+      gfc_error ("Intrinsic function %s at %L is not permitted for %s",
+		 op->symtree->n.sym->name, &op->where,
+		 is_co_reduce ? "CO_REDUCE" : "REDUCE");
       return false;
     }
 
@@ -2510,12 +2504,14 @@ gfc_check_co_reduce (gfc_expr *a, gfc_expr *op, gfc_expr *result_image,
 
   if (!gfc_compare_types (&a->ts, &sym->result->ts))
     {
-      gfc_error ("The A argument at %L has type %s but the function passed as "
-		 "OPERATION at %L returns %s",
+      gfc_error ("The %s argument at %L has type %s but the function passed "
+		 "as OPERATION at %L returns %s",
+		 is_co_reduce ? "A" : "ARRAY",
 		 &a->where, gfc_typename (a), &op->where,
 		 gfc_typename (&sym->result->ts));
       return false;
     }
+
   if (!gfc_compare_types (&a->ts, &formal->sym->ts)
       || !gfc_compare_types (&a->ts, &formal->next->sym->ts))
     {
@@ -2567,42 +2563,59 @@ gfc_check_co_reduce (gfc_expr *a, gfc_expr *op, gfc_expr *result_image,
 
   if (a->ts.type == BT_CHARACTER)
     {
-      gfc_charlen *cl;
       unsigned long actual_size, formal_size1, formal_size2, result_size;
 
-      cl = a->ts.u.cl;
-      actual_size = cl && cl->length && cl->length->expr_type == EXPR_CONSTANT
-		     ? mpz_get_ui (cl->length->value.integer) : 0;
-
-      cl = formal->sym->ts.u.cl;
-      formal_size1 = cl && cl->length && cl->length->expr_type == EXPR_CONSTANT
-		     ? mpz_get_ui (cl->length->value.integer) : 0;
-
-      cl = formal->next->sym->ts.u.cl;
-      formal_size2 = cl && cl->length && cl->length->expr_type == EXPR_CONSTANT
-		     ? mpz_get_ui (cl->length->value.integer) : 0;
-
-      cl = sym->ts.u.cl;
-      result_size = cl && cl->length && cl->length->expr_type == EXPR_CONSTANT
-		    ? mpz_get_ui (cl->length->value.integer) : 0;
+      actual_size = get_ul_from_cst_cl (a->ts.u.cl);
+      formal_size1 = get_ul_from_cst_cl (formal->sym->ts.u.cl);
+      formal_size2 = get_ul_from_cst_cl (formal->next->sym->ts.u.cl);
+      result_size = get_ul_from_cst_cl (sym->ts.u.cl);
 
       if (actual_size
 	  && ((formal_size1 && actual_size != formal_size1)
 	       || (formal_size2 && actual_size != formal_size2)))
 	{
-	  gfc_error ("The character length of the A argument at %L and of the "
-		     "arguments of the OPERATION at %L shall be the same",
-		     &a->where, &op->where);
+	  gfc_error ("The character length of the %s argument at %L and of "
+		     "the arguments of the OPERATION at %L shall be the same",
+		     is_co_reduce ? "A" : "ARRAY", &a->where, &op->where);
 	  return false;
 	}
+
       if (actual_size && result_size && actual_size != result_size)
 	{
-	  gfc_error ("The character length of the A argument at %L and of the "
-		     "function result of the OPERATION at %L shall be the same",
+	  gfc_error ("The character length of the %s argument at %L and of "
+		     "the function result of the OPERATION at %L shall be the "
+		     "same", is_co_reduce ? "A" : "ARRAY",
 		     &a->where, &op->where);
 	  return false;
 	}
     }
+  return true;
+}
+
+
+bool
+gfc_check_co_reduce (gfc_expr *a, gfc_expr *op, gfc_expr *result_image,
+		     gfc_expr *stat, gfc_expr *errmsg)
+{
+  if (a->ts.type == BT_CLASS)
+    {
+      gfc_error ("The A argument at %L of CO_REDUCE shall not be polymorphic",
+		 &a->where);
+      return false;
+    }
+
+  if (gfc_expr_attr (a).alloc_comp)
+    {
+      gfc_error ("Support for the A argument at %L with allocatable components"
+		 " is not yet implemented", &a->where);
+      return false;
+    }
+
+  if (!check_co_collective (a, result_image, stat, errmsg, true))
+    return false;
+
+  if (!check_operation (op, a, true))
+    return false;
 
   return true;
 }
@@ -5136,6 +5149,62 @@ gfc_check_real (gfc_expr *a, gfc_expr *kind)
 
 
 bool
+gfc_check_reduce (gfc_expr *array, gfc_expr *operation, gfc_expr *dim,
+		  gfc_expr *mask, gfc_expr *identity, gfc_expr *ordered)
+{
+  if (array->ts.type == BT_CLASS)
+    {
+      gfc_error ("The ARRAY argument at %L of REDUCE shall not be polymorphic",
+		 &array->where);
+      return false;
+    }
+
+  if (!check_operation (operation, array, false))
+    return false;
+
+  if (dim && (dim->rank || dim->ts.type != BT_INTEGER))
+    {
+      gfc_error ("The DIM argument at %L, if present, must be an integer "
+		 "scalar", &dim->where);
+      return false;
+    }
+
+  if (mask && (array->rank != mask->rank || mask->ts.type != BT_LOGICAL))
+    {
+      gfc_error ("The MASK argument at %L, if present, must be a logical "
+		 "array with the same rank as ARRAY", &mask->where);
+      return false;
+    }
+
+  if (mask
+      && !gfc_check_conformance (array, mask,
+				 _("arguments '%s' and '%s' for intrinsic %s"),
+				 "ARRAY", "MASK", "REDUCE"))
+    return false;
+
+  if (mask && !identity)
+    gfc_warning (0, "MASK present at %L without IDENTITY", &mask->where);
+
+  if (ordered && (ordered->rank || ordered->ts.type != BT_LOGICAL))
+    {
+      gfc_error ("The ORDERED argument at %L, if present, must be a logical "
+		 "scalar", &ordered->where);
+      return false;
+    }
+
+  if (identity && (identity->rank
+		   || !gfc_compare_types (&array->ts, &identity->ts)))
+    {
+      gfc_error ("The IDENTITY argument at %L, if present, must be a scalar "
+		 "with the same type as ARRAY", &identity->where);
+      return false;
+    }
+
+  return true;
+}
+
+
+bool
 gfc_check_rename (gfc_expr *path1, gfc_expr *path2)
 {
   if (!type_check (path1, 0, BT_CHARACTER))
diff --git a/gcc/fortran/gfortran.h b/gcc/fortran/gfortran.h
index 7c6e9b6..5ef7037 100644
--- a/gcc/fortran/gfortran.h
+++ b/gcc/fortran/gfortran.h
@@ -647,6 +647,7 @@ enum gfc_isym_id
   GFC_ISYM_RANK,
   GFC_ISYM_REAL,
   GFC_ISYM_REALPART,
+  GFC_ISYM_REDUCE,
   GFC_ISYM_RENAME,
   GFC_ISYM_REPEAT,
   GFC_ISYM_RESHAPE,
@@ -2543,6 +2544,8 @@ typedef union
 	    struct gfc_expr *);
   bool (*f5)(struct gfc_expr *, struct gfc_expr *, struct gfc_expr *,
 	    struct gfc_expr *, struct gfc_expr *);
+  bool (*f6)(struct gfc_expr *, struct gfc_expr *, struct gfc_expr *,
+	    struct gfc_expr *, struct gfc_expr *, struct gfc_expr *);
 }
 gfc_check_f;
 
diff --git a/gcc/fortran/intrinsic.cc b/gcc/fortran/intrinsic.cc
index 30f532b..d2ce74f 100644
--- a/gcc/fortran/intrinsic.cc
+++ b/gcc/fortran/intrinsic.cc
@@ -331,7 +331,7 @@ do_ts29113_check (gfc_intrinsic_sym *specific, gfc_actual_arglist *arg)
 static bool
 do_check (gfc_intrinsic_sym *specific, gfc_actual_arglist *arg)
 {
-  gfc_expr *a1, *a2, *a3, *a4, *a5;
+  gfc_expr *a1, *a2, *a3, *a4, *a5, *a6;
 
   if (arg == NULL)
     return (*specific->check.f0) ();
@@ -361,6 +361,11 @@ do_check (gfc_intrinsic_sym *specific, gfc_actual_arglist *arg)
   if (arg == NULL)
     return (*specific->check.f5) (a1, a2, a3, a4, a5);
 
+  a6 = arg->expr;
+  arg = arg->next;
+  if (arg == NULL)
+    return (*specific->check.f6) (a1, a2, a3, a4, a5, a6);
+
   gfc_internal_error ("do_check(): too many args");
 }
 
@@ -838,6 +843,44 @@ add_sym_6fl (const char *name, gfc_isym_id id, enum klass cl, int actual_ok,
 }
 
 
+/* Add a symbol to the function list where the function takes
+   6 arguments.  */
+
+static void
+add_sym_6 (const char *name, gfc_isym_id id, enum klass cl, int actual_ok,
+	   bt type, int kind, int standard,
+	   bool (*check) (gfc_expr *, gfc_expr *, gfc_expr *,
+			  gfc_expr *, gfc_expr *, gfc_expr *),
+	   gfc_expr *(*simplify) (gfc_expr *, gfc_expr *, gfc_expr *,
+				  gfc_expr *, gfc_expr *, gfc_expr *),
+	   void (*resolve) (gfc_expr *, gfc_expr *, gfc_expr *, gfc_expr *,
+			    gfc_expr *, gfc_expr *, gfc_expr *),
+	   const char *a1, bt type1, int kind1, int optional1,
+	   const char *a2, bt type2, int kind2, int optional2,
+	   const char *a3, bt type3, int kind3, int optional3,
+	   const char *a4, bt type4, int kind4, int optional4,
+	   const char *a5, bt type5, int kind5, int optional5,
+	   const char *a6, bt type6, int kind6, int optional6)
+{
+  gfc_check_f cf;
+  gfc_simplify_f sf;
+  gfc_resolve_f rf;
+
+  cf.f6 = check;
+  sf.f6 = simplify;
+  rf.f6 = resolve;
+
+  add_sym (name, id, cl, actual_ok, type, kind, standard, cf, sf, rf,
+	   a1, type1, kind1, optional1, INTENT_IN,
+	   a2, type2, kind2, optional2, INTENT_IN,
+	   a3, type3, kind3, optional3, INTENT_IN,
+	   a4, type4, kind4, optional4, INTENT_IN,
+	   a5, type5, kind5, optional5, INTENT_IN,
+	   a6, type6, kind6, optional6, INTENT_IN,
+	   (void *) 0);
+}
+
+
 /* MINVAL, MAXVAL, PRODUCT, and SUM also get special treatment because
    their argument also might have to be reordered.  */
 
@@ -1358,13 +1401,13 @@ add_functions (void)
     *c_ptr_2 = "c_ptr_2", *ca = "coarray", *com = "command",
     *dist = "distance", *dm = "dim", *f = "field", *failed="failed",
     *fs = "fsource", *han = "handler", *i = "i",
-    *image = "image", *j = "j", *kind = "kind",
+    *idy = "identity", *image = "image", *j = "j", *kind = "kind",
     *l = "l", *ln = "len", *level = "level", *m = "matrix", *ma = "matrix_a",
     *mb = "matrix_b", *md = "mode", *mo = "mold", *msk = "mask",
     *n = "n", *ncopies= "ncopies", *nm = "name", *num = "number",
-    *ord = "order", *p = "p", *p1 = "path1", *p2 = "path2",
-    *pad = "pad", *pid = "pid", *pos = "pos", *pt = "pointer",
-    *r = "r", *rd = "round",
+    *op = "operation", *ord = "order", *odd = "ordered", *p = "p",
+	*p1 = "path1", *p2 = "path2", *pad = "pad", *pid = "pid", *pos = "pos",
+	*pt = "pointer", *r = "r", *rd = "round",
     *s = "s", *set = "set", *sh = "shift", *shp = "shape",
     *sig = "sig", *src = "source", *ssg = "substring",
     *sta = "string_a", *stb = "string_b", *stg = "string",
@@ -2936,6 +2979,18 @@ add_functions (void)
 
   make_generic ("sngl", GFC_ISYM_SNGL, GFC_STD_F77);
 
+  add_sym_6 ("reduce", GFC_ISYM_REDUCE, CLASS_TRANSFORMATIONAL, ACTUAL_NO,
+	     BT_REAL, dr, GFC_STD_F2018,
+	     gfc_check_reduce, NULL, gfc_resolve_reduce,
+	     ar, BT_REAL, dr, REQUIRED,
+	     op, BT_REAL, dr, REQUIRED,
+	     dm, BT_INTEGER, di, OPTIONAL,
+	     msk, BT_LOGICAL, dl, OPTIONAL,
+	     idy, BT_REAL, dr, OPTIONAL,
+	     odd, BT_LOGICAL, dl, OPTIONAL);
+
+  make_generic ("reduce", GFC_ISYM_REDUCE, GFC_STD_F2018);
+
   add_sym_2 ("rename", GFC_ISYM_RENAME, CLASS_IMPURE, ACTUAL_NO, BT_INTEGER, di,
 	     GFC_STD_GNU, gfc_check_rename, NULL, gfc_resolve_rename,
 	     p1, BT_CHARACTER, dc, REQUIRED, p2, BT_CHARACTER, dc, REQUIRED);
diff --git a/gcc/fortran/intrinsic.h b/gcc/fortran/intrinsic.h
index 34a0248..fec1c24 100644
--- a/gcc/fortran/intrinsic.h
+++ b/gcc/fortran/intrinsic.h
@@ -144,6 +144,8 @@ bool gfc_check_rand (gfc_expr *);
 bool gfc_check_range (gfc_expr *);
 bool gfc_check_rank (gfc_expr *);
 bool gfc_check_real (gfc_expr *, gfc_expr *);
+bool gfc_check_reduce (gfc_expr *, gfc_expr *, gfc_expr *, gfc_expr *,
+		       gfc_expr *, gfc_expr *);
 bool gfc_check_rename (gfc_expr *, gfc_expr *);
 bool gfc_check_repeat (gfc_expr *, gfc_expr *);
 bool gfc_check_reshape (gfc_expr *, gfc_expr *, gfc_expr *, gfc_expr *);
@@ -589,6 +591,8 @@ void gfc_resolve_parity (gfc_expr *, gfc_expr *, gfc_expr *);
 void gfc_resolve_product (gfc_expr *, gfc_expr *, gfc_expr *, gfc_expr *);
 void gfc_resolve_real (gfc_expr *, gfc_expr *, gfc_expr *);
 void gfc_resolve_realpart (gfc_expr *, gfc_expr *);
+void gfc_resolve_reduce (gfc_expr *, gfc_expr *, gfc_expr *, gfc_expr *,
+			 gfc_expr *, gfc_expr *, gfc_expr *);
 void gfc_resolve_rename (gfc_expr *, gfc_expr *, gfc_expr *);
 void gfc_resolve_repeat (gfc_expr *, gfc_expr *, gfc_expr *);
 void gfc_resolve_reshape (gfc_expr *, gfc_expr *, gfc_expr *, gfc_expr *,
diff --git a/gcc/fortran/iresolve.cc b/gcc/fortran/iresolve.cc
index 55f7e19..8189d7a 100644
--- a/gcc/fortran/iresolve.cc
+++ b/gcc/fortran/iresolve.cc
@@ -2408,6 +2408,220 @@ gfc_resolve_realpart (gfc_expr *f, gfc_expr *a)
 }
 
 
+/* Generate a wrapper subroutine for the operation so that the library REDUCE
+   function can use pointer arithmetic for OPERATION and not be dependent on
+   knowledge of its type.  */
+static gfc_symtree *
+generate_reduce_op_wrapper (gfc_expr *op)
+{
+  gfc_symbol *operation = op->symtree->n.sym;
+  gfc_symbol *wrapper, *a, *b, *c;
+  gfc_symtree *st;
+  char tname[GFC_MAX_SYMBOL_LEN+1];
+  char *name;
+  gfc_namespace *ns;
+  gfc_expr *e;
+
+  /* Find the top-level namespace.  */
+  for (ns = gfc_current_ns; ns; ns = ns->parent)
+    if (!ns->parent)
+      break;
+
+  sprintf (tname, "%s_%s", operation->name,
+	   ns->proc_name ? ns->proc_name->name : "noname");
+  name = xasprintf ("__reduce_wrapper_%s", tname);
+
+  gfc_find_sym_tree (name, ns, 0, &st);
+
+  if (st && !strcmp (name, st->name))
+    {
+      free (name);
+      return st;
+    }
+
+  /* Create the wrapper namespace and contain it in 'ns'.  */
+  gfc_namespace *sub_ns = gfc_get_namespace (ns, 0);
+  sub_ns->sibling = ns->contained;
+  ns->contained = sub_ns;
+  sub_ns->resolved = 1;
+
+  /* Set up procedure symbol.  */
+  gfc_get_symbol (name, ns, &wrapper);
+  sub_ns->proc_name = wrapper;
+  wrapper->attr.flavor = FL_PROCEDURE;
+  wrapper->attr.subroutine = 1;
+  wrapper->attr.artificial = 1;
+  wrapper->attr.if_source = IFSRC_DECL;
+  if (ns->proc_name->attr.flavor == FL_MODULE)
+      wrapper->module = ns->proc_name->name;
+  gfc_set_sym_referenced (wrapper);
+
+  /* Set up formal argument for the argument 'a'.  */
+  gfc_get_symbol ("a", sub_ns, &a);
+  a->ts = operation->ts;
+  a->attr.flavor = FL_VARIABLE;
+  a->attr.dummy = 1;
+  a->attr.artificial = 1;
+  a->attr.intent = INTENT_INOUT;
+  wrapper->formal = gfc_get_formal_arglist ();
+  wrapper->formal->sym = a;
+  gfc_set_sym_referenced (a);
+
+  /* Set up formal argument for the argument 'b'.  This is optional.  When
+     present, the wrapped function is called, otherwise 'a' is assigned
+     to 'c'.  This way, deep copies are effected in the library.  */
+  gfc_get_symbol ("b", sub_ns, &b);
+  b->ts = operation->ts;
+  b->attr.flavor = FL_VARIABLE;
+  b->attr.dummy = 1;
+  b->attr.optional= 1;
+  b->attr.artificial = 1;
+  b->attr.intent = INTENT_INOUT;
+  wrapper->formal->next = gfc_get_formal_arglist ();
+  wrapper->formal->next->sym = b;
+  gfc_set_sym_referenced (b);
+
+  /* Set up formal argument for the argument 'c'.  */
+  gfc_get_symbol ("c", sub_ns, &c);
+  c->ts = operation->ts;
+  c->attr.flavor = FL_VARIABLE;
+  c->attr.dummy = 1;
+  c->attr.artificial = 1;
+  c->attr.intent = INTENT_INOUT;
+  wrapper->formal->next->next = gfc_get_formal_arglist ();
+  wrapper->formal->next->next->sym = c;
+  gfc_set_sym_referenced (c);
+
+/* The only code is:
+		if (present (b))
+		    c = operation (a, b)
+		else
+		    c = a
+		endif
+  A call with 'b' missing provides a convenient way for the library to do
+  an intrinsic assignment instead of a call to memcpy and, where allocatable
+  components are present, a deep copy.
+
+  Code for if (present (b))  */
+  sub_ns->code = gfc_get_code (EXEC_IF);
+  gfc_code *if_block = sub_ns->code;
+  if_block->block = gfc_get_code (EXEC_IF);
+  if_block->block->expr1 = gfc_get_expr ();
+  e = if_block->block->expr1;
+  e->expr_type = EXPR_FUNCTION;
+  e->where = gfc_current_locus;
+  gfc_get_sym_tree ("present", sub_ns, &e->symtree, false);
+  e->symtree->n.sym->attr.flavor = FL_PROCEDURE;
+  e->symtree->n.sym->attr.intrinsic = 1;
+  e->ts.type = BT_LOGICAL;
+  e->ts.kind = gfc_default_logical_kind;
+  e->value.function.isym = gfc_intrinsic_function_by_id (GFC_ISYM_PRESENT);
+  e->value.function.actual = gfc_get_actual_arglist ();
+  e->value.function.actual->expr = gfc_lval_expr_from_sym (b);
+
+/* Code for c = operation (a, b)  */
+  if_block->block->next = gfc_get_code (EXEC_ASSIGN);
+  if_block->block->next->expr1 = gfc_lval_expr_from_sym (c);
+  if_block->block->next->expr2 = gfc_get_expr ();
+  e = if_block->block->next->expr2;
+  e->expr_type = EXPR_FUNCTION;
+  e->where = gfc_current_locus;
+  if_block->block->next->expr2->ts = operation->ts;
+  gfc_get_sym_tree (operation->name, ns, &e->symtree, false);
+  e->value.function.esym = if_block->block->next->expr2->symtree->n.sym;
+  e->value.function.actual = gfc_get_actual_arglist ();
+  e->value.function.actual->expr = gfc_lval_expr_from_sym (a);
+  e->value.function.actual->next = gfc_get_actual_arglist ();
+  e->value.function.actual->next->expr = gfc_lval_expr_from_sym (b);
+
+  if_block->block->block = gfc_get_code (EXEC_IF);
+  if_block->block->block->next = gfc_get_code (EXEC_ASSIGN);
+  if_block->block->block->next->expr1 = gfc_lval_expr_from_sym (c);
+  if_block->block->block->next->expr2 = gfc_lval_expr_from_sym (a);
+
+  /* It is unexpected to have some symbols added at resolution.  Commit the
+     changes in order to keep a clean state.  */
+  gfc_commit_symbol (if_block->block->expr1->symtree->n.sym);
+  gfc_commit_symbol (wrapper);
+  gfc_commit_symbol (a);
+  gfc_commit_symbol (b);
+  gfc_commit_symbol (c);
+
+  gfc_find_sym_tree (name, ns, 0, &st);
+  free (name);
+
+  return st;
+}
+
+void
+gfc_resolve_reduce (gfc_expr *f, gfc_expr *array,
+		     gfc_expr *operation,
+		     gfc_expr *dim,
+		     gfc_expr *mask,
+		     gfc_expr *identity ATTRIBUTE_UNUSED,
+		     gfc_expr *ordered ATTRIBUTE_UNUSED)
+{
+  gfc_symtree *wrapper_symtree;
+  gfc_typespec ts;
+
+  gfc_resolve_expr (array);
+  if (array->ts.type == BT_CHARACTER && array->ref)
+    gfc_resolve_substring_charlen (array);
+
+  f->ts = array->ts;
+
+  /* Replace 'operation' with its subroutine wrapper so that pointers may be
+     used throughout the library function.  */
+  wrapper_symtree = generate_reduce_op_wrapper (operation);
+  gcc_assert (wrapper_symtree && wrapper_symtree->n.sym);
+  operation->symtree = wrapper_symtree;
+  operation->ts = operation->symtree->n.sym->ts;
+
+  /* The scalar library function converts the scalar result to a dimension
+     zero descriptor and then returns the data after the call.  */
+  if (f->ts.type == BT_CHARACTER)
+    {
+      if (dim && array->rank > 1)
+	{
+	  f->value.function.name = gfc_get_string (PREFIX ("reduce_c"));
+	  f->rank = array->rank - 1;
+	}
+      else
+	{
+	  f->value.function.name = gfc_get_string (PREFIX ("reduce_scalar_c"));
+	  f->rank = 0;
+	}
+    }
+  else
+    {
+      if (dim && array->rank > 1)
+	{
+	  f->value.function.name = gfc_get_string (PREFIX ("reduce"));
+	  f->rank = array->rank - 1;
+	}
+      else
+	{
+	  f->value.function.name = gfc_get_string (PREFIX ("reduce_scalar"));
+	  f->rank = 0;
+	}
+    }
+
+  if (dim)
+    {
+      ts = dim->ts;
+      ts.kind = 4;
+      gfc_convert_type_warn (dim, &ts, 1, 0);
+    }
+
+  if (mask)
+    {
+      ts = mask->ts;
+      ts.kind = 4;
+      gfc_convert_type_warn (mask, &ts, 1, 0);
+    }
+}
+
+
 void
 gfc_resolve_rename (gfc_expr *f, gfc_expr *p1 ATTRIBUTE_UNUSED,
 		    gfc_expr *p2 ATTRIBUTE_UNUSED)
diff --git a/gcc/fortran/trans-expr.cc b/gcc/fortran/trans-expr.cc
index 923d46c..4b90b06 100644
--- a/gcc/fortran/trans-expr.cc
+++ b/gcc/fortran/trans-expr.cc
@@ -6753,6 +6753,12 @@ gfc_conv_procedure_call (gfc_se * se, gfc_symbol * sym,
   gfc_intrinsic_sym *isym = expr && expr->rank ?
 			    expr->value.function.isym : NULL;
 
+  /* In order that the library function for intrinsic REDUCE be type and kind
+     agnostic, the result is passed by reference.  Allocatable components are
+     handled within the OPERATION wrapper.  */
+  bool reduce_scalar = expr && !expr->rank && expr->value.function.isym
+		       && expr->value.function.isym->id == GFC_ISYM_REDUCE;
+
   comp = gfc_get_proc_ptr_comp (expr);
 
   bool elemental_proc = (comp
@@ -8405,6 +8411,7 @@ gfc_conv_procedure_call (gfc_se * se, gfc_symbol * sym,
   byref = (comp && (comp->attr.dimension
 	   || (comp->ts.type == BT_CHARACTER && !sym->attr.is_bind_c)))
 	   || (!comp && gfc_return_by_reference (sym));
+
   if (byref)
     {
       if (se->direct_byref)
@@ -8589,6 +8596,17 @@ gfc_conv_procedure_call (gfc_se * se, gfc_symbol * sym,
       else if (ts.type == BT_CHARACTER)
 	vec_safe_push (retargs, len);
     }
+  else if (reduce_scalar)
+    {
+      /* In order that the library function for intrinsic REDUCE be type and
+	 kind agnostic, the result is passed by reference.  Allocatable
+	 components are handled within the OPERATION wrapper.  */
+      type = gfc_typenode_for_spec (&expr->ts);
+      result = gfc_create_var (type, "sr");
+      tmp =  gfc_build_addr_expr (pvoid_type_node, result);
+      vec_safe_push (retargs, tmp);
+    }
+
   gfc_free_interface_mapping (&mapping);
 
   /* We need to glom RETARGS + ARGLIST + STRINGARGS + APPEND_ARGS.  */
@@ -8773,10 +8791,12 @@ gfc_conv_procedure_call (gfc_se * se, gfc_symbol * sym,
 
       /* Transformational functions of derived types with allocatable
 	 components must have the result allocatable components copied when the
-	 argument is actually given.  */
+	 argument is actually given.  This is unnecessry for REDUCE because the
+	 wrapper for the OPERATION function takes care of this.  */
       arg = expr->value.function.actual;
       if (result && arg && expr->rank
 	  && isym && isym->transformational
+	  && isym->id != GFC_ISYM_REDUCE
 	  && arg->expr
 	  && arg->expr->ts.type == BT_DERIVED
 	  && arg->expr->ts.u.derived->attr.alloc_comp)
@@ -8801,6 +8821,14 @@ gfc_conv_procedure_call (gfc_se * se, gfc_symbol * sym,
 	  gfc_add_expr_to_block (&se->pre, tmp);
 	}
     }
+  else if (reduce_scalar)
+    {
+      /* Even though the REDUCE intrinsic library function returns the result
+	 by reference, the scalar call passes the result as se->expr.  */
+      gfc_add_expr_to_block (&se->pre, se->expr);
+      se->expr = result;
+      gfc_add_block_to_block (&se->post, &post);
+    }
   else
     {
       /* For a function with a class array result, save the result as
diff --git a/gcc/fortran/trans-intrinsic.cc b/gcc/fortran/trans-intrinsic.cc
index 373a067..6b55017 100644
--- a/gcc/fortran/trans-intrinsic.cc
+++ b/gcc/fortran/trans-intrinsic.cc
@@ -10806,6 +10806,7 @@ gfc_conv_intrinsic_function (gfc_se * se, gfc_expr * expr)
 	    case GFC_ISYM_EOSHIFT:
 	    case GFC_ISYM_PACK:
 	    case GFC_ISYM_RESHAPE:
+	    case GFC_ISYM_REDUCE:
 	      /* For all of those the first argument specifies the type and the
 		 third is optional.  */
 	      conv_generic_with_optional_char_arg (se, expr, 1, 3);
@@ -11478,6 +11479,7 @@ gfc_conv_intrinsic_function (gfc_se * se, gfc_expr * expr)
     case GFC_ISYM_MCLOCK:
     case GFC_ISYM_MCLOCK8:
     case GFC_ISYM_RAND:
+    case GFC_ISYM_REDUCE:
     case GFC_ISYM_RENAME:
     case GFC_ISYM_SECOND:
     case GFC_ISYM_SECNDS:
@@ -11934,6 +11936,7 @@ gfc_is_intrinsic_libcall (gfc_expr * expr)
     case GFC_ISYM_FAILED_IMAGES:
     case GFC_ISYM_STOPPED_IMAGES:
     case GFC_ISYM_PACK:
+    case GFC_ISYM_REDUCE:
     case GFC_ISYM_RESHAPE:
     case GFC_ISYM_UNPACK:
       /* Pass absent optional parameters.  */
diff --git a/gcc/testsuite/gfortran.dg/reduce_1.f90 b/gcc/testsuite/gfortran.dg/reduce_1.f90
new file mode 100644
index 0000000..585cad7
--- /dev/null
+++ b/gcc/testsuite/gfortran.dg/reduce_1.f90
@@ -0,0 +1,202 @@
+! { dg-do run }
+!
+! Test results from the F2018 intrinsic REDUCE
+!
+! Contributed by Paul Thomas  <pault@gcc.gnu.org>
+!
+
+module operations
+   type :: s
+      integer, allocatable :: i
+      integer :: j
+   end type s
+
+contains
+
+   pure function add(i,j) result(sum_ij)
+      integer, intent(in) :: i, j
+      integer             :: sum_ij
+      sum_ij = i + j
+   end function add
+!
+   pure function mult(i,j) result(prod_ij)
+      integer, intent(in) :: i, j
+      integer             :: prod_ij
+      prod_ij = i * j
+   end function mult
+
+   pure function mult_by_val(i,j) result(prod_ij)
+      integer, intent(in), value :: i, j
+      integer             :: prod_ij
+      prod_ij = i * j
+   end function mult_by_val
+
+   pure function non_com(i,j) result(nc_ij)
+      integer, intent(in) :: i, j
+      integer             :: nc_ij
+      if (i > j) then
+        nc_ij = i - j
+      else
+        nc_ij = i + j
+      endif
+   end function non_com
+
+   pure function c_op (i, j) result (ij)
+      character(8), intent(in) :: i, j
+      character(8) :: ij
+      integer :: n
+      ij = i
+      do n = 1, 8
+         if (i(n:n) .ne. j(n:n)) ij(n:n) = '!'
+      end do
+   end function c_op
+
+   pure function t_op (i, j) result (ij)
+      type(s), intent(in) :: i, j
+      type(s) :: ij
+      ij%i = non_com (i%i, j%i)
+      ij%j = non_com (j%j, i%j)
+   end function t_op
+
+   pure function t_add (i, j) result (ij)
+      type(s), intent(in) :: i, j
+      type(s) :: ij
+      ij%i = i%i + j%i
+      ij%j = j%j + i%j
+   end function t_add
+end module operations
+
+program test_reduce
+   use operations
+   implicit none
+   integer :: i
+   integer, parameter :: n = 3
+   integer, parameter :: vec(n) = [2, 5, 10]
+   integer, parameter :: mat(n,2) = reshape([vec,2*vec],shape=[size(vec),2])
+   integer :: res0
+   integer, dimension(:), allocatable :: res1
+   integer, dimension(:,:), allocatable :: res2
+   logical, parameter :: t = .true., f = .false.
+   LOGICAL, PARAMETER :: keep(n) = [t,f,t]
+   logical, parameter :: keepM(n,2) = reshape([keep,keep],shape=[n,2])
+   logical, parameter :: all_false(n,2) = reshape ([(f, i = 1,2*n)],[n,2])
+   character(*), parameter :: carray (4) = ['abctefgh', 'atcdefgh', &
+                                            'abcdefth', 'abcdtfgh']
+   character(:), allocatable :: cres0, cres1(:)
+   type(s), allocatable :: tres1(:)
+   type(s), allocatable :: tres2(:,:)
+   type(s) :: tres2_na(2, 4)
+   type(s), allocatable :: tarray(:,:,:)
+   type(s), allocatable :: tvec(:)
+   type(s), allocatable :: tres0
+   integer, allocatable :: ires(:)
+
+! Simple cases with and without DIM
+   res0 = reduce (vec, add, dim=1)
+   if (res0 /= 17) stop 1
+   res0 = reduce (vec, mult, 1)
+   if (res0 /= 100) stop 2
+   res1 = reduce (mat, add, 1)
+   if (any (res1 /= [17, 34])) stop 3
+   res1 = reduce (mat, mult, 1)
+   if (any (res1 /= [100, 800])) stop 4
+   res1 = reduce (mat, add, 2)
+   if (any (res1 /= [6, 15, 30])) stop 5
+   res1 = reduce (mat, mult, 2)
+   if (any (res1 /= [8, 50, 200])) stop 6
+   res0 = reduce (mat, add)
+   if (res0 /= 51) stop 7
+   res0 = reduce (mat, mult)
+   if (res0 /= 80000) stop 8
+! Repeat previous test with arguments passed by value to operation
+   res0 = reduce (mat, mult_by_val)
+   if (res0 /= 80000) stop 9
+
+! Using MASK and IDENTITY
+   res0 = reduce (vec,add, mask=keep, identity = 1)
+   if (res0 /= 12) stop 10
+   res0 = reduce (vec,mult, mask=keep, identity = 1)
+   if (res0 /= 20) stop 11
+   res0 = reduce (mat, add, mask=keepM, identity = 1)
+   if (res0 /= 36) stop 12
+   res0 = reduce (mat, mult, mask=keepM, identity = 1)
+   if (res0 /= 1600) stop 13
+   res0 = reduce (mat, mult, mask=all_false, identity = -1)
+   if (res0 /= -1) stop 14
+
+! 3-D ARRAYs with and without DIM and MASK
+   res0 = reduce (reshape ([(i, i=1,8)], [2,2,2]),mult)
+   if (res0 /= 40320) stop 15
+   res2 = reduce (reshape ([(i, i=1,8)], [2,2,2]),mult,dim=2)
+   if (any (res2 /= reshape ([3,8,35,48], [2,2]))) stop 16
+   res2 = reduce (reshape ([(i, i=1,8)], [2,2,2]),mult,dim=2, &
+                  mask=reshape ([t,f,t,f,t,f,t,f],[2,2,2]), identity=-1)
+   if (any (res2 /= reshape ([3,-1,35,-1], [2,2]))) stop 17
+   res2 = reduce (reshape([(i, i=1,16)], [2,4,2]), add, dim = 3, &
+                  mask=reshape([f,t,t,f,t,t,t,t,t,t,t,t,t,t,t,t],[2,4,2]), &
+                  identity=-1)
+   if (any (res2 /= reshape ([9,12,14,12,18,20,22,24], [2,4]))) stop 18
+   res1 = reduce (reshape([(i, i=1,16)], [4,4]),add, dim = 2, &
+                  mask=reshape([f,t,t,f,t,t,t,t,t,t,t,t,t,t,t,t],[4,4]), &
+                  identity=-1)
+   if (any (res1 /= [27,32,36,36])) stop 19
+
+! Verify that the library function treats non-comutative OPERATION in the
+! correct order. If this were incorrect,the result would be [9,8,8,12,8,8,8,8].
+   res2 = reduce (reshape([(i, i=1,16)], [2,4,2]), non_com, dim = 3, &
+                  mask=reshape([f,t,t,f,t,t,t,t,t,t,t,t,t,t,t,t],[2,4,2]), &
+                  identity=-1)
+   if (any (res2 /= reshape([9,12,14,12,18,20,22,24],shape(res2)))) stop 20
+
+! Character ARRAY and OPERATION
+   cres0 = reduce (carray, c_op); if (cres0 /= 'a!c!!f!h') stop 21
+   cres1 = reduce (reshape (carray, [2,2]), c_op, dim = 1)
+   if (any (cres1 /= ['a!c!efgh','abcd!f!h'])) stop 22
+
+! Derived type ARRAY and OPERATION - was checked for memory leaks of the
+! allocatable component.
+! tarray = reshape([(s(i, i), i = 1, 16)], [2,4,2]) leaks memory!
+   allocate (tvec(16))
+   do i = 1, 16
+     tvec(i)%i = i
+     tvec(i)%j = i
+   enddo
+   tarray = reshape(tvec, [2,4,2])
+
+   tres2 = reduce (tarray, t_op, dim = 3, &
+                   mask=reshape([t,t,t,f,t,t,t,t,t,f,t,t,t,t,t,t],[2,4,2]), &
+                   identity = s(NULL(),1))
+   ires = [10,2,14,12,18,20,22,24]
+   tres1 = reshape (tres2, [size (tres2, 1)* size (tres2, 2)])
+   do i = 1, size (tres2, 1)* size (tres2, 2)
+      if (tres1(i)%i /= ires(i)) stop 23
+   end do
+   if (any (tres2%j /= reshape([8,2,8,12,8,8,8,8],shape(tres2)))) stop 24
+
+! Check that the non-allocatable result with an allocatable component does not
+! leak memory from the allocatable component
+   tres2_na = reduce (tarray, t_op, dim = 3, &
+                      mask=reshape([t,t,t,f,t,t,t,t,t,f,t,t,t,t,t,t],[2,4,2]), &
+                      identity = s(NULL(),1))
+   tres1 = reshape (tres2_na, [size (tres2_na, 1)* size (tres2, 2)])
+   do i = 1, size (tres2_na, 1)* size (tres2_na, 2)
+      if (tres1(i)%i /= ires(i)) stop 25
+   end do
+   if (any (tres2_na%j /= reshape([8,2,8,12,8,8,8,8],shape(tres2_na)))) stop 26
+
+
+   tres0 = reduce (tarray, t_add)
+   if (tres0%i /= 136) stop 27
+   if (tres0%j /= 136) stop 28
+
+! Test array being a component of an array of derived types
+   i = reduce (tarray%j, add, &
+               mask=reshape([t,t,t,f,t,t,t,t,t,f,t,t,t,t,f,t],[2,4,2]), &
+               identity = 0)
+   if (i /= 107) stop 29
+
+
+! Deallocate the allocatable components and then the allocatable variables
+   tres2_na = reshape ([(s(NULL (), 0), i = 1, size (tres2_na))], shape (tres2_na))
+   deallocate (res1, res2, cres0, cres1, tarray, ires, tres0, tres1, tres2, tvec)
+end
diff --git a/gcc/testsuite/gfortran.dg/reduce_2.f90 b/gcc/testsuite/gfortran.dg/reduce_2.f90
new file mode 100644
index 0000000..52d7c68
--- /dev/null
+++ b/gcc/testsuite/gfortran.dg/reduce_2.f90
@@ -0,0 +1,145 @@
+! { dg-do compile }
+!
+! Test argument compliance for the F2018 intrinsic REDUCE
+!
+! Contributed by Paul Thomas  <pault@gcc.gnu.org>
+!
+  class (*), allocatable :: cstar (:)
+  integer, allocatable :: i(:,:,:)
+  integer :: n(2,2)
+  Logical :: l1(4), l2(2,3), l3(2,2)
+
+! The ARRAY argument at (1) of REDUCE shall not be polymorphic
+  print *, reduce (cstar, add) ! { dg-error "shall not be polymorphic" }
+
+! OPERATION argument at %L must be a PURE function
+  print *, reduce (i, iadd) ! { dg-error "must be a PURE function" }
+  print *, reduce (i, foo) ! { dg-error "must be a PURE function" }
+
+! The function passed as OPERATION at (1) shall have scalar nonallocatable
+! nonpointer arguments and return a nonallocatable nonpointer scalar
+  print *, reduce (i, vadd) ! { dg-error "return a nonallocatable nonpointer scalar" }
+
+! The function passed as OPERATION at (1) shall have two arguments
+  print *, reduce (i, add_1a) ! { dg-error "shall have two arguments" }
+  print *, reduce (i, add_3a) ! { dg-error "shall have two arguments" }
+
+!The ARRAY argument at (1) has type INTEGER(4) but the function passed as OPERATION at
+! (2) returns REAL(4)
+  print *, reduce (i, add_r) ! { dg-error "returns REAL" }
+
+! The function passed as OPERATION at (1) shall have scalar nonallocatable nonpointer
+! arguments and return a nonallocatable nonpointer scalar
+  print *, reduce (i, add_a) ! { dg-error "return a nonallocatable nonpointer scalar" }
+
+! The function passed as OPERATION at (1) shall have scalar nonallocatable nonpointer arguments and
+! return a nonallocatable nonpointer scalar
+  print *, reduce (i, add_array) ! { dg-error "scalar nonallocatable nonpointer arguments" }
+
+! The function passed as OPERATION at (1) shall not have the OPTIONAL attribute for either of the arguments
+  print *, reduce (i, add_optional) ! { dg-error "shall not have the OPTIONAL attribute" }
+
+! The function passed as OPERATION at (1) shall have the VALUE attribute either for none or both arguments
+  print *, reduce (i, add_one_value) ! { dg-error "VALUE attribute either for none or both arguments" }
+
+! The character length of the ARRAY argument at (1) and of the arguments of the OPERATION at (2)
+! shall be the same
+  print *, reduce ([character(4) :: 'abcd','efgh'], char_one) ! { dg-error "The character length of the ARRAY" }
+
+! The character length of the ARRAY argument at (1) and of the function result of the OPERATION
+! at (2) shall be the same
+  print *, reduce ([character(4) :: 'abcd','efgh'], char_two) ! { dg-error "function result of the OPERATION" }
+
+! The character length of the ARRAY argument at (1) and of the arguments of the OPERATION at
+! (2) shall be the same
+  print *, reduce ([character(4) :: 'abcd','efgh'], char_three) ! { dg-error "arguments of the OPERATION" }
+
+! The DIM argument at (1), if present, must be an integer scalar
+  print *, reduce (i, add, dim = 2.0) ! { dg-error "must be an integer scalar" }
+
+! The DIM argument at (1), if present, must be an integer scalar
+  print *, reduce (i, add, dim = [2]) ! { dg-error "must be an integer scalar" }
+
+! The MASK argument at (1), if present, must be a logical array with the same rank as ARRAY
+  print *, reduce (n, add, mask = l1) ! { dg-error "same rank as ARRAY" }
+  print *, reduce (n, add, mask = n) ! { dg-error "must be a logical array" }
+
+! Different shape for arguments 'ARRAY' and 'MASK' for intrinsic REDUCE at (1) on
+! dimension 2 (2 and 3)
+  print *, reduce (n, add, mask = l2) ! { dg-error "Different shape" }
+
+! The IDENTITY argument at (1), if present, must be a scalar with the same type as ARRAY
+  print *, reduce (n, add, mask = l3, identity = 1.0) ! { dg-error "same type as ARRAY" }
+  print *, reduce (n, add, mask = l3, identity = [1]) ! { dg-error "must be a scalar" }
+
+! MASK present at (1) without IDENTITY
+  print *, reduce (n, add, mask = l3) ! { dg-warning "without IDENTITY" }
+
+contains
+  pure function add(i,j) result(sum_ij)
+    integer, intent(in) :: i, j
+    integer             :: sum_ij
+    sum_ij = i + j
+  end function add
+  function iadd(i,j) result(sum_ij)
+    integer, intent(in) :: i, j
+    integer             :: sum_ij
+    sum_ij = i + j
+  end function iadd
+  pure function vadd(i,j) result(sum_ij)
+    integer, intent(in) :: i, j
+    integer             :: sum_ij(6)
+    sum_ij = i + j
+  end function vadd
+  pure function add_1a(i) result(sum_ij)
+    integer, intent(in) :: i
+    integer             :: sum_ij
+    sum_ij = 0
+  end function add_1a
+  pure function add_3a(i) result(sum_ij)
+    integer, intent(in) :: i
+    integer             :: sum_ij
+    sum_ij = 0
+  end function add_3a
+  pure function add_r(i, j) result(sum_ij)
+    integer, intent(in) :: i, j
+    real             :: sum_ij
+    sum_ij = 0.0
+  end function add_r
+  pure function add_a(i, j) result(sum_ij)
+    integer, intent(in) :: i, j
+    integer, allocatable :: sum_ij
+    sum_ij = 0
+  end function add_a
+  pure function add_array(i, j) result(sum_ij)
+    integer, intent(in), dimension(:) :: i, j
+    integer :: sum_ij
+    sum_ij = 0
+  end function add_array
+  pure function add_optional(i, j) result(sum_ij)
+    integer, intent(in), optional :: i, j
+    integer :: sum_ij
+    sum_ij = 0
+  end function add_optional
+  pure function add_one_value(i, j) result(sum_ij)
+    integer, intent(in), value :: i
+    integer, intent(in) :: j
+    integer :: sum_ij
+    sum_ij = 0
+  end function add_one_value
+  pure function char_one(i, j) result(sum_ij)
+    character(8), intent(in) :: i, j
+    character(8) :: sum_ij
+  end function char_one
+  pure function char_two(i, j) result(sum_ij)
+    character(4), intent(in) :: i, j
+    character(8) :: sum_ij
+  end function char_two
+  pure function char_three(i, j) result(sum_ij)
+    character(8), intent(in) :: i
+    character(4), intent(in) :: j
+    character(4) :: sum_ij
+  end function char_three
+  subroutine foo
+  end subroutine foo
+end