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

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

This file is part of the GNU Fortran 95 runtime library (libgfor).

Libgfortran is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 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 Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public
License along with libgfor; see the file COPYING.LIB.  If not,
write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
Boston, MA 02111-1307, USA.  */

#include "config.h"
#include <stdlib.h>
#include <assert.h>
#include "libgfortran.h"


extern void __sum_r8 (gfc_array_r8 *, gfc_array_r8 *, index_type *);
export_proto_np(__sum_r8);

void
__sum_r8 (gfc_array_r8 *retarray, gfc_array_r8 *array, index_type *pdim)
{
  index_type count[GFC_MAX_DIMENSIONS - 1];
  index_type extent[GFC_MAX_DIMENSIONS - 1];
  index_type sstride[GFC_MAX_DIMENSIONS - 1];
  index_type dstride[GFC_MAX_DIMENSIONS - 1];
  GFC_REAL_8 *base;
  GFC_REAL_8 *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;
  assert (rank == GFC_DESCRIPTOR_RANK (retarray));
  if (array->dim[0].stride == 0)
    array->dim[0].stride = 1;
  if (retarray->dim[0].stride == 0)
    retarray->dim[0].stride = 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;
    }
  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 (retarray->data == NULL)
    {
      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->data
	 = internal_malloc_size (sizeof (GFC_REAL_8)
		 		 * retarray->dim[rank-1].stride
				 * extent[rank-1]);
      retarray->base = 0;
    }
          
  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)
    {
      GFC_REAL_8 *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 proabably 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 *, gfc_array_r8 *, index_type *,
						gfc_array_l4 *);
export_proto_np(__msum_r8);

void
__msum_r8 (gfc_array_r8 * retarray, gfc_array_r8 * array, index_type *pdim, gfc_array_l4 * mask)
{
  index_type count[GFC_MAX_DIMENSIONS - 1];
  index_type extent[GFC_MAX_DIMENSIONS - 1];
  index_type sstride[GFC_MAX_DIMENSIONS - 1];
  index_type dstride[GFC_MAX_DIMENSIONS - 1];
  index_type mstride[GFC_MAX_DIMENSIONS - 1];
  GFC_REAL_8 *dest;
  GFC_REAL_8 *base;
  GFC_LOGICAL_4 *mbase;
  int rank;
  int dim;
  index_type n;
  index_type len;
  index_type delta;
  index_type mdelta;

  dim = (*pdim) - 1;
  rank = GFC_DESCRIPTOR_RANK (array) - 1;
  assert (rank == GFC_DESCRIPTOR_RANK (retarray));
  if (array->dim[0].stride == 0)
    array->dim[0].stride = 1;
  if (retarray->dim[0].stride == 0)
    retarray->dim[0].stride = 1;

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

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

  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;
  mbase = mask->data;

  if (GFC_DESCRIPTOR_SIZE (mask) != 4)
    {
      /* This allows the same loop to be used for all logical types.  */
      assert (GFC_DESCRIPTOR_SIZE (mask) == 8);
      for (n = 0; n < rank; n++)
        mstride[n] <<= 1;
      mdelta <<= 1;
      mbase = (GFOR_POINTER_L8_TO_L4 (mbase));
    }

  while (base)
    {
      GFC_REAL_8 *src;
      GFC_LOGICAL_4 *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 proabably 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];
            }
        }
    }
}
Commit	Line	Data
6de9cd9a DN	1	/* Implementation of the SUM intrinsic
	2	Copyright 2002 Free Software Foundation, Inc.
	3	Contributed by Paul Brook <paul@nowt.org>
	4
	5	This file is part of the GNU Fortran 95 runtime library (libgfor).
	6
	7	Libgfortran is free software; you can redistribute it and/or
	8	modify it under the terms of the GNU Lesser General Public
	9	License as published by the Free Software Foundation; either
	10	version 2.1 of the License, or (at your option) any later version.
	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
	15	GNU Lesser General Public License for more details.
	16
	17	You should have received a copy of the GNU Lesser General Public
	18	License along with libgfor; see the file COPYING.LIB. If not,
	19	write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
	20	Boston, MA 02111-1307, USA. */
	21
	22	#include "config.h"
	23	#include <stdlib.h>
	24	#include <assert.h>
	25	#include "libgfortran.h"
	26
