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

/* Implementation of the PRODUCT 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"

void
__product_c4 (gfc_array_c4 * retarray, gfc_array_c4 *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_COMPLEX_4 *base;
  GFC_COMPLEX_4 *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 (sizeof (GFC_COMPLEX_4) * 
                                        (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_COMPLEX_4 *src;
      GFC_COMPLEX_4 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 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];
            }
        }
    }
}

void
__mproduct_c4 (gfc_array_c4 * retarray, gfc_array_c4 * 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_COMPLEX_4 *dest;
  GFC_COMPLEX_4 *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_COMPLEX_4 *src;
      GFC_LOGICAL_4 *msrc;
      GFC_COMPLEX_4 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 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 PRODUCT 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
6de9cd9a DN	27	void
	28	__product_c4 (gfc_array_c4 * retarray, gfc_array_c4 array, index_type pdim)
	29	{
	30	index_type count[GFC_MAX_DIMENSIONS - 1];
	31	index_type extent[GFC_MAX_DIMENSIONS - 1];
	32	index_type sstride[GFC_MAX_DIMENSIONS - 1];
	33	index_type dstride[GFC_MAX_DIMENSIONS - 1];
	34	GFC_COMPLEX_4 *base;
	35	GFC_COMPLEX_4 *dest;
	36	index_type rank;
	37	index_type n;
	38	index_type len;
	39	index_type delta;
	40	index_type dim;
	41
	42	/* Make dim zero based to avoid confusion. */
	43	dim = (*pdim) - 1;
	44	rank = GFC_DESCRIPTOR_RANK (array) - 1;
	45	assert (rank == GFC_DESCRIPTOR_RANK (retarray));
	46	if (array->dim[0].stride == 0)
	47	array->dim[0].stride = 1;
	48	if (retarray->dim[0].stride == 0)
	49	retarray->dim[0].stride = 1;
	50
	51	len = array->dim[dim].ubound + 1 - array->dim[dim].lbound;
	52	delta = array->dim[dim].stride;
	53
	54	for (n = 0; n < dim; n++)
	55	{
	56	sstride[n] = array->dim[n].stride;
	57	extent[n] = array->dim[n].ubound + 1 - array->dim[n].lbound;
	58	}
	59	for (n = dim; n < rank; n++)
	60	{
	61	sstride[n] = array->dim[n + 1].stride;
	62	extent[n] =
	63	array->dim[n + 1].ubound + 1 - array->dim[n + 1].lbound;
	64	}
	65
6c167c45 VL	66	if (retarray->data == NULL)
	67	{
	68	for (n = 0; n < rank; n++)
	69	{
	70	retarray->dim[n].lbound = 0;
	71	retarray->dim[n].ubound = extent[n]-1;
	72	if (n == 0)
	73	retarray->dim[n].stride = 1;
	74	else
	75	retarray->dim[n].stride = retarray->dim[n-1].stride * extent[n-1];
	76	}
	77
	78	retarray->data = internal_malloc (sizeof (GFC_COMPLEX_4) *
	79	(retarray->dim[rank-1].stride * extent[rank-1]));
	80	retarray->base = 0;
	81	}
	82
6de9cd9a DN	83	for (n = 0; n < rank; n++)
	84	{
	85	count[n] = 0;
	86	dstride[n] = retarray->dim[n].stride;
	87	if (extent[n] <= 0)
	88	len = 0;
	89	}
	90
	91	base = array->data;
	92	dest = retarray->data;
	93
	94	while (base)
	95	{
	96	GFC_COMPLEX_4 *src;
	97	GFC_COMPLEX_4 result;
	98	src = base;
	99	{
	100
	101	result = 1;
	102	if (len <= 0)
	103	*dest = 1;
	104	else
	105	{
	106	for (n = 0; n < len; n++, src += delta)
	107	{
	108
	109	result = src;
	110	}
	111	*dest = result;
	112	}
	113	}
	114	/* Advance to the next element. */
	115	count[0]++;
	116	base += sstride[0];
	117	dest += dstride[0];
	118	n = 0;
	119	while (count[n] == extent[n])
	120	{
	121	/* When we get to the end of a dimension, reset it and increment
	122	the next dimension. */
	123	count[n] = 0;
	124	/* We could precalculate these products, but this is a less
	125	frequently used path so proabably not worth it. */
	126	base -= sstride[n] * extent[n];
	127	dest -= dstride[n] * extent[n];
	128	n++;
	129	if (n == rank)
	130	{
	131	/* Break out of the look. */
	132	base = NULL;
	133	break;
	134	}
	135	else
	136	{
	137	count[n]++;
	138	base += sstride[n];
	139	dest += dstride[n];
	140	}
	141	}
	142	}
	143	}
	144
	145	void
	146	__mproduct_c4 (gfc_array_c4 * retarray, gfc_array_c4 * array, index_type pdim, gfc_array_l4 mask)
