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

/* Implementation of the PRODUCT 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_COMPLEX_10) && defined (HAVE_GFC_COMPLEX_10)


extern void product_c10 (gfc_array_c10 * const restrict, 
	gfc_array_c10 * const restrict, const index_type * const restrict);
export_proto(product_c10);

void
product_c10 (gfc_array_c10 * const restrict retarray, 
	gfc_array_c10 * 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_COMPLEX_10 * restrict base;
  GFC_COMPLEX_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 = array->dim[dim].ubound + 1 - array->dim[dim].lbound;
  if (len < 0)
    len = 0;
  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_COMPLEX_10) * 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"
		       " PRODUCT 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"
			       " 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] = retarray->dim[n].stride;
      if (extent[n] <= 0)
        len = 0;
    }

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

  continue_loop = 1;
  while (continue_loop)
    {
      const GFC_COMPLEX_10 * restrict src;
      GFC_COMPLEX_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_c10 (gfc_array_c10 * const restrict, 
	gfc_array_c10 * const restrict, const index_type * const restrict,
	gfc_array_l1 * const restrict);
export_proto(mproduct_c10);

void
mproduct_c10 (gfc_array_c10 * const restrict retarray, 
	gfc_array_c10 * 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_COMPLEX_10 * restrict dest;
  const GFC_COMPLEX_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 = 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_COMPLEX_10) * 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 PRODUCT 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"
			       " 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++)
            {
              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"
			       " PRODUCT 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_COMPLEX_10 * restrict src;
      const GFC_LOGICAL_1 * restrict msrc;
      GFC_COMPLEX_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_c10 (gfc_array_c10 * const restrict, 
	gfc_array_c10 * const restrict, const index_type * const restrict,
	GFC_LOGICAL_4 *);
export_proto(sproduct_c10);

void
sproduct_c10 (gfc_array_c10 * const restrict retarray, 
	gfc_array_c10 * 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_COMPLEX_10 * restrict dest;
  index_type rank;
  index_type n;
  index_type dim;


  if (*mask)
    {
      product_c10 (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_COMPLEX_10) * 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"
		       " PRODUCT 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"
			       " 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] = retarray->dim[n].stride;
    }

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