1/* Implementation of the NORM2 intrinsic
2   Copyright (C) 2010-2020 Free Software Foundation, Inc.
3   Contributed by Tobias Burnus  <burnus@net-b.de>
4
5This file is part of the GNU Fortran runtime library (libgfortran).
6
7Libgfortran is free software; you can redistribute it and/or
8modify it under the terms of the GNU General Public
9License as published by the Free Software Foundation; either
10version 3 of the License, or (at your option) any later version.
11
12Libgfortran is distributed in the hope that it will be useful,
13but WITHOUT ANY WARRANTY; without even the implied warranty of
14MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15GNU General Public License for more details.
16
17Under Section 7 of GPL version 3, you are granted additional
18permissions described in the GCC Runtime Library Exception, version
193.1, as published by the Free Software Foundation.
20
21You should have received a copy of the GNU General Public License and
22a copy of the GCC Runtime Library Exception along with this program;
23see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
24<http://www.gnu.org/licenses/>.  */
25
26#include "libgfortran.h"
27
28
29
30#if defined (HAVE_GFC_REAL_16) && defined (HAVE_GFC_REAL_16) && (defined(GFC_REAL_16_IS_FLOAT128) || defined(HAVE_SQRTL)) && (defined(GFC_REAL_16_IS_FLOAT128) || defined(HAVE_FABSL))
31
32#if defined(GFC_REAL_16_IS_FLOAT128)
33#define MATHFUNC(funcname) funcname ## q
34#else
35#define MATHFUNC(funcname) funcname ## l
36#endif
37
38
39extern void norm2_r16 (gfc_array_r16 * const restrict,
40	gfc_array_r16 * const restrict, const index_type * const restrict);
41export_proto(norm2_r16);
42
43void
44norm2_r16 (gfc_array_r16 * const restrict retarray,
45	gfc_array_r16 * const restrict array,
46	const index_type * const restrict pdim)
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];
52  const GFC_REAL_16 * restrict base;
53  GFC_REAL_16 * restrict dest;
54  index_type rank;
55  index_type n;
56  index_type len;
57  index_type delta;
58  index_type dim;
59  int continue_loop;
60
61  /* Make dim zero based to avoid confusion.  */
62  rank = GFC_DESCRIPTOR_RANK (array) - 1;
63  dim = (*pdim) - 1;
64
65  if (unlikely (dim < 0 || dim > rank))
66    {
67      runtime_error ("Dim argument incorrect in NORM intrinsic: "
68 		     "is %ld, should be between 1 and %ld",
69		     (long int) dim + 1, (long int) rank + 1);
70    }
71
72  len = GFC_DESCRIPTOR_EXTENT(array,dim);
73  if (len < 0)
74    len = 0;
75  delta = GFC_DESCRIPTOR_STRIDE(array,dim);
76
77  for (n = 0; n < dim; n++)
78    {
79      sstride[n] = GFC_DESCRIPTOR_STRIDE(array,n);
80      extent[n] = GFC_DESCRIPTOR_EXTENT(array,n);
81
82      if (extent[n] < 0)
83	extent[n] = 0;
84    }
85  for (n = dim; n < rank; n++)
86    {
87      sstride[n] = GFC_DESCRIPTOR_STRIDE(array, n + 1);
88      extent[n] = GFC_DESCRIPTOR_EXTENT(array, n + 1);
89
90      if (extent[n] < 0)
91	extent[n] = 0;
92    }
93
94  if (retarray->base_addr == NULL)
95    {
96      size_t alloc_size, str;
97
98      for (n = 0; n < rank; n++)
99	{
100	  if (n == 0)
101	    str = 1;
102	  else
103	    str = GFC_DESCRIPTOR_STRIDE(retarray,n-1) * extent[n-1];
104
105	  GFC_DIMENSION_SET(retarray->dim[n], 0, extent[n] - 1, str);
106
107	}
108
109      retarray->offset = 0;
110      retarray->dtype.rank = rank;
111
112      alloc_size = GFC_DESCRIPTOR_STRIDE(retarray,rank-1) * extent[rank-1];
113
114      retarray->base_addr = xmallocarray (alloc_size, sizeof (GFC_REAL_16));
115      if (alloc_size == 0)
116	{
117	  /* Make sure we have a zero-sized array.  */
118	  GFC_DIMENSION_SET(retarray->dim[0], 0, -1, 1);
119	  return;
120
121	}
122    }
123  else
124    {
125      if (rank != GFC_DESCRIPTOR_RANK (retarray))
126	runtime_error ("rank of return array incorrect in"
127		       " NORM intrinsic: is %ld, should be %ld",
128		       (long int) (GFC_DESCRIPTOR_RANK (retarray)),
129		       (long int) rank);
130
131      if (unlikely (compile_options.bounds_check))
132	bounds_ifunction_return ((array_t *) retarray, extent,
133				 "return value", "NORM");
134    }
135
136  for (n = 0; n < rank; n++)
137    {
138      count[n] = 0;
139      dstride[n] = GFC_DESCRIPTOR_STRIDE(retarray,n);
140      if (extent[n] <= 0)
141	return;
142    }
143
144  base = array->base_addr;
145  dest = retarray->base_addr;
146
147  continue_loop = 1;
148  while (continue_loop)
149    {
150      const GFC_REAL_16 * restrict src;
151      GFC_REAL_16 result;
152      src = base;
153      {
154
155	GFC_REAL_16 scale;
156	result = 0;
157	scale = 1;
158	if (len <= 0)
159	  *dest = 0;
160	else
161	  {
162#if ! defined HAVE_BACK_ARG
163	    for (n = 0; n < len; n++, src += delta)
164	      {
165#endif
166
167	  if (*src != 0)
168	    {
169	      GFC_REAL_16 absX, val;
170	      absX = MATHFUNC(fabs) (*src);
171	      if (scale < absX)
172		{
173		  val = scale / absX;
174		  result = 1 + result * val * val;
175		  scale = absX;
176		}
177	      else
178		{
179		  val = absX / scale;
180		  result += val * val;
181		}
182	    }
183	      }
184	    result = scale * MATHFUNC(sqrt) (result);
185	    *dest = result;
186	  }
187      }
188      /* Advance to the next element.  */
189      count[0]++;
190      base += sstride[0];
191      dest += dstride[0];
192      n = 0;
193      while (count[n] == extent[n])
194	{
195	  /* When we get to the end of a dimension, reset it and increment
196	     the next dimension.  */
197	  count[n] = 0;
198	  /* We could precalculate these products, but this is a less
199	     frequently used path so probably not worth it.  */
200	  base -= sstride[n] * extent[n];
201	  dest -= dstride[n] * extent[n];
202	  n++;
203	  if (n >= rank)
204	    {
205	      /* Break out of the loop.  */
206	      continue_loop = 0;
207	      break;
208	    }
209	  else
210	    {
211	      count[n]++;
212	      base += sstride[n];
213	      dest += dstride[n];
214	    }
215	}
216    }
217}
218
219#endif
220