Hi All,

This is very much an early version of the F2018 REDUCE intrinsic. I am
posting it now because I have totally forgotten how to include new
functions in libgfortran.so. -static-libfortran works fine and the results
are the same as the other brands.

At present, it produces several of link warnings.
test_reduce.f90:23:2: warning: type of ‘_gfortran_reduce_scalar’ does not
match original declaration [-Wlto-type-mismatch]
   23 |   pure function add(i,j) result(sum_ij)
      |  ^
test_reduce.f90:23:2: note: return value type mismatch
test_reduce.f90:23:2: note: type ‘struct s’ should match type ‘int’
test_reduce.f90:23:2: note: ‘_gfortran_reduce_scalar’ was previously
declared here
test_reduce.f90:23:2: note: code may be misoptimized unless
‘-fno-strict-aliasing’ is used
/usr/bin/ld: warning: /tmp/ccfEUYXA.ltrans0.ltrans.o: requires executable
stack (because the .note.GNU-stack section is executable)

The last one is unavoidable because of the use of the wrapper for
'operation' that allows type agnostic use of pointer arithmetic in the
library functions. I am working on the type mismatch, which occurs when
different wrapper types appear in the same namespace.

Clearly there is a fair amount to do: clear the commented out
sections/lines, testcases and documentation.

Paul
diff --git a/gcc/fortran/check.cc b/gcc/fortran/check.cc
index 35458643835..ad617e177fb 100644
--- a/gcc/fortran/check.cc
+++ b/gcc/fortran/check.cc
@@ -5135,6 +5135,227 @@ 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 (formal->sym->attr.optional || formal->next->sym->attr.optional)
+    {
+      gfc_error ("The function passed as OPERATION at %L shall not have the "
+		 "OPTIONAL attribute for either of the arguments",
+		 &operation->where);
+      return false;
+    }
+
+  if (array->ts.type == BT_CHARACTER)
+    {
+      gfc_charlen *cl;
+      unsigned long actual_size, formal_size1, formal_size2, result_size;
+
+      cl = array->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;
+
+      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 INDENTITY", &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 5fe12764615..7df97e012c7 100644
--- a/gcc/fortran/gfortran.h
+++ b/gcc/fortran/gfortran.h
@@ -645,6 +645,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,
@@ -2519,6 +2520,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 dc60d98d51b..c561797d790 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..5671cdc7781 100644
--- a/gcc/fortran/iresolve.cc
+++ b/gcc/fortran/iresolve.cc
@@ -2408,6 +2408,182 @@ 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_gsymbol *gsym = NULL;
+
+  /* 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;
+  wrapper->formal = gfc_get_formal_arglist ();
+  wrapper->formal->sym = a;
+  gfc_set_sym_referenced (a);
+
+  /* Set up formal argument for the argument 'b'.  */
+  gfc_get_symbol ("b", sub_ns, &b);
+  b->ts = operation->ts;
+  b->attr.flavor = FL_VARIABLE;
+  b->attr.dummy = 1;
+  b->attr.artificial = 1;
+  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;
+  wrapper->formal->next->next = gfc_get_formal_arglist ();
+  wrapper->formal->next->next->sym = c;
+  gfc_set_sym_referenced (c);
+
+  /* The only code is the assignment c = operation (a, b).  */
+  sub_ns->code = gfc_get_code (EXEC_ASSIGN);
+  sub_ns->code->expr1 = gfc_lval_expr_from_sym (c);
+  sub_ns->code->expr2 = gfc_get_expr ();
+  sub_ns->code->expr2->expr_type = EXPR_FUNCTION;
+  sub_ns->code->expr2->where = gfc_current_locus;
+  sub_ns->code->expr2->ts = operation->ts;
+  gfc_get_sym_tree (operation->name, ns, &sub_ns->code->expr2->symtree, false);
+  sub_ns->code->expr2->value.function.esym = sub_ns->code->expr2->symtree->n.sym;
+  sub_ns->code->expr2->value.function.actual = gfc_get_actual_arglist ();
+  sub_ns->code->expr2->value.function.actual->expr
+					= gfc_lval_expr_from_sym (a);
+  sub_ns->code->expr2->value.function.actual->next = gfc_get_actual_arglist ();
+  sub_ns->code->expr2->value.function.actual->next->expr
+					= gfc_lval_expr_from_sym (b);
+
+  /* It is unexpected to have some symbols added at resolution. Commit the
+     changes in order to keep a clean state.  */
+  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 f923aeb9460..cf2149eb62b 100644
--- a/gcc/fortran/trans-expr.cc
+++ b/gcc/fortran/trans-expr.cc
@@ -8389,6 +8389,12 @@ 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));
+
+  /* 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 scalar_reduce = expr && !expr->rank && expr->value.function.isym
+		       && expr->value.function.isym->id == GFC_ISYM_REDUCE;
   if (byref)
     {
       if (se->direct_byref)
@@ -8573,6 +8579,16 @@ gfc_conv_procedure_call (gfc_se * se, gfc_symbol * sym,
       else if (ts.type == BT_CHARACTER)
 	vec_safe_push (retargs, len);
     }
+  else if (scalar_reduce)
+    {
+      type = gfc_typenode_for_spec (&expr->ts);
+      tmp = gfc_create_var (type, "sr");
+      se->expr = tmp;
+      result = tmp;
+      tmp =  gfc_build_addr_expr (pvoid_type_node, tmp);
+      vec_safe_push (retargs, tmp);
+    }
+
   gfc_free_interface_mapping (&mapping);
 
   /* We need to glom RETARGS + ARGLIST + STRINGARGS + APPEND_ARGS.  */
