libgfortran/m4/matmull.m4

   1 `/* Implementation of the MATMUL intrinsic
   2    Copyright 2002, 2005, 2006 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 sinclude(`matmul_asm_'rtype_code`.m4')dnl
 117 `
 118   if (GFC_DESCRIPTOR_RANK (retarray) == 1)
 119     {
 120       rxstride = retarray->dim[0].stride;
 121       rystride = rxstride;
 122     }
 123   else
 124     {
 125       rxstride = retarray->dim[0].stride;
 126       rystride = retarray->dim[1].stride;
 127     }
 128
 129   /* If we have rank 1 parameters, zero the absent stride, and set the size to
 130      one.  */
 131   if (GFC_DESCRIPTOR_RANK (a) == 1)
 132     {
 133       astride = a->dim[0].stride;
 134       count = a->dim[0].ubound + 1 - a->dim[0].lbound;
 135       xstride = 0;
 136       rxstride = 0;
 137       xcount = 1;
 138     }
 139   else
 140     {
 141       astride = a->dim[1].stride;
 142       count = a->dim[1].ubound + 1 - a->dim[1].lbound;
 143       xstride = a->dim[0].stride;
 144       xcount = a->dim[0].ubound + 1 - a->dim[0].lbound;
 145     }
 146   if (GFC_DESCRIPTOR_RANK (b) == 1)
 147     {
 148       bstride = b->dim[0].stride;
 149       assert(count == b->dim[0].ubound + 1 - b->dim[0].lbound);
 150       ystride = 0;
 151       rystride = 0;
 152       ycount = 1;
 153     }
 154   else
 155     {
 156       bstride = b->dim[0].stride;
 157       assert(count == b->dim[0].ubound + 1 - b->dim[0].lbound);
 158       ystride = b->dim[1].stride;
 159       ycount = b->dim[1].ubound + 1 - b->dim[1].lbound;
 160     }
 161
 162   for (y = 0; y < ycount; y++)
 163     {
 164       for (x = 0; x < xcount; x++)
 165         {
 166           /* Do the summation for this element.  For real and integer types
 167              this is the same as DOT_PRODUCT.  For complex types we use do
 168              a*b, not conjg(a)*b.  */
 169           pa = abase;
 170           pb = bbase;
 171           *dest = 0;
 172
 173           for (n = 0; n < count; n++)
 174             {
 175               if (*pa && *pb)
 176                 {
 177                   *dest = 1;
 178                   break;
 179                 }
 180               pa += astride;
 181               pb += bstride;
 182             }
 183
 184           dest += rxstride;
 185           abase += xstride;
 186         }
 187       abase -= xstride * xcount;
 188       bbase += ystride;
 189       dest += rystride - (rxstride * xcount);
 190     }
 191 }
 192
 193 #endif
 194 '