libgfortran/m4/matmull.m4

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