libgfortran/m4/matmull.m4

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