Hi All,

This version of the REDUCE intrinsic patch has evolved somewhat since the
posting on 2nd March. The most important changes are to the wrapper
function and the addition of two testsuite entries.

The wrapper function now effects:
    subroutine wrapper (a, b, c)
         type_of_ARRAY, intent(inout) :: a, c
         type_of_ARRAY, intent(inout), optional :: b
         if (present (b)) then
            c = OPERATION (a,b )
         else
            c = a
         end if
    end subroutine

The reason for wrapping OPERATION in a subroutine is to allow pointer
arithmetic to be used throughout in the library function. The only thing
that needs to be known about the type and kind of ARRAY is the element
size. The second branch in the wrapper allows deep copies to be done in the
library function, such that derived types with allocatable components do
not leak memory. This is needed at the final step of the algorithm to copy
the result from each iteration to the result and then nullify it.

This is undoubtedly a bit heavy going for intrinsic types and so, one day
soon I will possibly do a bit of M4ery. That said, the present version
works for all types of ARRAY and I worry a bit about how much this
intrinsic will be used. Thoughts?

A slight niggle is the linker error that comes up if compiled without any
optimization:
/usr/bin/ld: warning: /tmp/cc9cx9Rw.o: requires executable stack (because
the .note.GNU-stack section is executable)
I think that this is unlikely to present a security issue, however, since
it disappears at -O1, I went through each of the options triggered by -O1
but couldn't make it go away. Does anybody know why this is?

Regtests OK with FC41/x86_64 - OK for mainline?

Regards

Paul

Attachment: Change.Logs
Description: Binary data

diff --git a/gcc/fortran/check.cc b/gcc/fortran/check.cc
index 35458643835..0a5648a9d8f 100644
--- a/gcc/fortran/check.cc
+++ b/gcc/fortran/check.cc
@@ -2442,6 +2442,16 @@ gfc_check_co_broadcast (gfc_expr *a, gfc_expr *source_image, gfc_expr *stat,
 }
 
 
