libgfortran/generated/matmul_i8.c

   1 /* Implementation of the MATMUL intrinsic
   2    Copyright 2002 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 Lesser General Public
   9 License as published by the Free Software Foundation; either
  10 version 2.1 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 Lesser General Public License for more details.
  16
  17 You should have received a copy of the GNU Lesser General Public
  18 License along with libgfor; see the file COPYING.LIB.  If not,
  19 write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  20 Boston, MA 02111-1307, USA.  */
  21
  22 #include "config.h"
  23 #include <stdlib.h>
  24 #include <string.h>
  25 #include <assert.h>
  26 #include "libgfortran.h"
  27
  28 /* This is a C version of the following fortran pseudo-code. The key
  29    point is the loop order -- we access all arrays column-first, which
  30    improves the performance enough to boost galgel spec score by 50%.
  31
  32    DIMENSION A(M,COUNT), B(COUNT,N), C(M,N)
  33    C = 0
  34    DO J=1,N
  35      DO K=1,COUNT
  36        DO I=1,M
  37          C(I,J) = C(I,J)+A(I,K)*B(K,J)
  38 */
  39
  40 void
  41 __matmul_i8 (gfc_array_i8 * retarray, gfc_array_i8 * a, gfc_array_i8 * b)
  42 {
  43   GFC_INTEGER_8 *abase;
  44   GFC_INTEGER_8 *bbase;
  45   GFC_INTEGER_8 *dest;
  46
  47   index_type rxstride, rystride, axstride, aystride, bxstride, bystride;
  48   index_type x, y, n, count, xcount, ycount;
  49
  50   assert (GFC_DESCRIPTOR_RANK (a) == 2
  51           || GFC_DESCRIPTOR_RANK (b) == 2);
  52
  53 /* C[xcount,ycount] = A[xcount, count] * B[count,ycount]
  54
  55    Either A or B (but not both) can be rank 1:
  56
  57    o One-dimensional argument A is implicitly treated as a row matrix
  58      dimensioned [1,count], so xcount=1.
  59
  60    o One-dimensional argument B is implicitly treated as a column matrix
  61      dimensioned [count, 1], so ycount=1.
  62   */
  63
  64   if (retarray->data == NULL)
  65     {
  66       if (GFC_DESCRIPTOR_RANK (a) == 1)
  67         {
  68           retarray->dim[0].lbound = 0;
  69           retarray->dim[0].ubound = b->dim[1].ubound - b->dim[1].lbound;
  70           retarray->dim[0].stride = 1;
  71         }
  72       else if (GFC_DESCRIPTOR_RANK (b) == 1)
  73         {
  74           retarray->dim[0].lbound = 0;
  75           retarray->dim[0].ubound = a->dim[0].ubound - a->dim[0].lbound;
  76           retarray->dim[0].stride = 1;
  77         }
  78       else
  79         {
  80           retarray->dim[0].lbound = 0;
  81           retarray->dim[0].ubound = a->dim[0].ubound - a->dim[0].lbound;
  82           retarray->dim[0].stride = 1;
  83
  84           retarray->dim[1].lbound = 0;
  85           retarray->dim[1].ubound = b->dim[1].ubound - b->dim[1].lbound;
  86           retarray->dim[1].stride = retarray->dim[0].ubound+1;
  87         }
  88
  89       retarray->data
  90         = internal_malloc_size (sizeof (GFC_INTEGER_8) * size0 (retarray));
  91       retarray->base = 0;
  92     }
  93
  94   abase = a->data;
  95   bbase = b->data;
  96   dest = retarray->data;
  97
  98   if (retarray->dim[0].stride == 0)
  99     retarray->dim[0].stride = 1;
 100   if (a->dim[0].stride == 0)
 101     a->dim[0].stride = 1;
 102   if (b->dim[0].stride == 0)
 103     b->dim[0].stride = 1;
 104
 105
 106   if (GFC_DESCRIPTOR_RANK (retarray) == 1)
 107     {
 108       /* One-dimensional result may be addressed in the code below
 109          either as a row or a column matrix. We want both cases to
 110          work. */
 111       rxstride = rystride = retarray->dim[0].stride;
 112     }
 113   else
 114     {
 115       rxstride = retarray->dim[0].stride;
 116       rystride = retarray->dim[1].stride;
 117     }
 118
 119
 120   if (GFC_DESCRIPTOR_RANK (a) == 1)
 121     {
 122       /* Treat it as a a row matrix A[1,count]. */
 123       axstride = a->dim[0].stride;
 124       aystride = 1;
 125
 126       xcount = 1;
 127       count = a->dim[0].ubound + 1 - a->dim[0].lbound;
 128     }
 129   else
 130     {
 131       axstride = a->dim[0].stride;
 132       aystride = a->dim[1].stride;
 133
 134       count = a->dim[1].ubound + 1 - a->dim[1].lbound;
 135       xcount = a->dim[0].ubound + 1 - a->dim[0].lbound;
 136     }
 137
 138   assert(count == b->dim[0].ubound + 1 - b->dim[0].lbound);
 139
 140   if (GFC_DESCRIPTOR_RANK (b) == 1)
 141     {
 142       /* Treat it as a column matrix B[count,1] */
 143       bxstride = b->dim[0].stride;
 144
 145       /* bystride should never be used for 1-dimensional b.
 146          in case it is we want it to cause a segfault, rather than
 147          an incorrect result. */
 148       bystride = 0xDEADBEEF;
 149       ycount = 1;
 150     }
 151   else
 152     {
 153       bxstride = b->dim[0].stride;
 154       bystride = b->dim[1].stride;
 155       ycount = b->dim[1].ubound + 1 - b->dim[1].lbound;
 156     }
 157
 158   assert (a->base == 0);
 159   assert (b->base == 0);
 160   assert (retarray->base == 0);
 161
 162   abase = a->data;
 163   bbase = b->data;
 164   dest = retarray->data;
 165
 166   if (rxstride == 1 && axstride == 1 && bxstride == 1)
 167     {
 168       GFC_INTEGER_8 *bbase_y;
 169       GFC_INTEGER_8 *dest_y;
 170       GFC_INTEGER_8 *abase_n;
 171       GFC_INTEGER_8 bbase_yn;
 172
 173       memset (dest, 0, (sizeof (GFC_INTEGER_8) * size0(retarray)));
 174
 175       for (y = 0; y < ycount; y++)
 176         {
 177           bbase_y = bbase + y*bystride;
 178           dest_y = dest + y*rystride;
 179           for (n = 0; n < count; n++)
 180             {
 181               abase_n = abase + n*aystride;
 182               bbase_yn = bbase_y[n];
 183               for (x = 0; x < xcount; x++)
 184                 {
 185                   dest_y[x] += abase_n[x] * bbase_yn;
 186                 }
 187             }
 188         }
 189     }
 190   else
 191     {
 192       for (y = 0; y < ycount; y++)
 193         for (x = 0; x < xcount; x++)
 194           dest[x*rxstride + y*rystride] = (GFC_INTEGER_8)0;
 195
 196       for (y = 0; y < ycount; y++)
 197         for (n = 0; n < count; n++)
 198           for (x = 0; x < xcount; x++)
 199             /* dest[x,y] += a[x,n] * b[n,y] */
 200             dest[x*rxstride + y*rystride] += abase[x*axstride + n*aystride] * bbase[n*bxstride + y*bystride];
 201     }
 202 }