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

/* Implementation of the PRODUCT intrinsic
   Copyright 2002, 2007, 2009 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"
#include <stdlib.h>
#include <assert.h>


#if defined (HAVE_GFC_REAL_10) && defined (HAVE_GFC_REAL_10)


extern void product_r10 (gfc_array_r10 * const restrict, 
	gfc_array_r10 * const restrict, const index_type * const restrict);
export_proto(product_r10);

void
product_r10 (gfc_array_r10 * const restrict retarray, 
	gfc_array_r10 * 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_10 * restrict base;
  GFC_REAL_10 * restrict dest;
  index_type rank;
  index_type n;
  index_type len;
  index_type delta;
  index_type dim;
  int continue_loop;

  /* Make dim zero based to avoid confusion.  */
  dim = (*pdim) - 1;
  rank = GFC_DESCRIPTOR_RANK (array) - 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->data == 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;
      retarray->dtype = (array->dtype & ~GFC_DTYPE_RANK_MASK) | rank;

      alloc_size = sizeof (GFC_REAL_10) * 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->data = internal_malloc_size (alloc_size);
    }
  else
    {
      if (rank != GFC_DESCRIPTOR_RANK (retarray))
	runtime_error ("rank of return array incorrect in"
		       " PRODUCT 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", "PRODUCT");
    }

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

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

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

  result = 1;
        if (len <= 0)
	  *dest = 1;
	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.  */
	      continue_loop = 0;
	      break;
            }
          else
            {
              count[n]++;
              base += sstride[n];
              dest += dstride[n];
            }
        }
    }
}


extern void mproduct_r10 (gfc_array_r10 * const restrict, 
	gfc_array_r10 * const restrict, const index_type * const restrict,
	gfc_array_l1 * const restrict);
export_proto(mproduct_r10);

void
mproduct_r10 (gfc_array_r10 * const restrict retarray, 
	gfc_array_r10 * 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_10 * restrict dest;
  const GFC_REAL_10 * 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 = GFC_DESCRIPTOR_EXTENT(array,dim);
  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 = 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->data == 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 = sizeof (GFC_REAL_10) * GFC_DESCRIPTOR_STRIDE(retarray,rank-1)
    		   * 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.  */
	  GFC_DIMENSION_SET(retarray->dim[0], 0, -1, 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 PRODUCT intrinsic");

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

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

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

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

  result = 1;
        if (len <= 0)
	  *dest = 1;
	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 sproduct_r10 (gfc_array_r10 * const restrict, 
	gfc_array_r10 * const restrict, const index_type * const restrict,
	GFC_LOGICAL_4 *);
export_proto(sproduct_r10);

void
sproduct_r10 (gfc_array_r10 * const restrict retarray, 
	gfc_array_r10 * 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_10 * restrict dest;
  index_type rank;
  index_type n;
  index_type dim;


  if (*mask)
    {
      product_r10 (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] = 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->data == 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;
      retarray->dtype = (array->dtype & ~GFC_DTYPE_RANK_MASK) | rank;

      alloc_size = sizeof (GFC_REAL_10) * 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->data = internal_malloc_size (alloc_size);
    }
  else
    {
      if (rank != GFC_DESCRIPTOR_RANK (retarray))
	runtime_error ("rank of return array incorrect in"
		       " PRODUCT 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"
			       " PRODUCT 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->data;

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