[gcc.git] / libgfortran / generated / sum_r8.c

/* Implementation of the SUM intrinsic
   Copyright (C) 2002-2018 Free Software Foundation, Inc.
   Contributed by Paul Brook <paul@nowt.org>

This file is part of the GNU Fortran 95 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.

Libgfortran is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty 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"


#if defined (HAVE_GFC_REAL_8) && defined (HAVE_GFC_REAL_8)


extern void sum_r8 (gfc_array_r8 * const restrict, 
	gfc_array_r8 * const restrict, const index_type * const restrict);
export_proto(sum_r8);

void
sum_r8 (gfc_array_r8 * const restrict retarray, 
	gfc_array_r8 * const restrict array, 
	const index_type * const restrict pdim)
{
  index_type count[GFC_MAX_DIMENSIONS];
  index_type extent[GFC_MAX_DIMENSIONS];
  index_type sstride[GFC_MAX_DIMENSIONS];
  index_type dstride[GFC_MAX_DIMENSIONS];
  const GFC_REAL_8 * restrict base;
  GFC_REAL_8 * restrict dest;
  index_type rank;
  index_type n;
  index_type len;
  index_type delta;
  index_type dim;
  int continue_loop;

#ifdef HAVE_BACK_ARG
  assert(back == 0);
#endif

  /* Make dim zero based to avoid confusion.  */
  rank = GFC_DESCRIPTOR_RANK (array) - 1;
  dim = (*pdim) - 1;

  if (unlikely (dim < 0 || dim > rank))
    {
      runtime_error ("Dim argument incorrect in SUM intrinsic: "
 		     "is %ld, should be between 1 and %ld",
		     (long int) dim + 1, (long int) rank + 1);
    }

  len = GFC_DESCRIPTOR_EXTENT(array,dim);
  if (len < 0)
    len = 0;
  delta = GFC_DESCRIPTOR_STRIDE(array,dim);

  for (n = 0; n < dim; n++)
    {
      sstride[n] = GFC_DESCRIPTOR_STRIDE(array,n);
      extent[n] = GFC_DESCRIPTOR_EXTENT(array,n);

      if (extent[n] < 0)
	extent[n] = 0;
    }
  for (n = dim; n < rank; n++)
    {
      sstride[n] = GFC_DESCRIPTOR_STRIDE(array, n + 1);
      extent[n] = GFC_DESCRIPTOR_EXTENT(array, n + 1);

      if (extent[n] < 0)
	extent[n] = 0;
    }

  if (retarray->base_addr == NULL)
    {
      size_t alloc_size, str;

      for (n = 0; n < rank; n++)
	{
	  if (n == 0)
	    str = 1;
	  else
	    str = GFC_DESCRIPTOR_STRIDE(retarray,n-1) * extent[n-1];

	  GFC_DIMENSION_SET(retarray->dim[n], 0, extent[n] - 1, str);

	}

      retarray->offset = 0;
      GFC_DTYPE_COPY_SETRANK(retarray,array,rank);

      alloc_size = GFC_DESCRIPTOR_STRIDE(retarray,rank-1) * extent[rank-1];

      retarray->base_addr = xmallocarray (alloc_size, sizeof (GFC_REAL_8));
      if (alloc_size == 0)
	{
	  /* Make sure we have a zero-sized array.  */
	  GFC_DIMENSION_SET(retarray->dim[0], 0, -1, 1);
	  return;

	}
    }
  else
    {
      if (rank != GFC_DESCRIPTOR_RANK (retarray))
	runtime_error ("rank of return array incorrect in"
		       " SUM intrinsic: is %ld, should be %ld",
		       (long int) (GFC_DESCRIPTOR_RANK (retarray)),
		       (long int) rank);

      if (unlikely (compile_options.bounds_check))
	bounds_ifunction_return ((array_t *) retarray, extent,
				 "return value", "SUM");
    }

  for (n = 0; n < rank; n++)
    {
      count[n] = 0;
      dstride[n] = GFC_DESCRIPTOR_STRIDE(retarray,n);
      if (extent[n] <= 0)
	return;
    }

  base = array->base_addr;
  dest = retarray->base_addr;

  continue_loop = 1;
  while (continue_loop)
    {
      const GFC_REAL_8 * restrict src;
      GFC_REAL_8 result;
      src = base;
      {

  result = 0;
	if (len <= 0)
	  *dest = 0;
	else
	  {
	    for (n = 0; n < len; n++, src += delta)
	      {

  result += *src;
	      }
	    
	    *dest = result;
	  }
      }
      /* Advance to the next element.  */
      count[0]++;
      base += sstride[0];
      dest += dstride[0];
      n = 0;
      while (count[n] == extent[n])
	{
	  /* When we get to the end of a dimension, reset it and increment
	     the next dimension.  */
	  count[n] = 0;
	  /* We could precalculate these products, but this is a less
	     frequently used path so probably not worth it.  */
	  base -= sstride[n] * extent[n];
	  dest -= dstride[n] * extent[n];
	  n++;
	  if (n >= rank)
	    {
	      /* Break out of the loop.  */
	      continue_loop = 0;
	      break;
	    }
	  else
	    {
	      count[n]++;
	      base += sstride[n];
	      dest += dstride[n];
	    }
	}
    }
}