+/* 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;
+};
+
+
 bool
 gfc_check_co_reduce (gfc_expr *a, gfc_expr *op, gfc_expr *result_image,
 		     gfc_expr *stat, gfc_expr *errmsg)
@@ -2567,24 +2577,12 @@ 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)
@@ -5135,6 +5133,207 @@ 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)
+{
+  symbol_attribute attr;
+  gfc_formal_arglist *formal;
+  gfc_symbol *sym;
+
+  if (array->ts.type == BT_CLASS)
+    {
+      gfc_error ("The ARRAY argument at %L of REDUCE shall not be polymorphic",
+		 &array->where);
+      return false;
+    }
+
+  if (!gfc_resolve_expr (operation))
+    return false;
+
+  attr = gfc_expr_attr (operation);
+  if (!attr.pure || !attr.function)
+    {
+      gfc_error ("OPERATION argument at %L must be a PURE function",
+		 &operation->where);
+      return false;
+    }
+
+  if (attr.dimension)
+    {
+      gfc_error ("OPERATION argument at %L must be a SCALAR function",
+		 &operation->where);
+      return false;
+    }
+
+  if (attr.intrinsic)
+    {
+      /* 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 REDUCE",
+		 operation->symtree->n.sym->name, &operation->where);
+      return false;
+    }
+
+  if (gfc_is_proc_ptr_comp (operation))
+    {
+      gfc_component *comp = gfc_get_proc_ptr_comp (operation);
+      sym = comp->ts.interface;
+    }
+  else
+    sym = operation->symtree->n.sym;
+
+  formal = sym->formal;
+
+  if (!formal || !formal->next || formal->next->next)
+    {
+      gfc_error ("The function passed as OPERATION at %L shall have two "
+		 "arguments", &operation->where);
+      return false;
+    }
+
+  if (sym->result->ts.type == BT_UNKNOWN)
+    gfc_set_default_type (sym->result, 0, NULL);
+
+  if (!gfc_compare_types (&array->ts, &sym->result->ts))
+    {
+      gfc_error ("The ARRAY argument at %L has type %s but the function passed as "
+		 "OPERATION at %L returns %s",
+		 &array->where, gfc_typename (array), &operation->where,
+		 gfc_typename (&sym->result->ts));
+      return false;
+    }
+
+  if (!gfc_compare_types (&array->ts, &formal->sym->ts)
+      || !gfc_compare_types (&array->ts, &formal->next->sym->ts))
+    {
+      gfc_error ("The function passed as OPERATION at %L has arguments of type "
+		 "%s and %s but shall have type %s", &operation->where,
+		 gfc_typename (&formal->sym->ts),
+		 gfc_typename (&formal->next->sym->ts), gfc_typename (array));
+      return false;
+    }
+
+  if (operation->rank || attr.allocatable || attr.pointer)
+    {
+      gfc_error ("The function passed as OPERATION at %L shall be scalar, "
+		 "non-allocatable and non-pointer", &operation->where);
+      return false;
+    }
+
+  if (formal->sym->attr.dimension || formal->next->sym->attr.dimension)
+    {
+      gfc_error ("The function passed as OPERATION at %L shall have scalar "
+		 "arguments and return a scalar", &operation->where);
+      return false;
+    }
+
+  if (formal->sym->attr.allocatable || formal->sym->attr.allocatable
+      || formal->sym->attr.pointer || formal->sym->attr.pointer
+      || formal->sym->attr.optional || formal->sym->attr.optional
+      || formal->sym->ts.type == BT_CLASS || formal->sym->ts.type == BT_CLASS)
+    {
+      gfc_error ("Each argument of OPERATION at %L shall be a scalar, "
+		 "non-allocatable, non-pointer, non-polymorphic and "
+		 "nonoptional", &operation->where);
+      return false;
+    }
+
+  if (formal->sym->attr.value != formal->next->sym->attr.value)
+    {
+      gfc_error ("The function passed as OPERATION at %L shall have the VALUE "
+		 "attribute either for none or both arguments",
+		 &operation->where);
+      return false;
+    }
+
+  if (formal->sym->attr.target != formal->next->sym->attr.target)
+    {
+      gfc_error ("The function passed as OPERATION at %L shall have the TARGET "
+		 "attribute either for none or both arguments",
+		 &operation->where);
+      return false;
+    }
+
+  if (formal->sym->attr.asynchronous != formal->next->sym->attr.asynchronous)
+    {
+      gfc_error ("The function passed as OPERATION at %L shall have the "
+		 "ASYNCHRONOUS attribute either for none or both arguments",
+		 &operation->where);
+      return false;
+    }
+
+  if (array->ts.type == BT_CHARACTER)
+    {
+      unsigned long actual_size, formal_size1, formal_size2, result_size;
+
+      actual_size = get_ul_from_cst_cl (array->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 ARRAY argument at %L and of the "
+		     "arguments of the OPERATION at %L shall be the same",
+		     &array->where, &operation->where);
+	  return false;
+	}
+
+      if (actual_size && result_size && actual_size != result_size)
+	{
+	  gfc_error ("The character length of the ARRAY argument at %L and of the "
+		     "function result of the OPERATION at %L shall be the same",
+		     &array->where, &operation->where);
+	  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)
 {
diff --git a/gcc/fortran/gfortran.h b/gcc/fortran/gfortran.h
index 7c6e9b637db..5ef70378b1b 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 30f532b5766..678dc78b346 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,17 @@ 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 34a0248adbd..fec1c24a099 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 55f7e1911ec..3ba614dc9ad 100644
--- a/gcc/fortran/iresolve.cc
+++ b/gcc/fortran/iresolve.cc
@@ -2408,6 +2408,219 @@ 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;
+
+  /* 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 ();
+  if_block->block->expr1->expr_type = EXPR_FUNCTION;
+  if_block->block->expr1->where = gfc_current_locus;
+  gfc_get_sym_tree ("present", sub_ns, &if_block->block->expr1->symtree, false);
+  if_block->block->expr1->symtree->n.sym->attr.flavor = FL_PROCEDURE;
+  if_block->block->expr1->symtree->n.sym->attr.intrinsic = 1;
+  if_block->block->expr1->ts.type = BT_LOGICAL;
+  if_block->block->expr1->ts.kind = gfc_default_logical_kind;
+  if_block->block->expr1->value.function.isym = gfc_intrinsic_function_by_id (GFC_ISYM_PRESENT);
+  if_block->block->expr1->value.function.actual = gfc_get_actual_arglist ();
+  if_block->block->expr1->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 ();
+  if_block->block->next->expr2->expr_type = EXPR_FUNCTION;
+  if_block->block->next->expr2->where = gfc_current_locus;
+  if_block->block->next->expr2->ts = operation->ts;
+  gfc_get_sym_tree (operation->name, ns, &if_block->block->next->expr2->symtree, false);
+  if_block->block->next->expr2->value.function.esym = if_block->block->next->expr2->symtree->n.sym;
+  if_block->block->next->expr2->value.function.actual = gfc_get_actual_arglist ();
+  if_block->block->next->expr2->value.function.actual->expr
+					= gfc_lval_expr_from_sym (a);
+  if_block->block->next->expr2->value.function.actual->next = gfc_get_actual_arglist ();
+  if_block->block->next->expr2->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 = gfc_default_integer_kind;
+      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 d965539f11e..6fca3891477 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
@@ -8401,6 +8407,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)
@@ -8585,6 +8592,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.  */
@@ -8769,10 +8787,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)
@@ -8797,6 +8817,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 373a0678a2e..99a46dd28c1 100644
--- a/gcc/fortran/trans-intrinsic.cc
+++ b/gcc/fortran/trans-intrinsic.cc
@@ -3669,6 +3669,7 @@ gfc_get_symbol_for_expr (gfc_expr * expr, bool ignore_optional)
   sym->attr.proc = PROC_INTRINSIC;
   sym->attr.flavor = FL_PROCEDURE;
   sym->result = sym;
