libgfortran/m4/matmull.m4

   1 `/* Implementation of the MATMUL intrinsic
   2    Copyright (C) 2002-2024 Free Software Foundation, Inc.
   3    Contributed by Paul Brook <paul@nowt.org>
   4
   5 This file is part of the GNU Fortran runtime library (libgfortran).
   6
   7 Libgfortran is free software; you can redistribute it and/or
   8 modify it under the terms of the GNU General Public
   9 License as published by the Free Software Foundation; either
  10 version 3 of the License, or (at your option) any later version.
  11
  12 Libgfortran is distributed in the hope that it will be useful,
  13 but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 GNU General Public License for more details.
  16
  17 Under Section 7 of GPL version 3, you are granted additional
  18 permissions described in the GCC Runtime Library Exception, version
  19 3.1, as published by the Free Software Foundation.
  20
  21 You should have received a copy of the GNU General Public License and
  22 a copy of the GCC Runtime Library Exception along with this program;
  23 see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
  24 <http://www.gnu.org/licenses/>.  */
  25
  26 #include "libgfortran.h"
  27 #include <assert.h>'
  28
  29 include(iparm.m4)dnl
  30
  31 `#if defined (HAVE_'rtype_name`)
  32
  33 /* Dimensions: retarray(x,y) a(x, count) b(count,y).
  34    Either a or b can be rank 1.  In this case x or y is 1.  */
  35
  36 extern void matmul_'rtype_code` ('rtype` * const restrict,
  37         gfc_array_l1 * const restrict, gfc_array_l1 * const restrict);
  38 export_proto(matmul_'rtype_code`);
  39
  40 void
  41 matmul_'rtype_code` ('rtype` * const restrict retarray,
  42         gfc_array_l1 * const restrict a, gfc_array_l1 * const restrict b)
  43 {
  44   const GFC_LOGICAL_1 * restrict abase;
  45   const GFC_LOGICAL_1 * restrict bbase;
  46   'rtype_name` * restrict dest;
  47   index_type rxstride;
  48   index_type rystride;
  49   index_type xcount;
  50   index_type ycount;
  51   index_type xstride;
  52   index_type ystride;
  53   index_type x;
  54   index_type y;
  55   int a_kind;
  56   int b_kind;
  57
  58   const GFC_LOGICAL_1 * restrict pa;
  59   const GFC_LOGICAL_1 * restrict pb;
  60   index_type astride;
  61   index_type bstride;
  62   index_type count;
  63   index_type n;
  64
  65   assert (GFC_DESCRIPTOR_RANK (a) == 2
  66           || GFC_DESCRIPTOR_RANK (b) == 2);
  67
  68   if (retarray->base_addr == NULL)
  69     {
  70       if (GFC_DESCRIPTOR_RANK (a) == 1)
  71         {
  72           GFC_DIMENSION_SET(retarray->dim[0], 0,
  73                             GFC_DESCRIPTOR_EXTENT(b,1) - 1, 1);
  74         }
  75       else if (GFC_DESCRIPTOR_RANK (b) == 1)
  76         {
  77           GFC_DIMENSION_SET(retarray->dim[0], 0,
  78                             GFC_DESCRIPTOR_EXTENT(a,0) - 1, 1);
  79         }
  80       else
  81         {
  82           GFC_DIMENSION_SET(retarray->dim[0], 0,
  83                             GFC_DESCRIPTOR_EXTENT(a,0) - 1, 1);
  84
  85           GFC_DIMENSION_SET(retarray->dim[1], 0,
  86                             GFC_DESCRIPTOR_EXTENT(b,1) - 1,
  87                             GFC_DESCRIPTOR_EXTENT(retarray,0));
  88         }
  89
  90       retarray->base_addr
  91         = xmallocarray (size0 ((array_t *) retarray), sizeof ('rtype_name`));
  92       retarray->offset = 0;
  93     }
  94     else if (unlikely (compile_options.bounds_check))
  95       {
  96         index_type ret_extent, arg_extent;
  97
  98         if (GFC_DESCRIPTOR_RANK (a) == 1)
  99           {
 100             arg_extent = GFC_DESCRIPTOR_EXTENT(b,1);
 101             ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
 102             if (arg_extent != ret_extent)
 103               runtime_error ("Incorrect extent in return array in"
 104                              " MATMUL intrinsic: is %ld, should be %ld",
 105                              (long int) ret_extent, (long int) arg_extent);
 106           }
 107         else if (GFC_DESCRIPTOR_RANK (b) == 1)
 108           {
 109             arg_extent = GFC_DESCRIPTOR_EXTENT(a,0);
 110             ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
 111             if (arg_extent != ret_extent)
 112               runtime_error ("Incorrect extent in return array in"
 113                              " MATMUL intrinsic: is %ld, should be %ld",
 114                              (long int) ret_extent, (long int) arg_extent);
 115           }
 116         else
 117           {
 118             arg_extent = GFC_DESCRIPTOR_EXTENT(a,0);
 119             ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
 120             if (arg_extent != ret_extent)
 121               runtime_error ("Incorrect extent in return array in"
 122                              " MATMUL intrinsic for dimension 1:"
 123                              " is %ld, should be %ld",
 124                              (long int) ret_extent, (long int) arg_extent);
 125
 126             arg_extent = GFC_DESCRIPTOR_EXTENT(b,1);
 127             ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,1);
 128             if (arg_extent != ret_extent)
 129               runtime_error ("Incorrect extent in return array in"
 130                              " MATMUL intrinsic for dimension 2:"
 131                              " is %ld, should be %ld",
 132                              (long int) ret_extent, (long int) arg_extent);
 133           }
 134       }
 135
 136   abase = a->base_addr;
 137   a_kind = GFC_DESCRIPTOR_SIZE (a);
 138
 139   if (a_kind == 1 || a_kind == 2 || a_kind == 4 || a_kind == 8
 140 #ifdef HAVE_GFC_LOGICAL_16
 141      || a_kind == 16
 142 #endif
 143      )
 144     abase = GFOR_POINTER_TO_L1 (abase, a_kind);
 145   else
 146     internal_error (NULL, "Funny sized logical array");
 147
 148   bbase = b->base_addr;
 149   b_kind = GFC_DESCRIPTOR_SIZE (b);
 150
 151   if (b_kind == 1 || b_kind == 2 || b_kind == 4 || b_kind == 8
 152 #ifdef HAVE_GFC_LOGICAL_16
 153      || b_kind == 16
 154 #endif
 155      )
 156     bbase = GFOR_POINTER_TO_L1 (bbase, b_kind);
 157   else
 158     internal_error (NULL, "Funny sized logical array");
 159
 160   dest = retarray->base_addr;
 161 '
 162 sinclude(`matmul_asm_'rtype_code`.m4')dnl
 163 `
 164   if (GFC_DESCRIPTOR_RANK (retarray) == 1)
 165     {
 166       rxstride = GFC_DESCRIPTOR_STRIDE(retarray,0);
 167       rystride = rxstride;
 168     }
 169   else
 170     {
 171       rxstride = GFC_DESCRIPTOR_STRIDE(retarray,0);
 172       rystride = GFC_DESCRIPTOR_STRIDE(retarray,1);
 173     }
 174
 175   /* If we have rank 1 parameters, zero the absent stride, and set the size to
 176      one.  */
 177   if (GFC_DESCRIPTOR_RANK (a) == 1)
 178     {
 179       astride = GFC_DESCRIPTOR_STRIDE_BYTES(a,0);
 180       count = GFC_DESCRIPTOR_EXTENT(a,0);
 181       xstride = 0;
 182       rxstride = 0;
 183       xcount = 1;
 184     }
 185   else
 186     {
 187       astride = GFC_DESCRIPTOR_STRIDE_BYTES(a,1);
 188       count = GFC_DESCRIPTOR_EXTENT(a,1);
 189       xstride = GFC_DESCRIPTOR_STRIDE_BYTES(a,0);
 190       xcount = GFC_DESCRIPTOR_EXTENT(a,0);
 191     }
 192   if (GFC_DESCRIPTOR_RANK (b) == 1)
 193     {
 194       bstride = GFC_DESCRIPTOR_STRIDE_BYTES(b,0);
 195       assert(count == GFC_DESCRIPTOR_EXTENT(b,0));
 196       ystride = 0;
 197       rystride = 0;
 198       ycount = 1;
 199     }
 200   else
 201     {
 202       bstride = GFC_DESCRIPTOR_STRIDE_BYTES(b,0);
 203       assert(count == GFC_DESCRIPTOR_EXTENT(b,0));
 204       ystride = GFC_DESCRIPTOR_STRIDE_BYTES(b,1);
 205       ycount = GFC_DESCRIPTOR_EXTENT(b,1);
 206     }
 207
 208   for (y = 0; y < ycount; y++)
 209     {
 210       for (x = 0; x < xcount; x++)
 211         {
 212           /* Do the summation for this element.  For real and integer types
 213              this is the same as DOT_PRODUCT.  For complex types we use do
 214              a*b, not conjg(a)*b.  */
 215           pa = abase;
 216           pb = bbase;
 217           *dest = 0;
 218
 219           for (n = 0; n < count; n++)
 220             {
 221               if (*pa && *pb)
 222                 {
 223                   *dest = 1;
 224                   break;
 225                 }
 226               pa += astride;
 227               pb += bstride;
 228             }
 229
 230           dest += rxstride;
 231           abase += xstride;
 232         }
 233       abase -= xstride * xcount;
 234       bbase += ystride;
 235       dest += rystride - (rxstride * xcount);
 236     }
 237 }
 238
 239 #endif
 240 '