extern void msum_r8 (gfc_array_r8 * const restrict, 
	gfc_array_r8 * const restrict, const index_type * const restrict,
	gfc_array_l1 * const restrict);
export_proto(msum_r8);

void
msum_r8 (gfc_array_r8 * const restrict retarray, 
	gfc_array_r8 * const restrict array, 
	const index_type * const restrict pdim, 
	gfc_array_l1 * const restrict mask)
{
  index_type count[GFC_MAX_DIMENSIONS];
  index_type extent[GFC_MAX_DIMENSIONS];
  index_type sstride[GFC_MAX_DIMENSIONS];
  index_type dstride[GFC_MAX_DIMENSIONS];
  index_type mstride[GFC_MAX_DIMENSIONS];
  GFC_REAL_8 * restrict dest;
  const GFC_REAL_8 * restrict base;
  const GFC_LOGICAL_1 * restrict mbase;
  index_type rank;
  index_type dim;
  index_type n;
  index_type len;
  index_type delta;
  index_type mdelta;
  int mask_kind;

#ifdef HAVE_BACK_ARG
  assert (back == 0);
#endif
  dim = (*pdim) - 1;
  rank = GFC_DESCRIPTOR_RANK (array) - 1;


  if (unlikely (dim < 0 || dim > rank))
    {
      runtime_error ("Dim argument incorrect in SUM intrinsic: "
 		     "is %ld, should be between 1 and %ld",
		     (long int) dim + 1, (long int) rank + 1);
    }

  len = GFC_DESCRIPTOR_EXTENT(array,dim);
  if (len <= 0)
    return;

  mbase = mask->base_addr;

  mask_kind = GFC_DESCRIPTOR_SIZE (mask);

  if (mask_kind == 1 || mask_kind == 2 || mask_kind == 4 || mask_kind == 8
#ifdef HAVE_GFC_LOGICAL_16
      || mask_kind == 16
#endif
      )
    mbase = GFOR_POINTER_TO_L1 (mbase, mask_kind);
  else
    runtime_error ("Funny sized logical array");

  delta = GFC_DESCRIPTOR_STRIDE(array,dim);
  mdelta = GFC_DESCRIPTOR_STRIDE_BYTES(mask,dim);

  for (n = 0; n < dim; n++)
    {
      sstride[n] = GFC_DESCRIPTOR_STRIDE(array,n);
      mstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(mask,n);
      extent[n] = GFC_DESCRIPTOR_EXTENT(array,n);

      if (extent[n] < 0)
	extent[n] = 0;

    }
  for (n = dim; n < rank; n++)
    {
      sstride[n] = GFC_DESCRIPTOR_STRIDE(array,n + 1);
      mstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(mask, n + 1);
      extent[n] = GFC_DESCRIPTOR_EXTENT(array, n + 1);

      if (extent[n] < 0)
	extent[n] = 0;
    }

  if (retarray->base_addr == NULL)
    {
      size_t alloc_size, str;

      for (n = 0; n < rank; n++)
	{
	  if (n == 0)
	    str = 1;
	  else
	    str= GFC_DESCRIPTOR_STRIDE(retarray,n-1) * extent[n-1];

	  GFC_DIMENSION_SET(retarray->dim[n], 0, extent[n] - 1, str);

	}

      alloc_size = GFC_DESCRIPTOR_STRIDE(retarray,rank-1) * extent[rank-1];

      retarray->offset = 0;
      GFC_DTYPE_COPY_SETRANK(retarray,array,rank);

      if (alloc_size == 0)
	{
	  /* Make sure we have a zero-sized array.  */
	  GFC_DIMENSION_SET(retarray->dim[0], 0, -1, 1);
	  return;
	}
      else
	retarray->base_addr = xmallocarray (alloc_size, sizeof (GFC_REAL_8));

    }
  else
    {
      if (rank != GFC_DESCRIPTOR_RANK (retarray))
	runtime_error ("rank of return array incorrect in SUM intrinsic");

      if (unlikely (compile_options.bounds_check))
	{
	  bounds_ifunction_return ((array_t *) retarray, extent,
				   "return value", "SUM");
	  bounds_equal_extents ((array_t *) mask, (array_t *) array,
	  			"MASK argument", "SUM");
	}
    }

  for (n = 0; n < rank; n++)
    {
      count[n] = 0;
      dstride[n] = GFC_DESCRIPTOR_STRIDE(retarray,n);
      if (extent[n] <= 0)
	return;
    }

  dest = retarray->base_addr;
  base = array->base_addr;

  while (base)
    {
      const GFC_REAL_8 * restrict src;
      const GFC_LOGICAL_1 * restrict msrc;
      GFC_REAL_8 result;
      src = base;
      msrc = mbase;
      {

  result = 0;
	for (n = 0; n < len; n++, src += delta, msrc += mdelta)
	  {

  if (*msrc)
    result += *src;
	  }
	*dest = result;
      }
      /* Advance to the next element.  */
      count[0]++;
      base += sstride[0];
      mbase += mstride[0];
      dest += dstride[0];
      n = 0;
      while (count[n] == extent[n])
	{
	  /* When we get to the end of a dimension, reset it and increment
	     the next dimension.  */
	  count[n] = 0;
	  /* We could precalculate these products, but this is a less
	     frequently used path so probably not worth it.  */
	  base -= sstride[n] * extent[n];
	  mbase -= mstride[n] * extent[n];
	  dest -= dstride[n] * extent[n];
	  n++;
	  if (n >= rank)
	    {
	      /* Break out of the loop.  */
	      base = NULL;
	      break;
	    }
	  else
	    {
	      count[n]++;
	      base += sstride[n];
	      mbase += mstride[n];
	      dest += dstride[n];
	    }
	}
    }
}