7d7b8bfe RH	27
	28	extern void __sum_r8 (gfc_array_r8 , gfc_array_r8 , index_type *);
	29	export_proto_np(__sum_r8);
	30
6de9cd9a	31	void
7d7b8bfe	32	__sum_r8 (gfc_array_r8 retarray, gfc_array_r8 array, index_type *pdim)
6de9cd9a DN	33	{
	34	index_type count[GFC_MAX_DIMENSIONS - 1];
	35	index_type extent[GFC_MAX_DIMENSIONS - 1];
	36	index_type sstride[GFC_MAX_DIMENSIONS - 1];
	37	index_type dstride[GFC_MAX_DIMENSIONS - 1];
	38	GFC_REAL_8 *base;
	39	GFC_REAL_8 *dest;
	40	index_type rank;
	41	index_type n;
	42	index_type len;
	43	index_type delta;
	44	index_type dim;
	45
	46	/* Make dim zero based to avoid confusion. */
	47	dim = (*pdim) - 1;
	48	rank = GFC_DESCRIPTOR_RANK (array) - 1;
	49	assert (rank == GFC_DESCRIPTOR_RANK (retarray));
	50	if (array->dim[0].stride == 0)
	51	array->dim[0].stride = 1;
	52	if (retarray->dim[0].stride == 0)
	53	retarray->dim[0].stride = 1;
	54
	55	len = array->dim[dim].ubound + 1 - array->dim[dim].lbound;
	56	delta = array->dim[dim].stride;
	57
	58	for (n = 0; n < dim; n++)
	59	{
	60	sstride[n] = array->dim[n].stride;
	61	extent[n] = array->dim[n].ubound + 1 - array->dim[n].lbound;
	62	}
	63	for (n = dim; n < rank; n++)
	64	{
	65	sstride[n] = array->dim[n + 1].stride;
	66	extent[n] =
	67	array->dim[n + 1].ubound + 1 - array->dim[n + 1].lbound;
	68	}
	69
6c167c45 VL	70	if (retarray->data == NULL)
	71	{
	72	for (n = 0; n < rank; n++)
	73	{
	74	retarray->dim[n].lbound = 0;
	75	retarray->dim[n].ubound = extent[n]-1;
	76	if (n == 0)
	77	retarray->dim[n].stride = 1;
	78	else
	79	retarray->dim[n].stride = retarray->dim[n-1].stride * extent[n-1];
	80	}
	81
07d3cebe RH	82	retarray->data
	83	= internal_malloc_size (sizeof (GFC_REAL_8)
	84	* retarray->dim[rank-1].stride
	85	* extent[rank-1]);
6c167c45 VL	86	retarray->base = 0;
	87	}
	88
6de9cd9a DN	89	for (n = 0; n < rank; n++)
	90	{
	91	count[n] = 0;
	92	dstride[n] = retarray->dim[n].stride;
	93	if (extent[n] <= 0)
	94	len = 0;
	95	}
	96
	97	base = array->data;
	98	dest = retarray->data;
	99
	100	while (base)
	101	{
	102	GFC_REAL_8 *src;
	103	GFC_REAL_8 result;
	104	src = base;
	105	{
	106
	107	result = 0;
	108	if (len <= 0)
	109	*dest = 0;
	110	else
	111	{
	112	for (n = 0; n < len; n++, src += delta)
	113	{
	114
	115	result += *src;
	116	}
	117	*dest = result;
	118	}
	119	}
	120	/* Advance to the next element. */
	121	count[0]++;
	122	base += sstride[0];
	123	dest += dstride[0];
	124	n = 0;
	125	while (count[n] == extent[n])
	126	{
	127	/* When we get to the end of a dimension, reset it and increment
	128	the next dimension. */
	129	count[n] = 0;
	130	/* We could precalculate these products, but this is a less
	131	frequently used path so proabably not worth it. */
	132	base -= sstride[n] * extent[n];
	133	dest -= dstride[n] * extent[n];
	134	n++;
	135	if (n == rank)
	136	{
	137	/* Break out of the look. */
	138	base = NULL;
	139	break;
	140	}
	141	else
	142	{
	143	count[n]++;
	144	base += sstride[n];
	145	dest += dstride[n];
	146	}
	147	}
	148	}
	149	}
	150
7d7b8bfe RH	151
	152	extern void __msum_r8 (gfc_array_r8 , gfc_array_r8 , index_type *,
	153	gfc_array_l4 *);
	154	export_proto_np(__msum_r8);
	155
