libgfortran/m4/matmull.m4

   1 `/* Implementation of the MATMUL intrinsic
   2    Copyright 2002, 2005, 2006, 2007, 2009 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 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 <stdlib.h>
  28 #include <assert.h>'
  29
  30 include(iparm.m4)dnl
  31
  32 `#if defined (HAVE_'rtype_name`)
  33
  34 /* Dimensions: retarray(x,y) a(x, count) b(count,y).
  35    Either a or b can be rank 1.  In this case x or y is 1.  */
  36
  37 extern void matmul_'rtype_code` ('rtype` * const restrict,
  38         gfc_array_l1 * const restrict, gfc_array_l1 * const restrict);
  39 export_proto(matmul_'rtype_code`);
  40
  41 void
  42 matmul_'rtype_code` ('rtype` * const restrict retarray,
  43         gfc_array_l1 * const restrict a, gfc_array_l1 * const restrict b)
  44 {
  45   const GFC_LOGICAL_1 * restrict abase;
  46   const GFC_LOGICAL_1 * restrict bbase;
  47   'rtype_name` * restrict dest;
  48   index_type rxstride;
  49   index_type rystride;
  50   index_type xcount;
  51   index_type ycount;
  52   index_type xstride;
  53   index_type ystride;
  54   index_type x;
  55   index_type y;
  56   int a_kind;
  57   int b_kind;
  58
  59   const GFC_LOGICAL_1 * restrict pa;
  60   const GFC_LOGICAL_1 * restrict pb;
  61   index_type astride;
  62   index_type bstride;
  63   index_type count;
  64   index_type n;
  65
  66   assert (GFC_DESCRIPTOR_RANK (a) == 2
  67           || GFC_DESCRIPTOR_RANK (b) == 2);
  68
  69   if (retarray->data == NULL)
  70     {
  71       if (GFC_DESCRIPTOR_RANK (a) == 1)
  72         {
  73           retarray->dim[0].lbound = 0;
  74           retarray->dim[0].ubound = b->dim[1].ubound - b->dim[1].lbound;
  75           retarray->dim[0].stride = 1;
  76         }
  77       else if (GFC_DESCRIPTOR_RANK (b) == 1)
  78         {
  79           retarray->dim[0].lbound = 0;
  80           retarray->dim[0].ubound = a->dim[0].ubound - a->dim[0].lbound;
  81           retarray->dim[0].stride = 1;
  82         }
  83       else
  84         {
  85           retarray->dim[0].lbound = 0;
  86           retarray->dim[0].ubound = a->dim[0].ubound - a->dim[0].lbound;
  87           retarray->dim[0].stride = 1;
  88
  89           retarray->dim[1].lbound = 0;
  90           retarray->dim[1].ubound = b->dim[1].ubound - b->dim[1].lbound;
  91           retarray->dim[1].stride = retarray->dim[0].ubound+1;
  92         }
  93
  94       retarray->data
  95         = internal_malloc_size (sizeof ('rtype_name`) * size0 ((array_t *) retarray));
  96       retarray->offset = 0;
  97     }
  98     else if (unlikely (compile_options.bounds_check))
  99       {
 100         index_type ret_extent, arg_extent;
 101
 102         if (GFC_DESCRIPTOR_RANK (a) == 1)
 103           {
 104             arg_extent = b->dim[1].ubound + 1 - b->dim[1].lbound;
 105             ret_extent = retarray->dim[0].ubound + 1 - retarray->dim[0].lbound;
 106             if (arg_extent != ret_extent)
 107               runtime_error ("Incorrect extent in return array in"
 108                              " MATMUL intrinsic: is %ld, should be %ld",
 109                              (long int) ret_extent, (long int) arg_extent);
 110           }
 111         else if (GFC_DESCRIPTOR_RANK (b) == 1)
 112           {
 113             arg_extent = a->dim[0].ubound + 1 - a->dim[0].lbound;
 114             ret_extent = retarray->dim[0].ubound + 1 - retarray->dim[0].lbound;
 115             if (arg_extent != ret_extent)
 116               runtime_error ("Incorrect extent in return array in"
 117                              " MATMUL intrinsic: is %ld, should be %ld",
 118                              (long int) ret_extent, (long int) arg_extent);
 119           }
 120         else
 121           {
 122             arg_extent = a->dim[0].ubound + 1 - a->dim[0].lbound;
 123             ret_extent = retarray->dim[0].ubound + 1 - retarray->dim[0].lbound;
 124             if (arg_extent != ret_extent)
 125               runtime_error ("Incorrect extent in return array in"
 126                              " MATMUL intrinsic for dimension 1:"
 127                              " is %ld, should be %ld",
 128                              (long int) ret_extent, (long int) arg_extent);
 129
 130             arg_extent = b->dim[1].ubound + 1 - b->dim[1].lbound;
 131             ret_extent = retarray->dim[1].ubound + 1 - retarray->dim[1].lbound;
 132             if (arg_extent != ret_extent)
 133               runtime_error ("Incorrect extent in return array in"
 134                              " MATMUL intrinsic for dimension 2:"
 135                              " is %ld, should be %ld",
 136                              (long int) ret_extent, (long int) arg_extent);
 137           }
 138       }
 139
 140   abase = a->data;
 141   a_kind = GFC_DESCRIPTOR_SIZE (a);
 142
 143   if (a_kind == 1 || a_kind == 2 || a_kind == 4 || a_kind == 8
 144 #ifdef HAVE_GFC_LOGICAL_16
 145      || a_kind == 16
 146 #endif
 147      )
 148     abase = GFOR_POINTER_TO_L1 (abase, a_kind);
 149   else
 150     internal_error (NULL, "Funny sized logical array");
 151
 152   bbase = b->data;
 153   b_kind = GFC_DESCRIPTOR_SIZE (b);
 154
 155   if (b_kind == 1 || b_kind == 2 || b_kind == 4 || b_kind == 8
 156 #ifdef HAVE_GFC_LOGICAL_16
 157      || b_kind == 16
 158 #endif
 159      )
 160     bbase = GFOR_POINTER_TO_L1 (bbase, b_kind);
 161   else
 162     internal_error (NULL, "Funny sized logical array");
 163
 164   dest = retarray->data;
 165 '
 166 sinclude(`matmul_asm_'rtype_code`.m4')dnl
 167 `
 168   if (GFC_DESCRIPTOR_RANK (retarray) == 1)
 169     {
 170       rxstride = retarray->dim[0].stride;
 171       rystride = rxstride;
 172     }
 173   else
 174     {
 175       rxstride = retarray->dim[0].stride;
 176       rystride = retarray->dim[1].stride;
 177     }
 178
 179   /* If we have rank 1 parameters, zero the absent stride, and set the size to
 180      one.  */
 181   if (GFC_DESCRIPTOR_RANK (a) == 1)
 182     {
 183       astride = a->dim[0].stride * a_kind;
 184       count = a->dim[0].ubound + 1 - a->dim[0].lbound;
 185       xstride = 0;
 186       rxstride = 0;
 187       xcount = 1;
 188     }
 189   else
 190     {
 191       astride = a->dim[1].stride * a_kind;
 192       count = a->dim[1].ubound + 1 - a->dim[1].lbound;
 193       xstride = a->dim[0].stride * a_kind;
 194       xcount = a->dim[0].ubound + 1 - a->dim[0].lbound;
 195     }
 196   if (GFC_DESCRIPTOR_RANK (b) == 1)
 197     {
 198       bstride = b->dim[0].stride * b_kind;
 199       assert(count == b->dim[0].ubound + 1 - b->dim[0].lbound);
 200       ystride = 0;
 201       rystride = 0;
 202       ycount = 1;
 203     }
 204   else
 205     {
 206       bstride = b->dim[0].stride * b_kind;
 207       assert(count == b->dim[0].ubound + 1 - b->dim[0].lbound);
 208       ystride = b->dim[1].stride * b_kind;
 209       ycount = b->dim[1].ubound + 1 - b->dim[1].lbound;
 210     }
 211
 212   for (y = 0; y < ycount; y++)
 213     {
 214       for (x = 0; x < xcount; x++)
 215         {
 216           /* Do the summation for this element.  For real and integer types
 217              this is the same as DOT_PRODUCT.  For complex types we use do
 218              a*b, not conjg(a)*b.  */
 219           pa = abase;
 220           pb = bbase;
 221           *dest = 0;
 222
 223           for (n = 0; n < count; n++)
 224             {
 225               if (*pa && *pb)
 226                 {
 227                   *dest = 1;
 228                   break;
 229                 }
 230               pa += astride;
 231               pb += bstride;
 232             }
 233
 234           dest += rxstride;
 235           abase += xstride;
 236         }
 237       abase -= xstride * xcount;
 238       bbase += ystride;
 239       dest += rystride - (rxstride * xcount);
 240     }
 241 }
 242
 243 #endif
 244 '