libgfortran/generated/matmul_l16.c

   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
  31 #if defined (HAVE_GFC_LOGICAL_16)
  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_l16 (gfc_array_l16 * const restrict,
  37         gfc_array_l1 * const restrict, gfc_array_l1 * const restrict);
  38 export_proto(matmul_l16);
  39
  40 void
  41 matmul_l16 (gfc_array_l16 * 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   GFC_LOGICAL_16 * 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->data == NULL)
  69     {
  70       if (GFC_DESCRIPTOR_RANK (a) == 1)
  71         {
  72           retarray->dim[0].lbound = 0;
  73           retarray->dim[0].ubound = b->dim[1].ubound - b->dim[1].lbound;
  74           retarray->dim[0].stride = 1;
  75         }
  76       else if (GFC_DESCRIPTOR_RANK (b) == 1)
  77         {
  78           retarray->dim[0].lbound = 0;
  79           retarray->dim[0].ubound = a->dim[0].ubound - a->dim[0].lbound;
  80           retarray->dim[0].stride = 1;
  81         }
  82       else
  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           retarray->dim[1].lbound = 0;
  89           retarray->dim[1].ubound = b->dim[1].ubound - b->dim[1].lbound;
  90           retarray->dim[1].stride = retarray->dim[0].ubound+1;
  91         }
  92
  93       retarray->data
  94         = internal_malloc_size (sizeof (GFC_LOGICAL_16) * size0 ((array_t *) retarray));
  95       retarray->offset = 0;
  96     }
  97     else if (unlikely (compile_options.bounds_check))
  98       {
  99         index_type ret_extent, arg_extent;
 100
 101         if (GFC_DESCRIPTOR_RANK (a) == 1)
 102           {
 103             arg_extent = b->dim[1].ubound + 1 - b->dim[1].lbound;
 104             ret_extent = retarray->dim[0].ubound + 1 - retarray->dim[0].lbound;
 105             if (arg_extent != ret_extent)
 106               runtime_error ("Incorrect extent in return array in"
 107                              " MATMUL intrinsic: is %ld, should be %ld",
 108                              (long int) ret_extent, (long int) arg_extent);
 109           }
 110         else if (GFC_DESCRIPTOR_RANK (b) == 1)
 111           {
 112             arg_extent = a->dim[0].ubound + 1 - a->dim[0].lbound;
 113             ret_extent = retarray->dim[0].ubound + 1 - retarray->dim[0].lbound;
 114             if (arg_extent != ret_extent)
 115               runtime_error ("Incorrect extent in return array in"
 116                              " MATMUL intrinsic: is %ld, should be %ld",
 117                              (long int) ret_extent, (long int) arg_extent);
 118           }
 119         else
 120           {
 121             arg_extent = a->dim[0].ubound + 1 - a->dim[0].lbound;
 122             ret_extent = retarray->dim[0].ubound + 1 - retarray->dim[0].lbound;
 123             if (arg_extent != ret_extent)
 124               runtime_error ("Incorrect extent in return array in"
 125                              " MATMUL intrinsic for dimension 1:"
 126                              " is %ld, should be %ld",
 127                              (long int) ret_extent, (long int) arg_extent);
 128
 129             arg_extent = b->dim[1].ubound + 1 - b->dim[1].lbound;
 130             ret_extent = retarray->dim[1].ubound + 1 - retarray->dim[1].lbound;
 131             if (arg_extent != ret_extent)
 132               runtime_error ("Incorrect extent in return array in"
 133                              " MATMUL intrinsic for dimension 2:"
 134                              " is %ld, should be %ld",
 135                              (long int) ret_extent, (long int) arg_extent);
 136           }
 137       }
 138
 139   abase = a->data;
 140   a_kind = GFC_DESCRIPTOR_SIZE (a);
 141
 142   if (a_kind == 1 || a_kind == 2 || a_kind == 4 || a_kind == 8
 143 #ifdef HAVE_GFC_LOGICAL_16
 144      || a_kind == 16
 145 #endif
 146      )
 147     abase = GFOR_POINTER_TO_L1 (abase, a_kind);
 148   else
 149     internal_error (NULL, "Funny sized logical array");
 150
 151   bbase = b->data;
 152   b_kind = GFC_DESCRIPTOR_SIZE (b);
 153
 154   if (b_kind == 1 || b_kind == 2 || b_kind == 4 || b_kind == 8
 155 #ifdef HAVE_GFC_LOGICAL_16
 156      || b_kind == 16
 157 #endif
 158      )
 159     bbase = GFOR_POINTER_TO_L1 (bbase, b_kind);
 160   else
 161     internal_error (NULL, "Funny sized logical array");
 162
 163   dest = retarray->data;
 164
 165
 166   if (GFC_DESCRIPTOR_RANK (retarray) == 1)
 167     {
 168       rxstride = retarray->dim[0].stride;
 169       rystride = rxstride;
 170     }
 171   else
 172     {
 173       rxstride = retarray->dim[0].stride;
 174       rystride = retarray->dim[1].stride;
 175     }
 176
 177   /* If we have rank 1 parameters, zero the absent stride, and set the size to
 178      one.  */
 179   if (GFC_DESCRIPTOR_RANK (a) == 1)
 180     {
 181       astride = a->dim[0].stride * a_kind;
 182       count = a->dim[0].ubound + 1 - a->dim[0].lbound;
 183       xstride = 0;
 184       rxstride = 0;
 185       xcount = 1;
 186     }
 187   else
 188     {
 189       astride = a->dim[1].stride * a_kind;
 190       count = a->dim[1].ubound + 1 - a->dim[1].lbound;
 191       xstride = a->dim[0].stride * a_kind;
 192       xcount = a->dim[0].ubound + 1 - a->dim[0].lbound;
 193     }
 194   if (GFC_DESCRIPTOR_RANK (b) == 1)
 195     {
 196       bstride = b->dim[0].stride * b_kind;
 197       assert(count == b->dim[0].ubound + 1 - b->dim[0].lbound);
 198       ystride = 0;
 199       rystride = 0;
 200       ycount = 1;
 201     }
 202   else
 203     {
 204       bstride = b->dim[0].stride * b_kind;
 205       assert(count == b->dim[0].ubound + 1 - b->dim[0].lbound);
 206       ystride = b->dim[1].stride * b_kind;
 207       ycount = b->dim[1].ubound + 1 - b->dim[1].lbound;
 208     }
 209
 210   for (y = 0; y < ycount; y++)
 211     {
 212       for (x = 0; x < xcount; x++)
 213         {
 214           /* Do the summation for this element.  For real and integer types
 215              this is the same as DOT_PRODUCT.  For complex types we use do
 216              a*b, not conjg(a)*b.  */
 217           pa = abase;
 218           pb = bbase;
 219           *dest = 0;
 220
 221           for (n = 0; n < count; n++)
 222             {
 223               if (*pa && *pb)
 224                 {
 225                   *dest = 1;
 226                   break;
 227                 }
 228               pa += astride;
 229               pb += bstride;
 230             }
 231
 232           dest += rxstride;
 233           abase += xstride;
 234         }
 235       abase -= xstride * xcount;
 236       bbase += ystride;
 237       dest += rystride - (rxstride * xcount);
 238     }
 239 }
 240
 241 #endif
 242