[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 %d, should be %d",
		       GFC_DESCRIPTOR_RANK (retarray), 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 %d:"
			       " is %ld, should be %ld", 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 %d:"
			       " is %ld, should be %ld", 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 %d:"
			       " is %ld, should be %ld", 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 rank;
  index_type n;
  index_type dstride;
  GFC_REAL_8 *dest;

  if (*mask)
    {
      sum_r8 (retarray, array, pdim);
      return;
    }
    rank = GFC_DESCRIPTOR_RANK (array);
  if (rank <= 0)
    runtime_error ("Rank of array needs to be > 0");

  if (retarray->data == NULL)
    {
      retarray->dim[0].lbound = 0;
      retarray->dim[0].ubound = rank-1;
      retarray->dim[0].stride = 1;
      retarray->dtype = (retarray->dtype & ~GFC_DTYPE_RANK_MASK) | 1;
      retarray->offset = 0;
      retarray->data = internal_malloc_size (sizeof (GFC_REAL_8) * rank);
    }
  else
    {
      if (compile_options.bounds_check)
	{
	  int ret_rank;
	  index_type ret_extent;

	  ret_rank = GFC_DESCRIPTOR_RANK (retarray);
	  if (ret_rank != 1)
	    runtime_error ("rank of return array in SUM intrinsic"
			   " should be 1, is %d", ret_rank);

	  ret_extent = retarray->dim[0].ubound + 1 - retarray->dim[0].lbound;
	    if (ret_extent != rank)
	      runtime_error ("dimension of return array incorrect");
	}
    }
    dstride = retarray->dim[0].stride;
    dest = retarray->data;

    for (n = 0; n < rank; n++)
      dest[n * dstride] = 0 ;
}

#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 TK	118	runtime_error ("rank of return array incorrect in"
	119	" SUM intrinsic: is %d, should be %d",
	120	GFC_DESCRIPTOR_RANK (retarray), rank);
	121
	122	if (compile_options.bounds_check)
	123	{
	124	for (n=0; n < rank; n++)
	125	{
	126	index_type ret_extent;
	127
	128	ret_extent = retarray->dim[n].ubound + 1
	129	- retarray->dim[n].lbound;
	130	if (extent[n] != ret_extent)
	131	runtime_error ("Incorrect extent in return value of"
	132	" SUM intrinsic in dimension %d:"
	133	" is %ld, should be %ld", n + 1,
	134	(long int) ret_extent, (long int) extent[n]);
	135	}
	136	}
50dd63a9 TK	137	}
50dd63a9 TK	138
6de9cd9a DN	139	for (n = 0; n < rank; n++)
	140	{
	141	count[n] = 0;
	142	dstride[n] = retarray->dim[n].stride;
	143	if (extent[n] <= 0)
	144	len = 0;
	145	}
	146
	147	base = array->data;
	148	dest = retarray->data;
	149
	150	while (base)
	151	{
64acfd99	152	const GFC_REAL_8 * restrict src;
6de9cd9a DN	153	GFC_REAL_8 result;
	154	src = base;
	155	{
	156
	157	result = 0;
	158	if (len <= 0)
	159	*dest = 0;
	160	else
	161	{
	162	for (n = 0; n < len; n++, src += delta)
	163	{
	164
	165	result += *src;
	166	}
	167	*dest = result;
	168	}
	169	}
	170	/* Advance to the next element. */
	171	count[0]++;
	172	base += sstride[0];
	173	dest += dstride[0];
	174	n = 0;
	175	while (count[n] == extent[n])
	176	{
	177	/* When we get to the end of a dimension, reset it and increment
	178	the next dimension. */
	179	count[n] = 0;
	180	/* We could precalculate these products, but this is a less
5d7adf7a	181	frequently used path so probably not worth it. */
6de9cd9a DN	182	base -= sstride[n] * extent[n];
	183	dest -= dstride[n] * extent[n];
	184	n++;
	185	if (n == rank)
	186	{
	187	/* Break out of the look. */
	188	base = NULL;
	189	break;
	190	}
	191	else
	192	{
	193	count[n]++;
	194	base += sstride[n];
	195	dest += dstride[n];
	196	}
	197	}
	198	}
	199	}
	200
