libgfortran/m4/matmull.m4

   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 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., 59 Temple Place - Suite 330,
  29 Boston, MA 02111-1307, USA.  */
  30
  31 #include "config.h"
  32 #include <stdlib.h>
  33 #include <assert.h>
  34 #include "libgfortran.h"'
  35 include(iparm.m4)dnl
  36
  37 /* Dimensions: retarray(x,y) a(x, count) b(count,y).
  38    Either a or b can be rank 1.  In this case x or y is 1.  */
  39
  40 extern void matmul_`'rtype_code (rtype *, gfc_array_l4 *, gfc_array_l4 *);
  41 export_proto(matmul_`'rtype_code);
  42
  43 void
  44 matmul_`'rtype_code (rtype * retarray, gfc_array_l4 * a, gfc_array_l4 * b)
  45 {
  46   GFC_INTEGER_4 *abase;
  47   GFC_INTEGER_4 *bbase;
  48   rtype_name *dest;
  49   index_type rxstride;
  50   index_type rystride;
  51   index_type xcount;
  52   index_type ycount;
  53   index_type xstride;
  54   index_type ystride;
  55   index_type x;
  56   index_type y;
  57
  58   GFC_INTEGER_4 *pa;
  59   GFC_INTEGER_4 *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 (rtype_name) * size0 (retarray));
  95       retarray->base = 0;
  96     }
  97
  98   abase = a->data;
  99   if (GFC_DESCRIPTOR_SIZE (a) != 4)
 100     {
 101       assert (GFC_DESCRIPTOR_SIZE (a) == 8);
 102       abase = GFOR_POINTER_L8_TO_L4 (abase);
 103       astride <<= 1;
 104     }
 105   bbase = b->data;
 106   if (GFC_DESCRIPTOR_SIZE (b) != 4)
 107     {
 108       assert (GFC_DESCRIPTOR_SIZE (b) == 8);
 109       bbase = GFOR_POINTER_L8_TO_L4 (bbase);
 110       bstride <<= 1;
 111     }
 112   dest = retarray->data;
 113
 114   if (retarray->dim[0].stride == 0)
 115     retarray->dim[0].stride = 1;
 116   if (a->dim[0].stride == 0)
 117     a->dim[0].stride = 1;
 118   if (b->dim[0].stride == 0)
 119     b->dim[0].stride = 1;
 120
 121 sinclude(`matmul_asm_'rtype_code`.m4')dnl
 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 }