libgfortran/generated/matmul_l16.c

   1 /* Implementation of the MATMUL intrinsic
   2    Copyright 2002, 2005 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 "config.h"
  32 #include <stdlib.h>
  33 #include <assert.h>
  34 #include "libgfortran.h"
  35
  36 #if defined (HAVE_GFC_LOGICAL_16)
  37
  38 /* Dimensions: retarray(x,y) a(x, count) b(count,y).
  39    Either a or b can be rank 1.  In this case x or y is 1.  */
  40
  41 extern void matmul_l16 (gfc_array_l16 * const restrict,
  42         gfc_array_l4 * const restrict, gfc_array_l4 * const restrict);
  43 export_proto(matmul_l16);
  44
  45 void
  46 matmul_l16 (gfc_array_l16 * const restrict retarray,
  47         gfc_array_l4 * const restrict a, gfc_array_l4 * const restrict b)
  48 {
  49   const GFC_INTEGER_4 * restrict abase;
  50   const GFC_INTEGER_4 * restrict bbase;
  51   GFC_LOGICAL_16 * restrict dest;
  52   index_type rxstride;
  53   index_type rystride;
  54   index_type xcount;
  55   index_type ycount;
  56   index_type xstride;
  57   index_type ystride;
  58   index_type x;
  59   index_type y;
  60
  61   const GFC_INTEGER_4 * restrict pa;
  62   const GFC_INTEGER_4 * restrict pb;
  63   index_type astride;
  64   index_type bstride;
  65   index_type count;
  66   index_type n;
  67
  68   assert (GFC_DESCRIPTOR_RANK (a) == 2
  69           || GFC_DESCRIPTOR_RANK (b) == 2);
  70
  71   if (retarray->data == NULL)
  72     {
  73       if (GFC_DESCRIPTOR_RANK (a) == 1)
  74         {
  75           retarray->dim[0].lbound = 0;
  76           retarray->dim[0].ubound = b->dim[1].ubound - b->dim[1].lbound;
  77           retarray->dim[0].stride = 1;
  78         }
  79       else if (GFC_DESCRIPTOR_RANK (b) == 1)
  80         {
  81           retarray->dim[0].lbound = 0;
  82           retarray->dim[0].ubound = a->dim[0].ubound - a->dim[0].lbound;
  83           retarray->dim[0].stride = 1;
  84         }
  85       else
  86         {
  87           retarray->dim[0].lbound = 0;
  88           retarray->dim[0].ubound = a->dim[0].ubound - a->dim[0].lbound;
  89           retarray->dim[0].stride = 1;
  90
  91           retarray->dim[1].lbound = 0;
  92           retarray->dim[1].ubound = b->dim[1].ubound - b->dim[1].lbound;
  93           retarray->dim[1].stride = retarray->dim[0].ubound+1;
  94         }
  95
  96       retarray->data
  97         = internal_malloc_size (sizeof (GFC_LOGICAL_16) * size0 ((array_t *) retarray));
  98       retarray->offset = 0;
  99     }
 100
 101   abase = a->data;
 102   if (GFC_DESCRIPTOR_SIZE (a) != 4)
 103     {
 104       assert (GFC_DESCRIPTOR_SIZE (a) == 8);
 105       abase = GFOR_POINTER_L8_TO_L4 (abase);
 106     }
 107   bbase = b->data;
 108   if (GFC_DESCRIPTOR_SIZE (b) != 4)
 109     {
 110       assert (GFC_DESCRIPTOR_SIZE (b) == 8);
 111       bbase = GFOR_POINTER_L8_TO_L4 (bbase);
 112     }
 113   dest = retarray->data;
 114
 115   if (retarray->dim[0].stride == 0)
 116     retarray->dim[0].stride = 1;
 117   if (a->dim[0].stride == 0)
 118     a->dim[0].stride = 1;
 119   if (b->dim[0].stride == 0)
 120     b->dim[0].stride = 1;
 121
 122
 123   if (GFC_DESCRIPTOR_RANK (retarray) == 1)
 124     {
 125       rxstride = retarray->dim[0].stride;
 126       rystride = rxstride;
 127     }
 128   else
 129     {
 130       rxstride = retarray->dim[0].stride;
 131       rystride = retarray->dim[1].stride;
 132     }
 133
 134   /* If we have rank 1 parameters, zero the absent stride, and set the size to
 135      one.  */
 136   if (GFC_DESCRIPTOR_RANK (a) == 1)
 137     {
 138       astride = a->dim[0].stride;
 139       count = a->dim[0].ubound + 1 - a->dim[0].lbound;
 140       xstride = 0;
 141       rxstride = 0;
 142       xcount = 1;
 143     }
 144   else
 145     {
 146       astride = a->dim[1].stride;
 147       count = a->dim[1].ubound + 1 - a->dim[1].lbound;
 148       xstride = a->dim[0].stride;
 149       xcount = a->dim[0].ubound + 1 - a->dim[0].lbound;
 150     }
 151   if (GFC_DESCRIPTOR_RANK (b) == 1)
 152     {
 153       bstride = b->dim[0].stride;
 154       assert(count == b->dim[0].ubound + 1 - b->dim[0].lbound);
 155       ystride = 0;
 156       rystride = 0;
 157       ycount = 1;
 158     }
 159   else
 160     {
 161       bstride = b->dim[0].stride;
 162       assert(count == b->dim[0].ubound + 1 - b->dim[0].lbound);
 163       ystride = b->dim[1].stride;
 164       ycount = b->dim[1].ubound + 1 - b->dim[1].lbound;
 165     }
 166
 167   for (y = 0; y < ycount; y++)
 168     {
 169       for (x = 0; x < xcount; x++)
 170         {
 171           /* Do the summation for this element.  For real and integer types
 172              this is the same as DOT_PRODUCT.  For complex types we use do
 173              a*b, not conjg(a)*b.  */
 174           pa = abase;
 175           pb = bbase;
 176           *dest = 0;
 177
 178           for (n = 0; n < count; n++)
 179             {
 180               if (*pa && *pb)
 181                 {
 182                   *dest = 1;
 183                   break;
 184                 }
 185               pa += astride;
 186               pb += bstride;
 187             }
 188
 189           dest += rxstride;
 190           abase += xstride;
 191         }
 192       abase -= xstride * xcount;
 193       bbase += ystride;
 194       dest += rystride - (rxstride * xcount);
 195     }
 196 }
 197
 198 #endif