extern void ssum_r8 (gfc_array_r8 * const restrict, 
	gfc_array_r8 * const restrict, const index_type * const restrict,
	GFC_LOGICAL_4 *);
export_proto(ssum_r8);

void
ssum_r8 (gfc_array_r8 * const restrict retarray, 
	gfc_array_r8 * const restrict array, 
	const index_type * const restrict pdim, 
	GFC_LOGICAL_4 * mask)
{
  index_type count[GFC_MAX_DIMENSIONS];
  index_type extent[GFC_MAX_DIMENSIONS];
  index_type dstride[GFC_MAX_DIMENSIONS];
  GFC_REAL_8 * restrict dest;
  index_type rank;
  index_type n;
  index_type dim;


  if (*mask)
    {
#ifdef HAVE_BACK_ARG
      sum_r8 (retarray, array, pdim, back);
#else
      sum_r8 (retarray, array, pdim);
#endif
      return;
    }
  /* Make dim zero based to avoid confusion.  */
  dim = (*pdim) - 1;
  rank = GFC_DESCRIPTOR_RANK (array) - 1;

  if (unlikely (dim < 0 || dim > rank))
    {
      runtime_error ("Dim argument incorrect in SUM intrinsic: "
 		     "is %ld, should be between 1 and %ld",
		     (long int) dim + 1, (long int) rank + 1);
    }

  for (n = 0; n < dim; n++)
    {
      extent[n] = GFC_DESCRIPTOR_EXTENT(array,n);

      if (extent[n] <= 0)
	extent[n] = 0;
    }

  for (n = dim; n < rank; n++)
    {
      extent[n] =
	GFC_DESCRIPTOR_EXTENT(array,n + 1);

      if (extent[n] <= 0)
	extent[n] = 0;
    }

  if (retarray->base_addr == NULL)
    {
      size_t alloc_size, str;

      for (n = 0; n < rank; n++)
	{
	  if (n == 0)
	    str = 1;
	  else
	    str = GFC_DESCRIPTOR_STRIDE(retarray,n-1) * extent[n-1];

	  GFC_DIMENSION_SET(retarray->dim[n], 0, extent[n] - 1, str);

	}

      retarray->offset = 0;
      GFC_DTYPE_COPY_SETRANK(retarray,array,rank);

      alloc_size = GFC_DESCRIPTOR_STRIDE(retarray,rank-1) * extent[rank-1];

      if (alloc_size == 0)
	{
	  /* Make sure we have a zero-sized array.  */
	  GFC_DIMENSION_SET(retarray->dim[0], 0, -1, 1);
	  return;
	}
      else
	retarray->base_addr = xmallocarray (alloc_size, sizeof (GFC_REAL_8));
    }
  else
    {
      if (rank != GFC_DESCRIPTOR_RANK (retarray))
	runtime_error ("rank of return array incorrect in"
		       " SUM intrinsic: is %ld, should be %ld",
		       (long int) (GFC_DESCRIPTOR_RANK (retarray)),
		       (long int) rank);

      if (unlikely (compile_options.bounds_check))
	{
	  for (n=0; n < rank; n++)
	    {
	      index_type ret_extent;

	      ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,n);
	      if (extent[n] != ret_extent)
		runtime_error ("Incorrect extent in return value of"
			       " SUM intrinsic in dimension %ld:"
			       " is %ld, should be %ld", (long int) n + 1,
			       (long int) ret_extent, (long int) extent[n]);
	    }
	}
    }

  for (n = 0; n < rank; n++)
    {
      count[n] = 0;
      dstride[n] = GFC_DESCRIPTOR_STRIDE(retarray,n);
    }

  dest = retarray->base_addr;

  while(1)
    {
      *dest = 0;
      count[0]++;
      dest += dstride[0];
      n = 0;
      while (count[n] == extent[n])
	{
	  /* When we get to the end of a dimension, reset it and increment
	     the next dimension.  */
	  count[n] = 0;
	  /* We could precalculate these products, but this is a less
	     frequently used path so probably not worth it.  */
	  dest -= dstride[n] * extent[n];
	  n++;
	  if (n >= rank)
	    return;
	  else
	    {
	      count[n]++;
	      dest += dstride[n];
	    }
      	}
    }
}

