Subversion Repositories openrisc

/* Implementation of the NORM2 intrinsic

   Copyright 2010 Free Software Foundation, Inc.

   Contributed by Tobias Burnus  <burnus@net-b.de>

This file is part of the GNU Fortran 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 <math.h>

#include <assert.h>

#if defined (HAVE_GFC_REAL_8) && defined (HAVE_GFC_REAL_8) && defined (HAVE_SQRT) && defined (HAVE_FABS)

#define MATHFUNC(funcname) funcname

extern void norm2_r8 (gfc_array_r8 * const restrict,

        gfc_array_r8 * const restrict, const index_type * const restrict);

export_proto(norm2_r8);

void

norm2_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;

  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_8) * 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"

                       " NORM 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", "NORM");

    }

  for (n = 0; n < rank; n++)

    {

      count[n] = 0;

      dstride[n] = GFC_DESCRIPTOR_STRIDE(retarray,n);

      if (extent[n] <= 0)

        return;

    }

  base = array->data;

  dest = retarray->data;

  continue_loop = 1;

  while (continue_loop)

    {

      const GFC_REAL_8 * restrict src;

      GFC_REAL_8 result;

      src = base;

      {

        GFC_REAL_8 scale;

        result = 0;

        scale = 1;

        if (len <= 0)

          *dest = 0;

        else

          {

            for (n = 0; n < len; n++, src += delta)

              {

          if (*src != 0)

            {

              GFC_REAL_8 absX, val;

              absX = MATHFUNC(fabs) (*src);

              if (scale < absX)

                {

                  val = scale / absX;

                  result = 1 + result * val * val;

                  scale = absX;

                }

              else

                {

                  val = absX / scale;

                  result += val * val;

                }

            }

              }

            result = scale * MATHFUNC(sqrt) (result);

            *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];

            }

        }

    }

}

#endif