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