7d7b8bfe	201
64acfd99 JB	202	extern void msum_r8 (gfc_array_r8 * const restrict,
64acfd99 JB	203	gfc_array_r8 * const restrict, const index_type * const restrict,
28dc6b33	204	gfc_array_l1 * const restrict);
7f68c75f	205	export_proto(msum_r8);
7d7b8bfe	206
6de9cd9a	207	void
64acfd99 JB	208	msum_r8 (gfc_array_r8 * const restrict retarray,
	209	gfc_array_r8 * const restrict array,
	210	const index_type * const restrict pdim,
28dc6b33	211	gfc_array_l1 * const restrict mask)
6de9cd9a	212	{
e33e218b TK	213	index_type count[GFC_MAX_DIMENSIONS];
	214	index_type extent[GFC_MAX_DIMENSIONS];
	215	index_type sstride[GFC_MAX_DIMENSIONS];
	216	index_type dstride[GFC_MAX_DIMENSIONS];
	217	index_type mstride[GFC_MAX_DIMENSIONS];
64acfd99 JB	218	GFC_REAL_8 * restrict dest;
64acfd99 JB	219	const GFC_REAL_8 * restrict base;
28dc6b33	220	const GFC_LOGICAL_1 * restrict mbase;
6de9cd9a DN	221	int rank;
	222	int dim;
	223	index_type n;
	224	index_type len;
	225	index_type delta;
	226	index_type mdelta;
28dc6b33	227	int mask_kind;
6de9cd9a DN	228
	229	dim = (*pdim) - 1;
	230	rank = GFC_DESCRIPTOR_RANK (array) - 1;
e33e218b	231
6de9cd9a DN	232	len = array->dim[dim].ubound + 1 - array->dim[dim].lbound;
	233	if (len <= 0)
	234	return;
28dc6b33 TK	235
	236	mbase = mask->data;
	237
	238	mask_kind = GFC_DESCRIPTOR_SIZE (mask);
	239
	240	if (mask_kind == 1 \|\| mask_kind == 2 \|\| mask_kind == 4 \|\| mask_kind == 8
	241	#ifdef HAVE_GFC_LOGICAL_16
	242	\|\| mask_kind == 16
	243	#endif
	244	)
	245	mbase = GFOR_POINTER_TO_L1 (mbase, mask_kind);
	246	else
	247	runtime_error ("Funny sized logical array");
	248
6de9cd9a	249	delta = array->dim[dim].stride;
28dc6b33	250	mdelta = mask->dim[dim].stride * mask_kind;
6de9cd9a DN	251
	252	for (n = 0; n < dim; n++)
	253	{
	254	sstride[n] = array->dim[n].stride;
28dc6b33	255	mstride[n] = mask->dim[n].stride * mask_kind;
6de9cd9a	256	extent[n] = array->dim[n].ubound + 1 - array->dim[n].lbound;
80ee04b9 TK	257
	258	if (extent[n] < 0)
	259	extent[n] = 0;
	260
6de9cd9a DN	261	}
	262	for (n = dim; n < rank; n++)
	263	{
	264	sstride[n] = array->dim[n + 1].stride;
28dc6b33	265	mstride[n] = mask->dim[n + 1].stride * mask_kind;
6de9cd9a DN	266	extent[n] =
6de9cd9a DN	267	array->dim[n + 1].ubound + 1 - array->dim[n + 1].lbound;
80ee04b9 TK	268
	269	if (extent[n] < 0)
	270	extent[n] = 0;
6de9cd9a DN	271	}
6de9cd9a DN	272
50dd63a9 TK	273	if (retarray->data == NULL)
50dd63a9 TK	274	{
80ee04b9 TK	275	size_t alloc_size;
80ee04b9 TK	276
50dd63a9 TK	277	for (n = 0; n < rank; n++)
	278	{
	279	retarray->dim[n].lbound = 0;
	280	retarray->dim[n].ubound = extent[n]-1;
	281	if (n == 0)
	282	retarray->dim[n].stride = 1;
	283	else
	284	retarray->dim[n].stride = retarray->dim[n-1].stride * extent[n-1];
	285	}
	286
80ee04b9 TK	287	alloc_size = sizeof (GFC_REAL_8) * retarray->dim[rank-1].stride
	288	* extent[rank-1];
	289
efd4dc1a	290	retarray->offset = 0;
50dd63a9	291	retarray->dtype = (array->dtype & ~GFC_DTYPE_RANK_MASK) \| rank;
80ee04b9 TK	292
	293	if (alloc_size == 0)
	294	{
	295	/* Make sure we have a zero-sized array. */
	296	retarray->dim[0].lbound = 0;
	297	retarray->dim[0].ubound = -1;
	298	return;
	299	}
	300	else
	301	retarray->data = internal_malloc_size (alloc_size);
	302
50dd63a9 TK	303	}
	304	else
	305	{
50dd63a9	306	if (rank != GFC_DESCRIPTOR_RANK (retarray))
fd6590f8 TK	307	runtime_error ("rank of return array incorrect in SUM intrinsic");
	308
	309	if (compile_options.bounds_check)
	310	{
	311	for (n=0; n < rank; n++)
	312	{
	313	index_type ret_extent;
	314
	315	ret_extent = retarray->dim[n].ubound + 1
	316	- retarray->dim[n].lbound;
	317	if (extent[n] != ret_extent)
	318	runtime_error ("Incorrect extent in return value of"
	319	" SUM intrinsic in dimension %d:"
	320	" is %ld, should be %ld", n + 1,
	321	(long int) ret_extent, (long int) extent[n]);
	322	}
	323	for (n=0; n<= rank; n++)
	324	{
	325	index_type mask_extent, array_extent;
	326
	327	array_extent = array->dim[n].ubound + 1 - array->dim[n].lbound;
	328	mask_extent = mask->dim[n].ubound + 1 - mask->dim[n].lbound;
	329	if (array_extent != mask_extent)
	330	runtime_error ("Incorrect extent in MASK argument of"
	331	" SUM intrinsic in dimension %d:"
	332	" is %ld, should be %ld", n + 1,
	333	(long int) mask_extent, (long int) array_extent);
	334	}
	335	}
50dd63a9 TK	336	}
50dd63a9 TK	337
6de9cd9a DN	338	for (n = 0; n < rank; n++)
	339	{
	340	count[n] = 0;
	341	dstride[n] = retarray->dim[n].stride;
	342	if (extent[n] <= 0)
	343	return;
	344	}
	345
	346	dest = retarray->data;
	347	base = array->data;
6de9cd9a DN	348
	349	while (base)
	350	{
64acfd99	351	const GFC_REAL_8 * restrict src;
28dc6b33	352	const GFC_LOGICAL_1 * restrict msrc;
6de9cd9a DN	353	GFC_REAL_8 result;
	354	src = base;
	355	msrc = mbase;
	356	{
	357
	358	result = 0;
	359	if (len <= 0)
	360	*dest = 0;
	361	else
	362	{
	363	for (n = 0; n < len; n++, src += delta, msrc += mdelta)
	364	{
	365
	366	if (*msrc)
	367	result += *src;
	368	}
	369	*dest = result;
	370	}
	371	}
	372	/* Advance to the next element. */
	373	count[0]++;
	374	base += sstride[0];
	375	mbase += mstride[0];
	376	dest += dstride[0];
	377	n = 0;
	378	while (count[n] == extent[n])
	379	{
	380	/* When we get to the end of a dimension, reset it and increment
	381	the next dimension. */
	382	count[n] = 0;
	383	/* We could precalculate these products, but this is a less
5d7adf7a	384	frequently used path so probably not worth it. */
6de9cd9a DN	385	base -= sstride[n] * extent[n];
	386	mbase -= mstride[n] * extent[n];
	387	dest -= dstride[n] * extent[n];
	388	n++;
	389	if (n == rank)
	390	{
	391	/* Break out of the look. */
	392	base = NULL;
	393	break;
	394	}
	395	else
	396	{
	397	count[n]++;
	398	base += sstride[n];
	399	mbase += mstride[n];
	400	dest += dstride[n];
	401	}
	402	}
	403	}
	404	}
