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