147	{
148	index_type count[GFC_MAX_DIMENSIONS - 1];
149	index_type extent[GFC_MAX_DIMENSIONS - 1];
150	index_type sstride[GFC_MAX_DIMENSIONS - 1];
151	index_type dstride[GFC_MAX_DIMENSIONS - 1];
152	index_type mstride[GFC_MAX_DIMENSIONS - 1];
153	GFC_COMPLEX_4 *dest;
154	GFC_COMPLEX_4 *base;
155	GFC_LOGICAL_4 *mbase;
156	int rank;
157	int dim;
158	index_type n;
159	index_type len;
160	index_type delta;
161	index_type mdelta;
162
163	dim = (*pdim) - 1;
164	rank = GFC_DESCRIPTOR_RANK (array) - 1;
165	assert (rank == GFC_DESCRIPTOR_RANK (retarray));
166	if (array->dim[0].stride == 0)
167	array->dim[0].stride = 1;
168	if (retarray->dim[0].stride == 0)
169	retarray->dim[0].stride = 1;
170
171	len = array->dim[dim].ubound + 1 - array->dim[dim].lbound;
172	if (len <= 0)
173	return;
174	delta = array->dim[dim].stride;
175	mdelta = mask->dim[dim].stride;
176
177	for (n = 0; n < dim; n++)
178	{
179	sstride[n] = array->dim[n].stride;
180	mstride[n] = mask->dim[n].stride;
181	extent[n] = array->dim[n].ubound + 1 - array->dim[n].lbound;
182	}
183	for (n = dim; n < rank; n++)
184	{
185	sstride[n] = array->dim[n + 1].stride;
186	mstride[n] = mask->dim[n + 1].stride;
187	extent[n] =
188	array->dim[n + 1].ubound + 1 - array->dim[n + 1].lbound;
189	}
190
191	for (n = 0; n < rank; n++)
192	{
193	count[n] = 0;
194	dstride[n] = retarray->dim[n].stride;
195	if (extent[n] <= 0)
196	return;
197	}
198
199	dest = retarray->data;
200	base = array->data;
201	mbase = mask->data;
202
203	if (GFC_DESCRIPTOR_SIZE (mask) != 4)
204	{
205	/* This allows the same loop to be used for all logical types. */
206	assert (GFC_DESCRIPTOR_SIZE (mask) == 8);
207	for (n = 0; n < rank; n++)
208	mstride[n] <<= 1;
209	mdelta <<= 1;
210	mbase = (GFOR_POINTER_L8_TO_L4 (mbase));
211	}
212
213	while (base)
214	{
215	GFC_COMPLEX_4 *src;
216	GFC_LOGICAL_4 *msrc;
217	GFC_COMPLEX_4 result;
218	src = base;
219	msrc = mbase;
220	{
221
222	result = 1;
223	if (len <= 0)
224	*dest = 1;
225	else
226	{
227	for (n = 0; n < len; n++, src += delta, msrc += mdelta)
228	{
229
230	if (*msrc)
231	result = src;
232	}
233	*dest = result;
234	}
235	}
236	/* Advance to the next element. */
237	count[0]++;
238	base += sstride[0];
239	mbase += mstride[0];
240	dest += dstride[0];
241	n = 0;
242	while (count[n] == extent[n])
243	{
244	/* When we get to the end of a dimension, reset it and increment
245	the next dimension. */
246	count[n] = 0;
247	/* We could precalculate these products, but this is a less
248	frequently used path so proabably not worth it. */
249	base -= sstride[n] * extent[n];
250	mbase -= mstride[n] * extent[n];
251	dest -= dstride[n] * extent[n];
252	n++;
253	if (n == rank)
254	{
255	/* Break out of the look. */
256	base = NULL;
257	break;
258	}
259	else
260	{
261	count[n]++;
262	base += sstride[n];
263	mbase += mstride[n];
264	dest += dstride[n];
265	}
266	}
267	}
268	}
269