#endif
Commit	Line	Data
6de9cd9a	1	/* Implementation of the SUM intrinsic
85ec4feb	2	Copyright (C) 2002-2018 Free Software Foundation, Inc.
6de9cd9a DN	3	Contributed by Paul Brook <paul@nowt.org>
6de9cd9a DN	4
57dea9f6	5	This file is part of the GNU Fortran 95 runtime library (libgfortran).
6de9cd9a DN	6
6de9cd9a DN	7	Libgfortran is free software; you can redistribute it and/or
57dea9f6	8	modify it under the terms of the GNU General Public
6de9cd9a	9	License as published by the Free Software Foundation; either
748086b7	10	version 3 of the License, or (at your option) any later version.
6de9cd9a DN	11
	12	Libgfortran is distributed in the hope that it will be useful,
	13	but WITHOUT ANY WARRANTY; without even the implied warranty of
	14	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
57dea9f6	15	GNU General Public License for more details.
6de9cd9a	16
748086b7 JJ	17	Under Section 7 of GPL version 3, you are granted additional
	18	permissions described in the GCC Runtime Library Exception, version
	19	3.1, as published by the Free Software Foundation.
	20
	21	You should have received a copy of the GNU General Public License and
	22	a copy of the GCC Runtime Library Exception along with this program;
	23	see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
	24	<http://www.gnu.org/licenses/>. */
6de9cd9a	25
36ae8a61	26	#include "libgfortran.h"
6de9cd9a	27
7d7b8bfe	28
644cb69f FXC	29	#if defined (HAVE_GFC_REAL_8) && defined (HAVE_GFC_REAL_8)
	30
	31
64acfd99 JB	32	extern void sum_r8 (gfc_array_r8 * const restrict,
64acfd99 JB	33	gfc_array_r8 * const restrict, const index_type * const restrict);
7f68c75f	34	export_proto(sum_r8);
7d7b8bfe	35
6de9cd9a	36	void
64acfd99 JB	37	sum_r8 (gfc_array_r8 * const restrict retarray,
	38	gfc_array_r8 * const restrict array,
	39	const index_type * const restrict pdim)