+
   if (expr->rank > 0)
     {
       sym->attr.dimension = 1;
@@ -9296,6 +9297,8 @@ gfc_conv_intrinsic_sr_kind (gfc_se *se, gfc_expr *expr)
 }
 
 
+/* Generate code for REDUCE function with a scalar result. This is handled as a
+   special case within gfc_conv_procedure call.  */
 /* Generate code for TRIM (A) intrinsic function.  */
 
 static void
@@ -10806,6 +10809,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 +11482,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 +11939,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 00000000000..a2dbb646e55
--- /dev/null
+++ b/gcc/testsuite/gfortran.dg/reduce_1.f90
@@ -0,0 +1,203 @@
+! { dg-do run }
+!
+! Test results from the F2018 intrinsic REDUCE
+!
+! Contributed by Paul Thomas  <pa...@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, j
+   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 :: ivec(:)
+   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 00000000000..fb300da2e7f
--- /dev/null
+++ b/gcc/testsuite/gfortran.dg/reduce_2.f90
@@ -0,0 +1,142 @@
+! { dg-do compile }
+!
+! Test argument compliance for the F2018 intrinsic REDUCE
+!
+! Contributed by Paul Thomas  <pa...@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" }
+
+! OPERATION argument at (1) must be a SCALAR function
+  print *, reduce (i, vadd) ! { dg-error "must be a SCALAR function" }
+
+! 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 be scalar, non-allocatable and non-pointer
+  print *, reduce (i, add_a) ! { dg-error "shall be scalar, non-allocatable" }
+
+! The function passed as OPERATION at (1) shall have scalar arguments and return a scalar
+  print *, reduce (i, add_array) ! { dg-error "shall have scalar arguments" }
+
+! Each argument of OPERATION at (1) shall be a scalar, non-allocatable, non-pointer, 
+! non-polymorphic and nonoptional
+  print *, reduce (i, add_optional) ! { dg-error "shall be a scalar, non-allocatable, non-pointer" }
+
+! 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), pointer :: 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
diff --git a/libgfortran/Makefile.am b/libgfortran/Makefile.am
index 52bd6ea641f..21b35c76a06 100644
--- a/libgfortran/Makefile.am
+++ b/libgfortran/Makefile.am
@@ -149,6 +149,7 @@ intrinsics/spread_generic.c \
 intrinsics/string_intrinsics.c \
 intrinsics/rand.c \
 intrinsics/random.c \
+intrinsics/reduce.c \
 intrinsics/reshape_generic.c \
 intrinsics/reshape_packed.c \
 intrinsics/selected_int_kind.f90 \
diff --git a/libgfortran/Makefile.in b/libgfortran/Makefile.in
index 9b9ed2f86ba..6a63d8876b1 100644
--- a/libgfortran/Makefile.in
+++ b/libgfortran/Makefile.in
@@ -612,7 +612,7 @@ am__objects_59 = intrinsics/associated.lo intrinsics/abort.lo \
 	intrinsics/move_alloc.lo intrinsics/pack_generic.lo \
 	intrinsics/selected_char_kind.lo intrinsics/size.lo \
 	intrinsics/spread_generic.lo intrinsics/string_intrinsics.lo \
-	intrinsics/rand.lo intrinsics/random.lo \
+	intrinsics/rand.lo intrinsics/random.lo intrinsics/reduce.lo \
 	intrinsics/reshape_generic.lo intrinsics/reshape_packed.lo \
 	intrinsics/selected_int_kind.lo \
 	intrinsics/selected_real_kind.lo intrinsics/trigd.lo \