644cb69f	405
97a62038 TK	406
	407	extern void ssum_r8 (gfc_array_r8 * const restrict,
	408	gfc_array_r8 * const restrict, const index_type * const restrict,
	409	GFC_LOGICAL_4 *);
	410	export_proto(ssum_r8);
	411
	412	void
	413	ssum_r8 (gfc_array_r8 * const restrict retarray,
	414	gfc_array_r8 * const restrict array,
	415	const index_type * const restrict pdim,
	416	GFC_LOGICAL_4 * mask)
	417	{
	418	index_type rank;
	419	index_type n;
	420	index_type dstride;
	421	GFC_REAL_8 *dest;
	422
	423	if (*mask)
	424	{
	425	sum_r8 (retarray, array, pdim);
	426	return;
	427	}
	428	rank = GFC_DESCRIPTOR_RANK (array);
	429	if (rank <= 0)
	430	runtime_error ("Rank of array needs to be > 0");
	431
	432	if (retarray->data == NULL)
	433	{
	434	retarray->dim[0].lbound = 0;
	435	retarray->dim[0].ubound = rank-1;
	436	retarray->dim[0].stride = 1;
	437	retarray->dtype = (retarray->dtype & ~GFC_DTYPE_RANK_MASK) \| 1;
	438	retarray->offset = 0;
	439	retarray->data = internal_malloc_size (sizeof (GFC_REAL_8) * rank);
	440	}
	441	else
	442	{
fd6590f8 TK	443	if (compile_options.bounds_check)
	444	{
	445	int ret_rank;
	446	index_type ret_extent;
97a62038	447
fd6590f8 TK	448	ret_rank = GFC_DESCRIPTOR_RANK (retarray);
	449	if (ret_rank != 1)
	450	runtime_error ("rank of return array in SUM intrinsic"
	451	" should be 1, is %d", ret_rank);
97a62038	452
fd6590f8 TK	453	ret_extent = retarray->dim[0].ubound + 1 - retarray->dim[0].lbound;
	454	if (ret_extent != rank)
	455	runtime_error ("dimension of return array incorrect");
	456	}
	457	}
97a62038 TK	458	dstride = retarray->dim[0].stride;
	459	dest = retarray->data;
	460
	461	for (n = 0; n < rank; n++)
	462	dest[n * dstride] = 0 ;
	463	}
	464
644cb69f	465	#endif