6de9cd9a DN	156	void
	157	__msum_r8 (gfc_array_r8 * retarray, gfc_array_r8 * array, index_type pdim, gfc_array_l4 mask)
	158	{
	159	index_type count[GFC_MAX_DIMENSIONS - 1];
	160	index_type extent[GFC_MAX_DIMENSIONS - 1];
	161	index_type sstride[GFC_MAX_DIMENSIONS - 1];
	162	index_type dstride[GFC_MAX_DIMENSIONS - 1];
	163	index_type mstride[GFC_MAX_DIMENSIONS - 1];
	164	GFC_REAL_8 *dest;
	165	GFC_REAL_8 *base;
	166	GFC_LOGICAL_4 *mbase;
	167	int rank;
	168	int dim;
	169	index_type n;
	170	index_type len;
	171	index_type delta;
	172	index_type mdelta;
	173
	174	dim = (*pdim) - 1;
	175	rank = GFC_DESCRIPTOR_RANK (array) - 1;
	176	assert (rank == GFC_DESCRIPTOR_RANK (retarray));
	177	if (array->dim[0].stride == 0)
	178	array->dim[0].stride = 1;
	179	if (retarray->dim[0].stride == 0)
	180	retarray->dim[0].stride = 1;
	181
	182	len = array->dim[dim].ubound + 1 - array->dim[dim].lbound;
	183	if (len <= 0)
	184	return;
	185	delta = array->dim[dim].stride;
	186	mdelta = mask->dim[dim].stride;
	187
	188	for (n = 0; n < dim; n++)
	189	{
	190	sstride[n] = array->dim[n].stride;
	191	mstride[n] = mask->dim[n].stride;
	192	extent[n] = array->dim[n].ubound + 1 - array->dim[n].lbound;
	193	}
	194	for (n = dim; n < rank; n++)
	195	{
	196	sstride[n] = array->dim[n + 1].stride;
	197	mstride[n] = mask->dim[n + 1].stride;
	198	extent[n] =
	199	array->dim[n + 1].ubound + 1 - array->dim[n + 1].lbound;
	200	}
	201
	202	for (n = 0; n < rank; n++)
	203	{
	204	count[n] = 0;
	205	dstride[n] = retarray->dim[n].stride;
	206	if (extent[n] <= 0)
	207	return;
	208	}
	209
	210	dest = retarray->data;
	211	base = array->data;
	212	mbase = mask->data;
	213
	214	if (GFC_DESCRIPTOR_SIZE (mask) != 4)
	215	{
	216	/* This allows the same loop to be used for all logical types. */
	217	assert (GFC_DESCRIPTOR_SIZE (mask) == 8);
	218	for (n = 0; n < rank; n++)
	219	mstride[n] <<= 1;
220	mdelta <<= 1;
221	mbase = (GFOR_POINTER_L8_TO_L4 (mbase));
222	}
223
224	while (base)
225	{
226	GFC_REAL_8 *src;
227	GFC_LOGICAL_4 *msrc;
228	GFC_REAL_8 result;
229	src = base;
230	msrc = mbase;
231	{
232
233	result = 0;
234	if (len <= 0)
235	*dest = 0;
236	else
237	{
238	for (n = 0; n < len; n++, src += delta, msrc += mdelta)
239	{
240
241	if (*msrc)
242	result += *src;
243	}
244	*dest = result;
245	}
246	}
247	/* Advance to the next element. */
248	count[0]++;
249	base += sstride[0];
250	mbase += mstride[0];
251	dest += dstride[0];
252	n = 0;
253	while (count[n] == extent[n])
254	{
255	/* When we get to the end of a dimension, reset it and increment
256	the next dimension. */
257	count[n] = 0;
258	/* We could precalculate these products, but this is a less
259	frequently used path so proabably not worth it. */
260	base -= sstride[n] * extent[n];
261	mbase -= mstride[n] * extent[n];
262	dest -= dstride[n] * extent[n];
263	n++;
264	if (n == rank)
265	{
266	/* Break out of the look. */
267	base = NULL;
268	break;
269	}
270	else
271	{
272	count[n]++;
273	base += sstride[n];
274	mbase += mstride[n];
275	dest += dstride[n];
276	}
277	}
278	}
279	}