@@ -1028,7 +1028,7 @@ gfor_helper_src = intrinsics/associated.c intrinsics/abort.c \
 	intrinsics/move_alloc.c intrinsics/pack_generic.c \
 	intrinsics/selected_char_kind.c intrinsics/size.c \
 	intrinsics/spread_generic.c intrinsics/string_intrinsics.c \
-	intrinsics/rand.c intrinsics/random.c \
+	intrinsics/rand.c intrinsics/random.c intrinsics/reduce.c \
 	intrinsics/reshape_generic.c intrinsics/reshape_packed.c \
 	intrinsics/selected_int_kind.f90 \
 	intrinsics/selected_real_kind.f90 intrinsics/trigd.c \
@@ -3476,6 +3476,8 @@ intrinsics/rand.lo: intrinsics/$(am__dirstamp) \
 	intrinsics/$(DEPDIR)/$(am__dirstamp)
 intrinsics/random.lo: intrinsics/$(am__dirstamp) \
 	intrinsics/$(DEPDIR)/$(am__dirstamp)
+intrinsics/reduce.lo: intrinsics/$(am__dirstamp) \
+	intrinsics/$(DEPDIR)/$(am__dirstamp)
 intrinsics/reshape_generic.lo: intrinsics/$(am__dirstamp) \
 	intrinsics/$(DEPDIR)/$(am__dirstamp)
 intrinsics/reshape_packed.lo: intrinsics/$(am__dirstamp) \
@@ -4591,6 +4593,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@intrinsics/$(DEPDIR)/perror.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@intrinsics/$(DEPDIR)/rand.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@intrinsics/$(DEPDIR)/random.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@intrinsics/$(DEPDIR)/reduce.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@intrinsics/$(DEPDIR)/rename.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@intrinsics/$(DEPDIR)/reshape_generic.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@intrinsics/$(DEPDIR)/reshape_packed.Plo@am__quote@
diff --git a/libgfortran/gfortran.map b/libgfortran/gfortran.map
index 159a862e90a..7725e125882 100644
--- a/libgfortran/gfortran.map
+++ b/libgfortran/gfortran.map
@@ -2023,4 +2023,8 @@ GFORTRAN_15 {
     _gfortran_pow_m16_m4;
     _gfortran_pow_m16_m8;
     _gfortran_pow_m16_m16;
+    _gfortran_reduce;
+    _gfortran_reduce_scalar;
+    _gfortran_reduce_c;
+    _gfortran_reduce_scalar_c;
 } GFORTRAN_14;