@@ -8785,6 +8801,12 @@ gfc_conv_procedure_call (gfc_se * se, gfc_symbol * sym,
 	  gfc_add_expr_to_block (&se->pre, tmp);
 	}
     }
+  else if (scalar_reduce)
+    {
+      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 51237d0d3be..d3e6d626afa 100644
--- a/gcc/fortran/trans-intrinsic.cc
+++ b/gcc/fortran/trans-intrinsic.cc
@@ -4250,6 +4250,20 @@ 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 0
+  if (expr->value.function.isym
+      && expr->value.function.isym->id == GFC_ISYM_REDUCE)
+    {
+      if (expr->value.function.actual
+	  && expr->value.function.actual->next
+	  && expr->value.function.actual->next->next
+	  && expr->value.function.actual->next->next->expr == NULL)
+	expr->rank = 0;
+      else if (expr->value.function.actual
+	       && expr->value.function.actual->expr)
+	expr->rank = expr->value.function.actual->expr->rank - 1;
+    }
+#endif
   if (expr->rank > 0)
     {
       sym->attr.dimension = 1;
@@ -9941,6 +9955,16 @@ 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.  */
+static void
+gfc_conv_intrinsic_scalar_reduce (gfc_se *se, gfc_expr *expr)
+{
+  gcc_assert (expr->rank == 0);
+  gfc_conv_intrinsic_funcall (se, expr);
+}
+
+
 /* Generate code for TRIM (A) intrinsic function.  */
 
 static void
@@ -11451,6 +11475,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);
@@ -11482,6 +11507,10 @@ gfc_conv_intrinsic_function (gfc_se * se, gfc_expr * expr)
     case GFC_ISYM_NONE:
       gcc_unreachable ();
 
+    case GFC_ISYM_REDUCE:
+      gfc_conv_intrinsic_scalar_reduce (se, expr);
+      break;
+
     case GFC_ISYM_REPEAT:
       gfc_conv_intrinsic_repeat (se, expr);
       break;
@@ -12575,6 +12604,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/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..4f47f3e7110 100644
--- a/libgfortran/Makefile.in
+++ b/libgfortran/Makefile.in
@@ -613,6 +613,7 @@ am__objects_59 = intrinsics/associated.lo intrinsics/abort.lo \
 	intrinsics/selected_char_kind.lo intrinsics/size.lo \
 	intrinsics/spread_generic.lo intrinsics/string_intrinsics.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 \
@@ -1029,6 +1030,7 @@ gfor_helper_src = intrinsics/associated.c intrinsics/abort.c \
 	intrinsics/selected_char_kind.c intrinsics/size.c \
 	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 \
 	intrinsics/selected_real_kind.f90 intrinsics/trigd.c \
@@ -3476,6 +3478,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) \
diff --git a/libgfortran/gfortran.map b/libgfortran/gfortran.map
index 159a862e90a..4948e578bb5 100644
--- a/libgfortran/gfortran.map
+++ b/libgfortran/gfortran.map
@@ -2023,4 +2023,6 @@ GFORTRAN_15 {
     _gfortran_pow_m16_m4;
     _gfortran_pow_m16_m8;
     _gfortran_pow_m16_m16;
+    _gfortran_reduce;
+    _gfortran_reduce_scalar;
 } GFORTRAN_14;
