libgfortran/generated/matmul_c4.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_c4 (gfc_array_c4 * retarray, gfc_array_c4 * a, gfc_array_c4 * b)
  42 {
  43   GFC_COMPLEX_4 *abase;
  44   GFC_COMPLEX_4 *bbase;
  45   GFC_COMPLEX_4 *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 = internal_malloc (sizeof (GFC_COMPLEX_4) * size0 (retarray));
  90       retarray->base = 0;
  91     }
  92
  93   abase = a->data;
  94   bbase = b->data;
  95   dest = retarray->data;
  96
  97   if (retarray->dim[0].stride == 0)
  98     retarray->dim[0].stride = 1;
  99   if (a->dim[0].stride == 0)
 100     a->dim[0].stride = 1;
 101   if (b->dim[0].stride == 0)
 102     b->dim[0].stride = 1;
 103
 104
 105   if (GFC_DESCRIPTOR_RANK (retarray) == 1)
 106     {
 107       /* One-dimensional result may be addressed in the code below
 108          either as a row or a column matrix. We want both cases to
 109          work. */
 110       rxstride = rystride = retarray->dim[0].stride;
 111     }
 112   else
 113     {
 114       rxstride = retarray->dim[0].stride;
 115       rystride = retarray->dim[1].stride;
 116     }
 117
 118
 119   if (GFC_DESCRIPTOR_RANK (a) == 1)
 120     {
 121       /* Treat it as a a row matrix A[1,count]. */
 122       axstride = a->dim[0].stride;
 123       aystride = 1;
 124
 125       xcount = 1;
 126       count = a->dim[0].ubound + 1 - a->dim[0].lbound;
 127     }
 128   else
 129     {
 130       axstride = a->dim[0].stride;
 131       aystride = a->dim[1].stride;
 132
 133       count = a->dim[1].ubound + 1 - a->dim[1].lbound;
 134       xcount = a->dim[0].ubound + 1 - a->dim[0].lbound;
 135     }
 136
 137   assert(count == b->dim[0].ubound + 1 - b->dim[0].lbound);
 138
 139   if (GFC_DESCRIPTOR_RANK (b) == 1)
 140     {
 141       /* Treat it as a column matrix B[count,1] */
 142       bxstride = b->dim[0].stride;
 143
 144       /* bystride should never be used for 1-dimensional b.
 145          in case it is we want it to cause a segfault, rather than
 146          an incorrect result. */
 147       bystride = 0xDEADBEEF;
 148       ycount = 1;
 149     }
 150   else
 151     {
 152       bxstride = b->dim[0].stride;
 153       bystride = b->dim[1].stride;
 154       ycount = b->dim[1].ubound + 1 - b->dim[1].lbound;
 155     }
 156
 157   assert (a->base == 0);
 158   assert (b->base == 0);
 159   assert (retarray->base == 0);
 160
 161   abase = a->data;
 162   bbase = b->data;
 163   dest = retarray->data;
 164
 165   if (rxstride == 1 && axstride == 1 && bxstride == 1)
 166     {
 167       GFC_COMPLEX_4 *bbase_y;
 168       GFC_COMPLEX_4 *dest_y;
 169       GFC_COMPLEX_4 *abase_n;
 170       GFC_COMPLEX_4 bbase_yn;
 171
 172       memset (dest, 0, (sizeof (GFC_COMPLEX_4) * size0(retarray)));
 173
 174       for (y = 0; y < ycount; y++)
 175         {
 176           bbase_y = bbase + y*bystride;
 177           dest_y = dest + y*rystride;
 178           for (n = 0; n < count; n++)
 179             {
 180               abase_n = abase + n*aystride;
 181               bbase_yn = bbase_y[n];
 182               for (x = 0; x < xcount; x++)
 183                 {
 184                   dest_y[x] += abase_n[x] * bbase_yn;
 185                 }
 186             }
 187         }
 188     }
 189   else
 190     {
 191       for (y = 0; y < ycount; y++)
 192         for (x = 0; x < xcount; x++)
 193           dest[x*rxstride + y*rystride] = (GFC_COMPLEX_4)0;
 194
 195       for (y = 0; y < ycount; y++)
 196         for (n = 0; n < count; n++)
 197           for (x = 0; x < xcount; x++)
 198             /* dest[x,y] += a[x,n] * b[n,y] */
 199             dest[x*rxstride + y*rystride] += abase[x*axstride + n*aystride] * bbase[n*bxstride + y*bystride];
 200     }
 201 }
 202