diff --git a/libgfortran/intrinsics/reduce.c b/libgfortran/intrinsics/reduce.c
new file mode 100644
index 00000000000..ae80a22d992
--- /dev/null
+++ b/libgfortran/intrinsics/reduce.c
@@ -0,0 +1,292 @@
+/* Generic implementation of the reduce intrinsic
+   Copyright (C) 2002-2025 Free Software Foundation, Inc.
+   Contributed by Paul Thomas <pa...@gcc.gnu.org>
+
+This file is part of the GNU Fortran runtime library (libgfortran).
+
+Libgfortran is free software; you can redistribute it and/or
+modify it under the terms of the GNU General Public
+License as published by the Free Software Foundation; either
+version 3 of the License, or (at your option) any later version.
+
+Ligbfortran is distributed in the hope that it will be useful,
+but WITHOUT ANY WarrayANTY; without even the implied warrayanty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"
+#include <string.h>
+#include <stdio.h>
+
+typedef GFC_FULL_ARRAY_DESCRIPTOR(GFC_MAX_DIMENSIONS, char) parray;
+
+extern void reduce (parray *, parray *, void (*operation) (void *, void *, void *),
+		    int *, gfc_array_l4 *, void *, void *);
+export_proto(reduce);
+
+void
+reduce (parray *ret,
+	parray *array,
+	void (*operation) (void *, void *, void *),
+	GFC_INTEGER_4 *dim,
+	gfc_array_l4 *mask,
+	void *identity,
+	void *ordered __attribute__ ((unused)))
+{
+  GFC_LOGICAL_4 maskR = 0;
+  void *array_ptr;
+  void *buffer_ptr;
+  void *zero;
+  void *buffer;
+  void *res;
+  index_type ext0, ext1, ext2;
+  index_type str0, str1, str2;
+  index_type idx0, idx1, idx2;
+  index_type dimen, dimen_m1, ldx;
+  bool started;
+  bool masked = false;
+  bool dim_present = dim != NULL;
+  bool mask_present = mask != NULL;
+  bool identity_present = identity != NULL;
+  bool scalar_result;
+  int i;
+  int array_rank = (int)GFC_DESCRIPTOR_RANK (array);
+  size_t elem_len = GFC_DESCRIPTOR_SIZE (array);
+
+/* First part of the algorithm: Set up the indexing by reshaping array as
+   [product (extent[1:dim-1]), extent[dim], product (extent[dim+1:rank])].
+   If the products are empty (eg. dim = rank or scalar result), then they
+   default to unity.
+   The result, an intermediate result buffer and 'zero' are then allocated.
+
+   Second part of the algorithm: Do the reduction itself.
+     1] cycle over the first and third dimension
+     2] do the iteration along the second dimension, starting at 2nd element
+     3] iteration is controlled by the mask, if present.
+	(i) until the first true is encountered, 'buffer_ptr' is updated to
+	   point to the previous array element
+	(ii) at the second true the operation is performed
+	(iii) buffer_ptr is rest to point to the result buffer
+     4] otherwise iteration proceeds as if mask is all trues
+     5] effectively, at each call to the operation
+	*buffer_ptr = operation (*buffer_ptr, array(idx1, idx2, idx3))
+     6] the operation wrapper renders this as a pointer operation by putting
+	the assignment within a subroutine, thereby making the iteration type
+	agnostic.
+     7] the operation wrapper has a second mode, when the second argument is
+	NULL, which effects *c = *a.
+     8] in both modes, allocatable components are corectly taken care of, which
+	is why the second mode is used to fill the result elements.  */
+
+  /* First part  */
+  if (dim_present)
+    {
+      if ((*dim < 1) || (*dim > array_rank))
+	runtime_error ("DIM in REDUCE intrinsic is less than 0 or greater than "
+		       "the rank of ARRAY");
+      dimen = (index_type) *dim;
+    }
+  else
+    dimen = 1;
+  dimen_m1 = dimen -1;
+
+  /* Set up the shape and strides for the reduction. This is made relatively
+     painless by the use of pointer arithmetic throughout (except for MASK,
+     whose type is known.  */
+  ext0 = ext1 = ext2 = 1;
+  str0 = str1 = str2 = 1;
+
+  scalar_result = (!dim_present && array_rank > 1) || array_rank == 1;
+
+  for (i = 0; i < array_rank; i++)
+    {
+      /* The shape of the reshaped array  */
+      index_type ext = GFC_DESCRIPTOR_EXTENT (array,i);
+      index_type str = GFC_DESCRIPTOR_STRIDE (array,i);
+
+      if (masked && (ext != GFC_DESCRIPTOR_EXTENT (mask, i)))
+	runtime_error ("shape mismatch between ARRAY and MASK in REDUCE "
+		       "intrinsic");
+
+      if (scalar_result)
+	{
+	  ext1 *= ext;
+	  continue;
+	}
+      else if (i < dimen_m1)
+	ext0 *= ext;
+      else if (i == dimen_m1)
+	ext1 = ext;
+      else
+	ext2 *= ext;
+
+      /* The dimensions of the return array.  */
+      if (i < (int)(dimen - 1))
+	GFC_DIMENSION_SET(ret->dim[i], 0, ext - 1, str);
+      else if (i < array_rank - 1)
+	GFC_DIMENSION_SET (ret->dim[i], 0, ext - 1, str);
+    }
+
+  if (!scalar_result)
+    {
+      str1 = GFC_DESCRIPTOR_STRIDE (array, dimen_m1);
+      if (dimen < array_rank)
+	str2 = GFC_DESCRIPTOR_STRIDE (array, dimen);
+      else
+	str2 = 1;
+    }
+
+  /* Allocate the result data, the result buffer and zero.  */
+  if (ret->base_addr == NULL)
+    ret->base_addr = calloc ((size_t)(ext0 * ext2), elem_len);
+  buffer = calloc (1, elem_len);
+  zero = calloc (1, elem_len);
+
+  /* Second part  */
+  for (idx0 = 0; idx0 < ext0; idx0++)
+    {
+      for (idx2 = 0; idx2 < ext2; idx2++)
+	{
+	  ldx = idx0 + idx2 * str2;
+	  if (mask_present)
+	    maskR = mask->base_addr[ldx];
+
+	  started = (mask_present && maskR) || !mask_present;
+
+	  buffer_ptr = array->base_addr
+			+ (size_t)((idx0 + idx2 * str2) * elem_len);
+
+	  /* Start the iteration.  */
+	  for (idx1 = 1; idx1 < ext1; idx1++)
+	    {
+	      /* If masked, cycle until after first element that is not masked
+		 out. Then set 'started' and cycle so that this becomes the first
+		 element in the reduction.  */
+	      ldx = idx0 + idx1 * str1 + idx2 * str2;
+	      if (mask_present)
+		maskR = mask->base_addr[ldx];
+
+	      array_ptr = array->base_addr
+			  + (size_t)((idx0 + idx1 * str1
+				      + idx2 * str2) * elem_len);
+	      if (!started)
+		{
+		  if (mask_present && maskR)
+		    started = true;
+		  buffer_ptr = array_ptr;
+		  continue;
+		}
+
+	      /* Call the operation, if not masked out, with the previous
+		 element or the buffer and current element as arguments. The
+		 result is stored in the buffer and the buffer_ptr set to
+		 point to buffer, instead of the previous array element.  */
+	      if ((mask_present && maskR) || !mask_present)
+		{
+		  operation (buffer_ptr, array_ptr, buffer);
+		  buffer_ptr = buffer;
+		}
+	    }
+
+	  /* If this result element is empty emit an error or, if available,
+	     set to identity.  */
+	  res = ret->base_addr + (size_t)((idx0 + idx2 * str1) * elem_len);
+	  if (started)
+	    {
+	      operation (buffer_ptr, NULL, res);
+	      operation (zero, NULL, buffer);
+	    }
+	  else
+	    {
+	      if (!identity_present)
+		runtime_error ("Empty column and no IDENTITY in REDUCE "
+			       "intrinsic");
+	      else
+		operation (identity, NULL, res);
+	    }
+	}
+    }
+  free (zero);
+  free (buffer);
+}
+
+
+extern void reduce_scalar (void *, parray *, void (*operation) (void *, void *, void *),
+			    int *, gfc_array_l4 *, void *, void *);
+export_proto(reduce_scalar);
+
+void
+reduce_scalar (void *res,
+	       parray *array,
+	       void (*operation) (void *, void *, void *),
+	       int *dim,
+	       gfc_array_l4 *mask,
+	       void *identity,
+	       void *ordered)
+{
+  parray ret;
+  ret.base_addr = NULL;
+  ret.dtype.rank = 0;
+  reduce (&ret, array, operation, dim, mask, identity, ordered);
+  memcpy (res, ret.base_addr, GFC_DESCRIPTOR_SIZE(array));
+  if (ret.base_addr) free (ret.base_addr);
+}
+
+extern void reduce_c (parray *, index_type, parray *,
+		      void (*operation) (void *, void *, void *),
+		      int *, gfc_array_l4 *, void *, void *,
+		      index_type, index_type);
+export_proto(reduce_c);
+
+void
+reduce_c (parray *ret,
+	  index_type ret_strlen,
+	  parray *array,
+	  void (*operation) (void *, void *, void *),
+	  GFC_INTEGER_4 *dim,
+	  gfc_array_l4 *mask,
+	  void *identity,
+	  void *ordered,
+	  index_type array_strlen __attribute__ ((unused)),
+	  index_type identity_strlen __attribute__ ((unused)))
+{
+  reduce (ret, array, operation, dim, mask, identity, ordered);
+}
+
+
+extern void reduce_scalar_c (void *, index_type, parray *,
+		      void (*operation) (void *, void *, void *),
+		      int *, gfc_array_l4 *, void *, void *,
+		      index_type, index_type);
+export_proto(reduce_scalar_c);
+
+
+void
+reduce_scalar_c (void *res,
+		 index_type res_strlen,
+		 parray *array,
+		 void (*operation) (void *, void *, void *),
+		 int *dim,
+		 gfc_array_l4 *mask,
+		 void *identity,
+		 void *ordered,
+		 index_type array_strlen __attribute__ ((unused)),
+		 index_type identity_strlen __attribute__ ((unused)))
+{
+  parray ret;
+  ret.base_addr = NULL;
+  ret.dtype.rank = 0;
+  reduce (&ret, array, operation, dim, mask, identity, ordered);
+  memcpy (res, ret.base_addr, GFC_DESCRIPTOR_SIZE(array));
+  if (ret.base_addr) free (ret.base_addr);
+}
+

Reply via email to