diff --git a/libgfortran/intrinsics/reduce.c b/libgfortran/intrinsics/reduce.c
new file mode 100644
index 00000000000..43d945f9b49
--- /dev/null
+++ b/libgfortran/intrinsics/reduce.c
@@ -0,0 +1,266 @@
+/* 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;
+
+static index_type
+product (index_type arr[], index_type start, index_type end)
+{
+  index_type i;
+  index_type prod;
+  prod = 1;
+  /* Retained fortran indexing for start and end!  */
+  for (i = start - 1; i < end; i++)
+    prod = prod * arr[i];
+  return prod;
+}
+
+
+extern void reduce (parray *, parray *, void (*operator) (void *, void *, void *),
+		    int *, gfc_array_l4 *, void *, void *);
+export_proto(reduce);
+
+void
+reduce (parray *ret,
+	parray *array,
+	void (*operator) (void *, void *, void *),
+	GFC_INTEGER_4 *dim,
+	gfc_array_l4 *mask,
+	void *identity,
+	void *ordered __attribute__ ((unused)))
+{
+  GFC_LOGICAL_4 maskR = 0;
+  index_type array_shape[GFC_MAX_DIMENSIONS];
+  void *ptr;
+  void *accumulator;
+  void *res;
+  index_type idx[3], idx_stride[3];
+  index_type dimen, dimen_m1, idx0, idx1, idx2, ldx;
+  bool started;
+  bool masked = false;
+  bool dim_present = dim != NULL;
+  bool mask_present = mask != NULL;
+  bool identity_present = identity != NULL;
+  int i;
+  int array_rank = (int)GFC_DESCRIPTOR_RANK (array);
+  size_t elem_len = GFC_DESCRIPTOR_SIZE (array);
+
+/* First part of the algorithm  */
+  for (i = 0; i < array_rank; i++)
+    array_shape[i] = GFC_DESCRIPTOR_EXTENT (array,i);
+
+  if (dim_present)
+    dimen = (index_type) *dim;
+  else
+    dimen = 1;
+  dimen_m1 = dimen -1;
+
+  for (i = 0; i < array_rank; i++)
+    {
+      array_shape[i] = GFC_DESCRIPTOR_EXTENT (array,i);
+      if (masked && (array_shape[i] != GFC_DESCRIPTOR_EXTENT (mask, i)))
+	runtime_error ("shape mismatch between ARRAY and MASK in REDUCE "
+		       "intrinsic");
+    }
+
+  /* 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.  */
+  idx[0] = idx[1] = idx[2] = 1;
+  idx_stride[0] = idx_stride[1] = idx_stride[2] = 1;
+
+  if ((!dim_present && array_rank > 1) || array_rank == 1)
+    idx[1] = product (array_shape, 1, (index_type)array_rank);
+  else
+    {
+      idx[0] = product (array_shape, 1 , dimen_m1);
+      idx[1] = product (array_shape, dimen , dimen);
+      idx[2] = product (array_shape, (dimen + 1), (index_type)array_rank);
+
+      idx_stride[1] = GFC_DESCRIPTOR_STRIDE (array, dimen_m1);
+      if (dimen < array_rank)
+        idx_stride[2] = GFC_DESCRIPTOR_STRIDE (array, dimen);
+      else
+	idx_stride[2] = 1;
+
+      for (i = 0; i < (array_rank - 1); i++)
+	{
+	  if (i < (int)(dimen - 1))
+	    GFC_DIMENSION_SET(ret->dim[i], 0,
+			      GFC_DESCRIPTOR_EXTENT (array,i) - 1,
+			      GFC_DESCRIPTOR_STRIDE (array,i));
+	  else
+	    GFC_DIMENSION_SET (ret->dim[i], 0,
+			       GFC_DESCRIPTOR_EXTENT (array,i + 1) - 1,
+			       GFC_DESCRIPTOR_STRIDE (array,i));
+	}
+    }
+
+  /* Allocate the stash for the iterated dimension and the result data  */
+  accumulator = malloc ((size_t)(elem_len * idx[1]));
+  if (ret->base_addr == NULL)
+    ret->base_addr = malloc ((size_t)(elem_len * idx[0] * idx[2]));
+
+
+  /* Do the reduction  */
+  for (idx0 = 0; idx0 < idx[0]; idx0++)
+    {
+      for (idx2 = 0; idx2 < idx[2]; idx2++)
+        {
+	  ldx = idx0 + idx2 * idx_stride[2];
+	  if (mask_present)
+	    maskR = mask->base_addr[ldx];
+
+	  started = (mask_present && maskR) || !mask_present;
+
+	  /* Fill the stash along the iterated dimension.  */
+	  for (idx1 = 0; idx1 < idx[1]; idx1++)
+	    {
+	      ptr = array->base_addr + (size_t)((idx0 + idx1 * idx_stride[1]
+						 + idx2 * idx_stride[2]) * elem_len);
+	      memcpy (accumulator + (size_t)(idx1 * elem_len), ptr, elem_len);
+	    }
+
+	  /* Start the iteration.  */
+	  for (idx1 = 1; idx1 < idx[1]; 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 * idx_stride[1] + idx2 * idx_stride[2];
+	      if (mask_present)
+		maskR = mask->base_addr[ldx];
+
+	      ptr = accumulator + (size_t)(idx1 * elem_len);
+	      if (!started)
+		{
+		  if (mask_present && maskR)
+		    started = true;
+		  continue;
+		}
+
+	      /* Call the operation, if not masked out, with the previous and
+		 current elements as arguments. Place the result in the
+		 previous element.  */
+	      if ((mask_present && maskR) || !mask_present)
+		operator (ptr - elem_len, ptr, ptr - elem_len);
+
+	      /* Unconditionally carry over the previous element into this one
+		 so that the accumulated result is taken forward.  */
+	      memcpy (ptr, ptr - elem_len, elem_len);
+	    }
+
+	  ptr = accumulator + (size_t)((idx[1] - 1) * elem_len);
+	  res = ret->base_addr + (size_t)((idx0 + idx2 * idx_stride[1]) * elem_len);
+
+	  /* If this result element is empty emit an error or, if available,
+	     set to identity.  */
+	  if (started)
+	    memcpy (res, ptr, elem_len);
+	  else
+	    {
+	      if (!identity_present)
+		runtime_error ("Empty column and no IDENTITY in REDUCE intrinsic");
+	      else
+		memcpy (res, identity, elem_len);
+	    }
+	}
+    }
+  free (accumulator);
+}
+
+
+extern void reduce_scalar (void *, parray *, void (*operator) (void *, void *, void *),
+			    int *, gfc_array_l4 *, void *, void *);
+export_proto(reduce_scalar);
+
+void
+reduce_scalar (void *res,
+	       parray *array,
+	       void (*operator) (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, operator, dim, mask, identity, ordered);
+  memcpy (res, ret.base_addr, GFC_DESCRIPTOR_SIZE(array));
+}
+
+extern void reduce_c (parray *, index_type, parray *,
+		      void (*operator) (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 (*operator) (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, operator, dim, mask, identity, ordered);
+}
+
+
+extern void reduce_scalar_c (void *, index_type, parray *,
+		      void (*operator) (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 (*operator) (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, operator, dim, mask, identity, ordered);
+  memcpy (res, ret.base_addr, GFC_DESCRIPTOR_SIZE(array));
+}
+

Attachment: Change.Logs
Description: Binary data

program test_reduce ! intrinsic reduce works with nagfor and Intel Fortran, not yet gfortran
implicit none
integer, parameter :: n = 3, vec(n) = [2, 5, 10], &
         mat(n,2) = reshape([vec,2*vec],shape=[size(vec),2])
integer, dimension(:,:), allocatable :: res2
logical, parameter :: t = .true., f = .false., keep(size(vec)) = [t,f,t]
logical, parameter :: keepM(n,2) = reshape([keep,keep],shape=[n,2])
integer :: i

print *, "my reduce ", reduce (vec, add, dim=1), "should be 17"

print *, "my reduce ", reduce (vec, mult, 1), "should be 100"

print *, "my reduce ", reduce (mat, add, 1), "should be 17 34"

print *, "my reduce ", reduce (mat, mult, 1), "should be 100 800"

print *, "my reduce ", reduce (mat, add, 2), "should be 6 15 30"

print *, "my reduce ", reduce (mat, mult, 2), "should be 8 50 200"

print *, "my reduce ", reduce (mat, add), "should be 51"

print *, "my reduce ", reduce (mat, mult), "should be 80000"

print *, "with mask"

print *, "my reduce ", reduce (vec,add, mask=keep), "should be 12"

print *, "my reduce ", reduce (vec,mult, mask=keep), "should be 20"

print *, "my reduce ", reduce (mat, add, mask=keepM), "should be 36"

print *, "my reduce ", reduce (mat, mult, mask=keepM), "should be 1600"

Print *, "higher dimension"

print *, "my reduce ", reduce (reshape ([(i, i=1,8)], [2,2,2]),mult), "should be 40320"

print *, "my reduce ", reduce (reshape ([(i, i=1,8)], [2,2,2]),mult,dim=2), "should be 3 8 35 48"

print *, "my reduce ", 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), "should be 3 -1 35 -1"

print *, "my reduce ", 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, ordered = f), "should be 9 12 14 12 18 20 22 24"
print *,  "my reduce ", 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, ordered = f), "should be 27 32 36 36"
!print*,mat(:,1) ! 2  5 10
!print*,mat(:,2) ! 4 10 20
! In lines below the two expressions give the same results.
! reduce applied to vec(:)
!print*,reduce(vec,add),sum(vec) ! 17
!print*,reduce(vec,mult),product(vec) ! 100
!print*,reduce(vec,add,mask=keep),sum(vec,mask=keep) ! 12
!print*,reduce(vec,mult,mask=keep),product(vec,mask=keep) ! 20
! reduce applied to mat(:,:)
!print*,reduce(mat,add),sum(mat) ! 51
!print*,reduce(mat,mult),product(mat) ! 80000
!print*,reduce(mat,add,dim=1),sum(mat,dim=1) ! 17 34
!print*,reduce(mat,mult,dim=1),product(mat,dim=1) ! 100 800
!print*,reduce(mat,add,dim=2),sum(mat,dim=2) ! 6 15 30
!print*,reduce(mat,mult,dim=2),product(mat,dim=2) ! 8 50 200

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
end

Reply via email to