6de9cd9a	40	{
e33e218b TK	41	index_type count[GFC_MAX_DIMENSIONS];
	42	index_type extent[GFC_MAX_DIMENSIONS];
	43	index_type sstride[GFC_MAX_DIMENSIONS];
	44	index_type dstride[GFC_MAX_DIMENSIONS];
64acfd99 JB	45	const GFC_REAL_8 * restrict base;
64acfd99 JB	46	GFC_REAL_8 * restrict dest;
6de9cd9a DN	47	index_type rank;
	48	index_type n;
	49	index_type len;
	50	index_type delta;
	51	index_type dim;
da96f5ab	52	int continue_loop;
6de9cd9a	53
64b1806b TK	54	#ifdef HAVE_BACK_ARG
	55	assert(back == 0);
	56	#endif
	57
6de9cd9a	58	/* Make dim zero based to avoid confusion. */
6de9cd9a	59	rank = GFC_DESCRIPTOR_RANK (array) - 1;
cfdf6ff6 TK	60	dim = (*pdim) - 1;
	61
	62	if (unlikely (dim < 0 \|\| dim > rank))
	63	{
	64	runtime_error ("Dim argument incorrect in SUM intrinsic: "
	65	"is %ld, should be between 1 and %ld",
	66	(long int) dim + 1, (long int) rank + 1);
	67	}
e33e218b	68
dfb55fdc	69	len = GFC_DESCRIPTOR_EXTENT(array,dim);
da96f5ab TK	70	if (len < 0)
da96f5ab TK	71	len = 0;
dfb55fdc	72	delta = GFC_DESCRIPTOR_STRIDE(array,dim);
6de9cd9a DN	73
	74	for (n = 0; n < dim; n++)
	75	{
dfb55fdc TK	76	sstride[n] = GFC_DESCRIPTOR_STRIDE(array,n);
dfb55fdc TK	77	extent[n] = GFC_DESCRIPTOR_EXTENT(array,n);
80ee04b9 TK	78
	79	if (extent[n] < 0)
	80	extent[n] = 0;
6de9cd9a DN	81	}
	82	for (n = dim; n < rank; n++)
	83	{
dfb55fdc TK	84	sstride[n] = GFC_DESCRIPTOR_STRIDE(array, n + 1);
dfb55fdc TK	85	extent[n] = GFC_DESCRIPTOR_EXTENT(array, n + 1);
80ee04b9 TK	86
	87	if (extent[n] < 0)
	88	extent[n] = 0;
6de9cd9a DN	89	}
6de9cd9a DN	90
21d1335b	91	if (retarray->base_addr == NULL)
6c167c45	92	{
dfb55fdc	93	size_t alloc_size, str;
80ee04b9	94
6c167c45	95	for (n = 0; n < rank; n++)
80927a56 JJ	96	{
80927a56 JJ	97	if (n == 0)
dfb55fdc	98	str = 1;
80927a56 JJ	99	else
80927a56 JJ	100	str = GFC_DESCRIPTOR_STRIDE(retarray,n-1) * extent[n-1];
dfb55fdc TK	101
	102	GFC_DIMENSION_SET(retarray->dim[n], 0, extent[n] - 1, str);
	103
80927a56	104	}
6c167c45	105
efd4dc1a	106	retarray->offset = 0;
fa3c4d47	107	GFC_DTYPE_COPY_SETRANK(retarray,array,rank);
80ee04b9	108
92e6f3a4	109	alloc_size = GFC_DESCRIPTOR_STRIDE(retarray,rank-1) * extent[rank-1];
80ee04b9	110
92e6f3a4	111	retarray->base_addr = xmallocarray (alloc_size, sizeof (GFC_REAL_8));
80ee04b9 TK	112	if (alloc_size == 0)
	113	{
	114	/* Make sure we have a zero-sized array. */
dfb55fdc	115	GFC_DIMENSION_SET(retarray->dim[0], 0, -1, 1);
80ee04b9	116	return;
dfb55fdc	117
80ee04b9	118	}
6c167c45	119	}
50dd63a9 TK	120	else
50dd63a9 TK	121	{
50dd63a9	122	if (rank != GFC_DESCRIPTOR_RANK (retarray))
fd6590f8	123	runtime_error ("rank of return array incorrect in"
ccacefc7 TK	124	" SUM intrinsic: is %ld, should be %ld",
	125	(long int) (GFC_DESCRIPTOR_RANK (retarray)),
	126	(long int) rank);
fd6590f8	127
9731c4a3	128	if (unlikely (compile_options.bounds_check))
16bff921 TK	129	bounds_ifunction_return ((array_t *) retarray, extent,
16bff921 TK	130	"return value", "SUM");
50dd63a9 TK	131	}
50dd63a9 TK	132
6de9cd9a DN	133	for (n = 0; n < rank; n++)
	134	{
	135	count[n] = 0;
dfb55fdc	136	dstride[n] = GFC_DESCRIPTOR_STRIDE(retarray,n);
6de9cd9a	137	if (extent[n] <= 0)
facc1285	138	return;
6de9cd9a DN	139	}
6de9cd9a DN	140
21d1335b TB	141	base = array->base_addr;
21d1335b TB	142	dest = retarray->base_addr;
6de9cd9a	143
da96f5ab TK	144	continue_loop = 1;
da96f5ab TK	145	while (continue_loop)
6de9cd9a	146	{
64acfd99	147	const GFC_REAL_8 * restrict src;
6de9cd9a DN	148	GFC_REAL_8 result;
	149	src = base;
	150	{
	151
	152	result = 0;
80927a56	153	if (len <= 0)
6de9cd9a DN	154	*dest = 0;
	155	else
	156	{
	157	for (n = 0; n < len; n++, src += delta)
	158	{
	159
	160	result += *src;
80927a56	161	}
0cd0559e	162
6de9cd9a DN	163	*dest = result;
	164	}
	165	}
	166	/* Advance to the next element. */
	167	count[0]++;
	168	base += sstride[0];
	169	dest += dstride[0];
	170	n = 0;
	171	while (count[n] == extent[n])
80927a56 JJ	172	{
	173	/* When we get to the end of a dimension, reset it and increment
	174	the next dimension. */
	175	count[n] = 0;
	176	/* We could precalculate these products, but this is a less
	177	frequently used path so probably not worth it. */
	178	base -= sstride[n] * extent[n];
	179	dest -= dstride[n] * extent[n];
	180	n++;
80dd631f	181	if (n >= rank)
80927a56	182	{
80dd631f	183	/* Break out of the loop. */
da96f5ab TK	184	continue_loop = 0;
da96f5ab TK	185	break;
80927a56 JJ	186	}
	187	else
	188	{
	189	count[n]++;
	190	base += sstride[n];
	191	dest += dstride[n];
	192	}
	193	}
6de9cd9a DN	194	}
	195	}
	196
