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

/* Implementation of the SUM intrinsic
   Copyright 2002, 2007 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 2 of the License, or (at your option) any later version.

In addition to the permissions in the GNU General Public License, the
Free Software Foundation gives you unlimited permission to link the
compiled version of this file into combinations with other programs,
and to distribute those combinations without any restriction coming
from the use of this file.  (The General Public License restrictions
do apply in other respects; for example, they cover modification of
the file, and distribution when not linked into a combine
executable.)

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.

You should have received a copy of the GNU General Public
License along with libgfortran; see the file COPYING.  If not,
write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
Boston, MA 02110-1301, USA.  */

#include "libgfortran.h"
#include <stdlib.h>
#include <assert.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;

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

  len = array->dim[dim].ubound + 1 - array->dim[dim].lbound;
  delta = array->dim[dim].stride;

  for (n = 0; n < dim; n++)
    {
      sstride[n] = array->dim[n].stride;
      extent[n] = array->dim[n].ubound + 1 - array->dim[n].lbound;

      if (extent[n] < 0)
	extent[n] = 0;
    }
  for (n = dim; n < rank; n++)
    {
      sstride[n] = array->dim[n + 1].stride;
      extent[n] =
        array->dim[n + 1].ubound + 1 - array->dim[n + 1].lbound;

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

  if (retarray->data == NULL)
    {
      size_t alloc_size;

      for (n = 0; n < rank; n++)
        {
          retarray->dim[n].lbound = 0;
          retarray->dim[n].ubound = extent[n]-1;
          if (n == 0)
            retarray->dim[n].stride = 1;
          else
            retarray->dim[n].stride = retarray->dim[n-1].stride * extent[n-1];
        }

      retarray->offset = 0;
      retarray->dtype = (array->dtype & ~GFC_DTYPE_RANK_MASK) | rank;

      alloc_size = sizeof (GFC_REAL_8) * retarray->dim[rank-1].stride
    		   * extent[rank-1];

      if (alloc_size == 0)
	{
	  /* Make sure we have a zero-sized array.  */
	  retarray->dim[0].lbound = 0;
	  retarray->dim[0].ubound = -1;
	  return;
	}
      else
	retarray->data = internal_malloc_size (alloc_size);
    }
  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 (compile_options.bounds_check)
	{
	  for (n=0; n < rank; n++)
	    {
	      index_type ret_extent;

	      ret_extent = retarray->dim[n].ubound + 1
		- retarray->dim[n].lbound;
	      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] = retarray->dim[n].stride;
      if (extent[n] <= 0)
        len = 0;
    }

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

  while (base)
    {
      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 look.  */
              base = NULL;
              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;
  int rank;
  int dim;
  index_type n;
  index_type len;
  index_type delta;
  index_type mdelta;
  int mask_kind;

  dim = (*pdim) - 1;
  rank = GFC_DESCRIPTOR_RANK (array) - 1;

  len = array->dim[dim].ubound + 1 - array->dim[dim].lbound;
  if (len <= 0)
    return;

  mbase = mask->data;

  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 = array->dim[dim].stride;
  mdelta = mask->dim[dim].stride * mask_kind;

  for (n = 0; n < dim; n++)
    {
      sstride[n] = array->dim[n].stride;
      mstride[n] = mask->dim[n].stride * mask_kind;
      extent[n] = array->dim[n].ubound + 1 - array->dim[n].lbound;

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

    }
  for (n = dim; n < rank; n++)
    {
      sstride[n] = array->dim[n + 1].stride;
      mstride[n] = mask->dim[n + 1].stride * mask_kind;
      extent[n] =
        array->dim[n + 1].ubound + 1 - array->dim[n + 1].lbound;

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

  if (retarray->data == NULL)
    {
      size_t alloc_size;

      for (n = 0; n < rank; n++)
        {
          retarray->dim[n].lbound = 0;
          retarray->dim[n].ubound = extent[n]-1;
          if (n == 0)
            retarray->dim[n].stride = 1;
          else
            retarray->dim[n].stride = retarray->dim[n-1].stride * extent[n-1];
        }

      alloc_size = sizeof (GFC_REAL_8) * retarray->dim[rank-1].stride
    		   * extent[rank-1];

      retarray->offset = 0;
      retarray->dtype = (array->dtype & ~GFC_DTYPE_RANK_MASK) | rank;

      if (alloc_size == 0)
	{
	  /* Make sure we have a zero-sized array.  */
	  retarray->dim[0].lbound = 0;
	  retarray->dim[0].ubound = -1;
	  return;
	}
      else
	retarray->data = internal_malloc_size (alloc_size);

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

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

	      ret_extent = retarray->dim[n].ubound + 1
		- retarray->dim[n].lbound;
	      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++)
            {
              index_type mask_extent, array_extent;

	      array_extent = array->dim[n].ubound + 1 - array->dim[n].lbound;
	      mask_extent = mask->dim[n].ubound + 1 - mask->dim[n].lbound;
	      if (array_extent != mask_extent)
		runtime_error ("Incorrect extent in MASK argument of"
			       " SUM intrinsic in dimension %ld:"
			       " is %ld, should be %ld", (long int) n + 1,
			       (long int) mask_extent, (long int) array_extent);
	    }
	}
    }

  for (n = 0; n < rank; n++)
    {
      count[n] = 0;
      dstride[n] = retarray->dim[n].stride;
      if (extent[n] <= 0)
        return;
    }

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

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

  result = 0;
        if (len <= 0)
	  *dest = 0;
	else
	  {
	    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 look.  */
              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 sstride[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)
    {
      sum_r8 (retarray, array, pdim);
      return;
    }
  /* Make dim zero based to avoid confusion.  */
  dim = (*pdim) - 1;
  rank = GFC_DESCRIPTOR_RANK (array) - 1;

  for (n = 0; n < dim; n++)
    {
      sstride[n] = array->dim[n].stride;
      extent[n] = array->dim[n].ubound + 1 - array->dim[n].lbound;

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

  for (n = dim; n < rank; n++)
    {
      sstride[n] = array->dim[n + 1].stride;
      extent[n] =
        array->dim[n + 1].ubound + 1 - array->dim[n + 1].lbound;

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

  if (retarray->data == NULL)
    {
      size_t alloc_size;

      for (n = 0; n < rank; n++)
        {
          retarray->dim[n].lbound = 0;
          retarray->dim[n].ubound = extent[n]-1;
          if (n == 0)
            retarray->dim[n].stride = 1;
          else
            retarray->dim[n].stride = retarray->dim[n-1].stride * extent[n-1];
        }

      retarray->offset = 0;
      retarray->dtype = (array->dtype & ~GFC_DTYPE_RANK_MASK) | rank;

      alloc_size = sizeof (GFC_REAL_8) * retarray->dim[rank-1].stride
    		   * extent[rank-1];

      if (alloc_size == 0)
	{
	  /* Make sure we have a zero-sized array.  */
	  retarray->dim[0].lbound = 0;
	  retarray->dim[0].ubound = -1;
	  return;
	}
      else
	retarray->data = internal_malloc_size (alloc_size);
    }
  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 (compile_options.bounds_check)
	{
	  for (n=0; n < rank; n++)
	    {
	      index_type ret_extent;

	      ret_extent = retarray->dim[n].ubound + 1
		- retarray->dim[n].lbound;
	      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] = retarray->dim[n].stride;
    }

  dest = retarray->data;

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