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