1/* Implementation of the RESHAPE intrinsic
2   Copyright (C) 2002-2022 Free Software Foundation, Inc.
3   Contributed by Paul Brook <paul@nowt.org>
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#if defined (HAVE_GFC_COMPLEX_4)
30
31typedef GFC_FULL_ARRAY_DESCRIPTOR(1, index_type) shape_type;
32
33
34extern void reshape_c4 (gfc_array_c4 * const restrict,
35	gfc_array_c4 * const restrict,
36	shape_type * const restrict,
37	gfc_array_c4 * const restrict,
38	shape_type * const restrict);
39export_proto(reshape_c4);
40
41void
42reshape_c4 (gfc_array_c4 * const restrict ret,
43	gfc_array_c4 * const restrict source,
44	shape_type * const restrict shape,
45	gfc_array_c4 * const restrict pad,
46	shape_type * const restrict order)
47{
48  /* r.* indicates the return array.  */
49  index_type rcount[GFC_MAX_DIMENSIONS];
50  index_type rextent[GFC_MAX_DIMENSIONS];
51  index_type rstride[GFC_MAX_DIMENSIONS];
52  index_type rstride0;
53  index_type rdim;
54  index_type rsize;
55  index_type rs;
56  index_type rex;
57  GFC_COMPLEX_4 *rptr;
58  /* s.* indicates the source array.  */
59  index_type scount[GFC_MAX_DIMENSIONS];
60  index_type sextent[GFC_MAX_DIMENSIONS];
61  index_type sstride[GFC_MAX_DIMENSIONS];
62  index_type sstride0;
63  index_type sdim;
64  index_type ssize;
65  const GFC_COMPLEX_4 *sptr;
66  /* p.* indicates the pad array.  */
67  index_type pcount[GFC_MAX_DIMENSIONS];
68  index_type pextent[GFC_MAX_DIMENSIONS];
69  index_type pstride[GFC_MAX_DIMENSIONS];
70  index_type pdim;
71  index_type psize;
72  const GFC_COMPLEX_4 *pptr;
73
74  const GFC_COMPLEX_4 *src;
75  int sempty, pempty, shape_empty;
76  index_type shape_data[GFC_MAX_DIMENSIONS];
77
78  rdim = GFC_DESCRIPTOR_EXTENT(shape,0);
79  /* rdim is always > 0; this lets the compiler optimize more and
80   avoids a potential warning.  */
81  GFC_ASSERT(rdim>0);
82
83  if (rdim != GFC_DESCRIPTOR_RANK(ret))
84    runtime_error("rank of return array incorrect in RESHAPE intrinsic");
85
86  shape_empty = 0;
87
88  for (index_type n = 0; n < rdim; n++)
89    {
90      shape_data[n] = shape->base_addr[n * GFC_DESCRIPTOR_STRIDE(shape,0)];
91      if (shape_data[n] <= 0)
92      {
93        shape_data[n] = 0;
94	shape_empty = 1;
95      }
96    }
97
98  if (ret->base_addr == NULL)
99    {
100      index_type alloc_size;
101
102      rs = 1;
103      for (index_type n = 0; n < rdim; n++)
104	{
105	  rex = shape_data[n];
106
107	  GFC_DIMENSION_SET(ret->dim[n], 0, rex - 1, rs);
108
109	  rs *= rex;
110	}
111      ret->offset = 0;
112
113      if (unlikely (rs < 1))
114        alloc_size = 0;
115      else
116        alloc_size = rs;
117
118      ret->base_addr = xmallocarray (alloc_size, sizeof (GFC_COMPLEX_4));
119      ret->dtype.rank = rdim;
120    }
121
122  if (shape_empty)
123    return;
124
125  if (pad)
126    {
127      pdim = GFC_DESCRIPTOR_RANK (pad);
128      psize = 1;
129      pempty = 0;
130      for (index_type n = 0; n < pdim; n++)
131        {
132          pcount[n] = 0;
133          pstride[n] = GFC_DESCRIPTOR_STRIDE(pad,n);
134          pextent[n] = GFC_DESCRIPTOR_EXTENT(pad,n);
135          if (pextent[n] <= 0)
136	    {
137	      pempty = 1;
138	      pextent[n] = 0;
139	    }
140
141          if (psize == pstride[n])
142            psize *= pextent[n];
143          else
144            psize = 0;
145        }
146      pptr = pad->base_addr;
147    }
148  else
149    {
150      pdim = 0;
151      psize = 1;
152      pempty = 1;
153      pptr = NULL;
154    }
155
156  if (unlikely (compile_options.bounds_check))
157    {
158      index_type ret_extent, source_extent;
159
160      rs = 1;
161      for (index_type n = 0; n < rdim; n++)
162	{
163	  rs *= shape_data[n];
164	  ret_extent = GFC_DESCRIPTOR_EXTENT(ret,n);
165	  if (ret_extent != shape_data[n])
166	    runtime_error("Incorrect extent in return value of RESHAPE"
167			  " intrinsic in dimension %ld: is %ld,"
168			  " should be %ld", (long int) n+1,
169			  (long int) ret_extent, (long int) shape_data[n]);
170	}
171
172      source_extent = 1;
173      sdim = GFC_DESCRIPTOR_RANK (source);
174      for (index_type n = 0; n < sdim; n++)
175	{
176	  index_type se;
177	  se = GFC_DESCRIPTOR_EXTENT(source,n);
178	  source_extent *= se > 0 ? se : 0;
179	}
180
181      if (rs > source_extent && (!pad || pempty))
182	runtime_error("Incorrect size in SOURCE argument to RESHAPE"
183		      " intrinsic: is %ld, should be %ld",
184		      (long int) source_extent, (long int) rs);
185
186      if (order)
187	{
188	  int seen[GFC_MAX_DIMENSIONS];
189	  index_type v;
190
191	  for (index_type n = 0; n < rdim; n++)
192	    seen[n] = 0;
193
194	  for (index_type n = 0; n < rdim; n++)
195	    {
196	      v = order->base_addr[n * GFC_DESCRIPTOR_STRIDE(order,0)] - 1;
197
198	      if (v < 0 || v >= rdim)
199		runtime_error("Value %ld out of range in ORDER argument"
200			      " to RESHAPE intrinsic", (long int) v + 1);
201
202	      if (seen[v] != 0)
203		runtime_error("Duplicate value %ld in ORDER argument to"
204			      " RESHAPE intrinsic", (long int) v + 1);
205
206	      seen[v] = 1;
207	    }
208	}
209    }
210
211  rsize = 1;
212  for (index_type n = 0; n < rdim; n++)
213    {
214      index_type dim;
215      if (order)
216        dim = order->base_addr[n * GFC_DESCRIPTOR_STRIDE(order,0)] - 1;
217      else
218        dim = n;
219
220      rcount[n] = 0;
221      rstride[n] = GFC_DESCRIPTOR_STRIDE(ret,dim);
222      rextent[n] = GFC_DESCRIPTOR_EXTENT(ret,dim);
223      if (rextent[n] < 0)
224        rextent[n] = 0;
225
226      if (rextent[n] != shape_data[dim])
227        runtime_error ("shape and target do not conform");
228
229      if (rsize == rstride[n])
230        rsize *= rextent[n];
231      else
232        rsize = 0;
233      if (rextent[n] <= 0)
234        return;
235    }
236
237  sdim = GFC_DESCRIPTOR_RANK (source);
238
239  /* sdim is always > 0; this lets the compiler optimize more and
240   avoids a warning.  */
241  GFC_ASSERT(sdim>0);
242
243  ssize = 1;
244  sempty = 0;
245  for (index_type n = 0; n < sdim; n++)
246    {
247      scount[n] = 0;
248      sstride[n] = GFC_DESCRIPTOR_STRIDE(source,n);
249      sextent[n] = GFC_DESCRIPTOR_EXTENT(source,n);
250      if (sextent[n] <= 0)
251	{
252	  sempty = 1;
253	  sextent[n] = 0;
254	}
255
256      if (ssize == sstride[n])
257        ssize *= sextent[n];
258      else
259        ssize = 0;
260    }
261
262  if (rsize != 0 && ssize != 0 && psize != 0)
263    {
264      rsize *= sizeof (GFC_COMPLEX_4);
265      ssize *= sizeof (GFC_COMPLEX_4);
266      psize *= sizeof (GFC_COMPLEX_4);
267      reshape_packed ((char *)ret->base_addr, rsize, (char *)source->base_addr,
268		      ssize, pad ? (char *)pad->base_addr : NULL, psize);
269      return;
270    }
271  rptr = ret->base_addr;
272  src = sptr = source->base_addr;
273  rstride0 = rstride[0];
274  sstride0 = sstride[0];
275
276  if (sempty && pempty)
277    abort ();
278
279  if (sempty)
280    {
281      /* Pretend we are using the pad array the first time around, too.  */
282      src = pptr;
283      sptr = pptr;
284      sdim = pdim;
285      for (index_type dim = 0; dim < pdim; dim++)
286	{
287	  scount[dim] = pcount[dim];
288	  sextent[dim] = pextent[dim];
289	  sstride[dim] = pstride[dim];
290	  sstride0 = pstride[0];
291	}
292    }
293
294  while (rptr)
295    {
296      /* Select between the source and pad arrays.  */
297      *rptr = *src;
298      /* Advance to the next element.  */
299      rptr += rstride0;
300      src += sstride0;
301      rcount[0]++;
302      scount[0]++;
303
304      /* Advance to the next destination element.  */
305      index_type n = 0;
306      while (rcount[n] == rextent[n])
307        {
308          /* When we get to the end of a dimension, reset it and increment
309             the next dimension.  */
310          rcount[n] = 0;
311          /* We could precalculate these products, but this is a less
312             frequently used path so probably not worth it.  */
313          rptr -= rstride[n] * rextent[n];
314          n++;
315          if (n == rdim)
316            {
317              /* Break out of the loop.  */
318              rptr = NULL;
319              break;
320            }
321          else
322            {
323              rcount[n]++;
324              rptr += rstride[n];
325            }
326        }
327      /* Advance to the next source element.  */
328      n = 0;
329      while (scount[n] == sextent[n])
330        {
331          /* When we get to the end of a dimension, reset it and increment
332             the next dimension.  */
333          scount[n] = 0;
334          /* We could precalculate these products, but this is a less
335             frequently used path so probably not worth it.  */
336          src -= sstride[n] * sextent[n];
337          n++;
338          if (n == sdim)
339            {
340              if (sptr && pad)
341                {
342                  /* Switch to the pad array.  */
343                  sptr = NULL;
344                  sdim = pdim;
345                  for (index_type dim = 0; dim < pdim; dim++)
346                    {
347                      scount[dim] = pcount[dim];
348                      sextent[dim] = pextent[dim];
349                      sstride[dim] = pstride[dim];
350                      sstride0 = sstride[0];
351                    }
352                }
353              /* We now start again from the beginning of the pad array.  */
354              src = pptr;
355              break;
356            }
357          else
358            {
359              scount[n]++;
360              src += sstride[n];
361            }
362        }
363    }
364}
365
366#endif
367