7d7b8bfe	197
64acfd99 JB	198	extern void msum_r8 (gfc_array_r8 * const restrict,
64acfd99 JB	199	gfc_array_r8 * const restrict, const index_type * const restrict,
28dc6b33	200	gfc_array_l1 * const restrict);
7f68c75f	201	export_proto(msum_r8);
7d7b8bfe	202
6de9cd9a	203	void
64acfd99 JB	204	msum_r8 (gfc_array_r8 * const restrict retarray,
	205	gfc_array_r8 * const restrict array,
	206	const index_type * const restrict pdim,
28dc6b33	207	gfc_array_l1 * const restrict mask)
6de9cd9a	208	{
e33e218b TK	209	index_type count[GFC_MAX_DIMENSIONS];
	210	index_type extent[GFC_MAX_DIMENSIONS];
	211	index_type sstride[GFC_MAX_DIMENSIONS];
	212	index_type dstride[GFC_MAX_DIMENSIONS];
	213	index_type mstride[GFC_MAX_DIMENSIONS];
64acfd99 JB	214	GFC_REAL_8 * restrict dest;
64acfd99 JB	215	const GFC_REAL_8 * restrict base;
28dc6b33	216	const GFC_LOGICAL_1 * restrict mbase;
cfdf6ff6 TK	217	index_type rank;
cfdf6ff6 TK	218	index_type dim;
6de9cd9a DN	219	index_type n;
	220	index_type len;
	221	index_type delta;
	222	index_type mdelta;
28dc6b33	223	int mask_kind;
6de9cd9a	224
64b1806b TK	225	#ifdef HAVE_BACK_ARG
	226	assert (back == 0);
	227	#endif
6de9cd9a DN	228	dim = (*pdim) - 1;
6de9cd9a DN	229	rank = GFC_DESCRIPTOR_RANK (array) - 1;
e33e218b	230
cfdf6ff6 TK	231
	232	if (unlikely (dim < 0 \|\| dim > rank))
	233	{
	234	runtime_error ("Dim argument incorrect in SUM intrinsic: "
	235	"is %ld, should be between 1 and %ld",
	236	(long int) dim + 1, (long int) rank + 1);
	237	}
	238
dfb55fdc	239	len = GFC_DESCRIPTOR_EXTENT(array,dim);
6de9cd9a DN	240	if (len <= 0)
6de9cd9a DN	241	return;
28dc6b33	242
21d1335b	243	mbase = mask->base_addr;
28dc6b33 TK	244
	245	mask_kind = GFC_DESCRIPTOR_SIZE (mask);
	246
	247	if (mask_kind == 1 \|\| mask_kind == 2 \|\| mask_kind == 4 \|\| mask_kind == 8
	248	#ifdef HAVE_GFC_LOGICAL_16
	249	\|\| mask_kind == 16
	250	#endif
	251	)
	252	mbase = GFOR_POINTER_TO_L1 (mbase, mask_kind);
	253	else
	254	runtime_error ("Funny sized logical array");
	255
dfb55fdc TK	256	delta = GFC_DESCRIPTOR_STRIDE(array,dim);
dfb55fdc TK	257	mdelta = GFC_DESCRIPTOR_STRIDE_BYTES(mask,dim);
6de9cd9a DN	258
	259	for (n = 0; n < dim; n++)
	260	{
dfb55fdc TK	261	sstride[n] = GFC_DESCRIPTOR_STRIDE(array,n);
	262	mstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(mask,n);
	263	extent[n] = GFC_DESCRIPTOR_EXTENT(array,n);
80ee04b9 TK	264
	265	if (extent[n] < 0)
	266	extent[n] = 0;
	267
6de9cd9a DN	268	}
	269	for (n = dim; n < rank; n++)
	270	{
dfb55fdc TK	271	sstride[n] = GFC_DESCRIPTOR_STRIDE(array,n + 1);
	272	mstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(mask, n + 1);
	273	extent[n] = GFC_DESCRIPTOR_EXTENT(array, n + 1);
80ee04b9 TK	274
	275	if (extent[n] < 0)
	276	extent[n] = 0;
6de9cd9a DN	277	}
6de9cd9a DN	278
21d1335b	279	if (retarray->base_addr == NULL)
50dd63a9	280	{
dfb55fdc	281	size_t alloc_size, str;
80ee04b9	282
50dd63a9	283	for (n = 0; n < rank; n++)
80927a56 JJ	284	{
	285	if (n == 0)
	286	str = 1;
	287	else
	288	str= GFC_DESCRIPTOR_STRIDE(retarray,n-1) * extent[n-1];
dfb55fdc TK	289
	290	GFC_DIMENSION_SET(retarray->dim[n], 0, extent[n] - 1, str);
	291
80927a56	292	}
50dd63a9	293
92e6f3a4	294	alloc_size = GFC_DESCRIPTOR_STRIDE(retarray,rank-1) * extent[rank-1];
80ee04b9	295
efd4dc1a	296	retarray->offset = 0;
fa3c4d47	297	GFC_DTYPE_COPY_SETRANK(retarray,array,rank);
80ee04b9 TK	298
	299	if (alloc_size == 0)
	300	{
	301	/* Make sure we have a zero-sized array. */
dfb55fdc	302	GFC_DIMENSION_SET(retarray->dim[0], 0, -1, 1);
80ee04b9 TK	303	return;
	304	}
	305	else
92e6f3a4	306	retarray->base_addr = xmallocarray (alloc_size, sizeof (GFC_REAL_8));
80ee04b9	307
50dd63a9 TK	308	}
	309	else
	310	{
50dd63a9	311	if (rank != GFC_DESCRIPTOR_RANK (retarray))
fd6590f8 TK	312	runtime_error ("rank of return array incorrect in SUM intrinsic");
fd6590f8 TK	313
9731c4a3	314	if (unlikely (compile_options.bounds_check))
fd6590f8	315	{
16bff921 TK	316	bounds_ifunction_return ((array_t *) retarray, extent,
	317	"return value", "SUM");
	318	bounds_equal_extents ((array_t ) mask, (array_t ) array,
	319	"MASK argument", "SUM");
fd6590f8	320	}
50dd63a9 TK	321	}
50dd63a9 TK	322
6de9cd9a DN	323	for (n = 0; n < rank; n++)
	324	{
	325	count[n] = 0;
dfb55fdc	326	dstride[n] = GFC_DESCRIPTOR_STRIDE(retarray,n);
6de9cd9a	327	if (extent[n] <= 0)
80927a56	328	return;
6de9cd9a DN	329	}
6de9cd9a DN	330
21d1335b TB	331	dest = retarray->base_addr;
21d1335b TB	332	base = array->base_addr;
6de9cd9a DN	333
	334	while (base)
	335	{
64acfd99	336	const GFC_REAL_8 * restrict src;
28dc6b33	337	const GFC_LOGICAL_1 * restrict msrc;
6de9cd9a DN	338	GFC_REAL_8 result;
	339	src = base;
	340	msrc = mbase;
	341	{
	342
	343	result = 0;
036e1775	344	for (n = 0; n < len; n++, src += delta, msrc += mdelta)
6de9cd9a	345	{
6de9cd9a DN	346
	347	if (*msrc)
	348	result += *src;
6de9cd9a	349	}
036e1775	350	*dest = result;
6de9cd9a DN	351	}
	352	/* Advance to the next element. */
	353	count[0]++;
	354	base += sstride[0];
	355	mbase += mstride[0];
	356	dest += dstride[0];
	357	n = 0;
	358	while (count[n] == extent[n])
80927a56 JJ	359	{
	360	/* When we get to the end of a dimension, reset it and increment
	361	the next dimension. */
	362	count[n] = 0;
	363	/* We could precalculate these products, but this is a less
	364	frequently used path so probably not worth it. */
	365	base -= sstride[n] * extent[n];
	366	mbase -= mstride[n] * extent[n];
	367	dest -= dstride[n] * extent[n];
	368	n++;
80dd631f	369	if (n >= rank)
80927a56	370	{
80dd631f	371	/* Break out of the loop. */
80927a56 JJ	372	base = NULL;
	373	break;
	374	}
	375	else
	376	{
	377	count[n]++;
	378	base += sstride[n];
	379	mbase += mstride[n];
	380	dest += dstride[n];
	381	}
	382	}
6de9cd9a DN	383	}
6de9cd9a DN	384	}
644cb69f	385
97a62038 TK	386
	387	extern void ssum_r8 (gfc_array_r8 * const restrict,
	388	gfc_array_r8 * const restrict, const index_type * const restrict,
	389	GFC_LOGICAL_4 *);
	390	export_proto(ssum_r8);
	391
	392	void
	393	ssum_r8 (gfc_array_r8 * const restrict retarray,
	394	gfc_array_r8 * const restrict array,
	395	const index_type * const restrict pdim,
	396	GFC_LOGICAL_4 * mask)
	397	{
802367d7 TK	398	index_type count[GFC_MAX_DIMENSIONS];
802367d7 TK	399	index_type extent[GFC_MAX_DIMENSIONS];
802367d7 TK	400	index_type dstride[GFC_MAX_DIMENSIONS];
802367d7 TK	401	GFC_REAL_8 * restrict dest;
97a62038 TK	402	index_type rank;
97a62038 TK	403	index_type n;
802367d7 TK	404	index_type dim;
802367d7 TK	405
97a62038 TK	406
	407	if (*mask)
	408	{
64b1806b TK	409	#ifdef HAVE_BACK_ARG
	410	sum_r8 (retarray, array, pdim, back);
	411	#else
97a62038	412	sum_r8 (retarray, array, pdim);
64b1806b	413	#endif
97a62038 TK	414	return;
97a62038 TK	415	}
802367d7 TK	416	/* Make dim zero based to avoid confusion. */
	417	dim = (*pdim) - 1;
	418	rank = GFC_DESCRIPTOR_RANK (array) - 1;
	419
cfdf6ff6 TK	420	if (unlikely (dim < 0 \|\| dim > rank))
	421	{
	422	runtime_error ("Dim argument incorrect in SUM intrinsic: "
	423	"is %ld, should be between 1 and %ld",
	424	(long int) dim + 1, (long int) rank + 1);
	425	}
	426
802367d7 TK	427	for (n = 0; n < dim; n++)
802367d7 TK	428	{
dfb55fdc	429	extent[n] = GFC_DESCRIPTOR_EXTENT(array,n);
802367d7 TK	430
	431	if (extent[n] <= 0)
	432	extent[n] = 0;
	433	}
	434
	435	for (n = dim; n < rank; n++)
	436	{
802367d7	437	extent[n] =
80927a56	438	GFC_DESCRIPTOR_EXTENT(array,n + 1);
802367d7 TK	439
802367d7 TK	440	if (extent[n] <= 0)
80927a56	441	extent[n] = 0;
802367d7	442	}
97a62038	443
21d1335b	444	if (retarray->base_addr == NULL)
97a62038	445	{
dfb55fdc	446	size_t alloc_size, str;
802367d7 TK	447
802367d7 TK	448	for (n = 0; n < rank; n++)
80927a56 JJ	449	{
	450	if (n == 0)
	451	str = 1;
	452	else
	453	str = GFC_DESCRIPTOR_STRIDE(retarray,n-1) * extent[n-1];
dfb55fdc TK	454
	455	GFC_DIMENSION_SET(retarray->dim[n], 0, extent[n] - 1, str);
	456
80927a56	457	}
802367d7	458
97a62038	459	retarray->offset = 0;
fa3c4d47	460	GFC_DTYPE_COPY_SETRANK(retarray,array,rank);
802367d7	461
92e6f3a4	462	alloc_size = GFC_DESCRIPTOR_STRIDE(retarray,rank-1) * extent[rank-1];
802367d7 TK	463
	464	if (alloc_size == 0)
	465	{
	466	/* Make sure we have a zero-sized array. */
dfb55fdc	467	GFC_DIMENSION_SET(retarray->dim[0], 0, -1, 1);
802367d7 TK	468	return;
	469	}
	470	else
92e6f3a4	471	retarray->base_addr = xmallocarray (alloc_size, sizeof (GFC_REAL_8));
97a62038 TK	472	}
	473	else
	474	{
802367d7 TK	475	if (rank != GFC_DESCRIPTOR_RANK (retarray))
	476	runtime_error ("rank of return array incorrect in"
	477	" SUM intrinsic: is %ld, should be %ld",
	478	(long int) (GFC_DESCRIPTOR_RANK (retarray)),
	479	(long int) rank);
	480
9731c4a3	481	if (unlikely (compile_options.bounds_check))
fd6590f8	482	{
802367d7 TK	483	for (n=0; n < rank; n++)
	484	{
	485	index_type ret_extent;
97a62038	486
dfb55fdc	487	ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,n);
802367d7 TK	488	if (extent[n] != ret_extent)
	489	runtime_error ("Incorrect extent in return value of"
	490	" SUM intrinsic in dimension %ld:"
	491	" is %ld, should be %ld", (long int) n + 1,
	492	(long int) ret_extent, (long int) extent[n]);
	493	}
fd6590f8 TK	494	}
fd6590f8 TK	495	}
97a62038	496
802367d7 TK	497	for (n = 0; n < rank; n++)
	498	{
	499	count[n] = 0;
dfb55fdc	500	dstride[n] = GFC_DESCRIPTOR_STRIDE(retarray,n);
802367d7 TK	501	}
802367d7 TK	502
21d1335b	503	dest = retarray->base_addr;
802367d7 TK	504
	505	while(1)
	506	{
	507	*dest = 0;
	508	count[0]++;
	509	dest += dstride[0];
	510	n = 0;
	511	while (count[n] == extent[n])
80927a56	512	{
802367d7	513	/* When we get to the end of a dimension, reset it and increment
80927a56 JJ	514	the next dimension. */
	515	count[n] = 0;
	516	/* We could precalculate these products, but this is a less
	517	frequently used path so probably not worth it. */
	518	dest -= dstride[n] * extent[n];
	519	n++;
80dd631f	520	if (n >= rank)
802367d7	521	return;
80927a56 JJ	522	else
	523	{
	524	count[n]++;
	525	dest += dstride[n];
	526	}
802367d7 TK	527	}
802367d7 TK	528	}
97a62038 TK	529	}
97a62038 TK	530
644cb69f	531	#endif