1`/* Implementation of the MATMUL intrinsic
2   Copyright (C) 2002-2018 Free Software Foundation, Inc.
3   Contributed by Paul Brook <paul@nowt.org>
4
5This file is part of the GNU Fortran runtime library (libgfortran).
6
7Libgfortran is free software; you can redistribute it and/or
8modify it under the terms of the GNU General Public
9License as published by the Free Software Foundation; either
10version 3 of the License, or (at your option) any later version.
11
12Libgfortran is distributed in the hope that it will be useful,
13but WITHOUT ANY WARRANTY; without even the implied warranty of
14MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15GNU General Public License for more details.
16
17Under Section 7 of GPL version 3, you are granted additional
18permissions described in the GCC Runtime Library Exception, version
193.1, as published by the Free Software Foundation.
20
21You should have received a copy of the GNU General Public License and
22a copy of the GCC Runtime Library Exception along with this program;
23see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
24<http://www.gnu.org/licenses/>.  */
25
26#include "libgfortran.h"
27#include <assert.h>'
28
29include(iparm.m4)dnl
30
31`#if defined (HAVE_'rtype_name`)
32
33/* Dimensions: retarray(x,y) a(x, count) b(count,y).
34   Either a or b can be rank 1.  In this case x or y is 1.  */
35
36extern void matmul_'rtype_code` ('rtype` * const restrict,
37	gfc_array_l1 * const restrict, gfc_array_l1 * const restrict);
38export_proto(matmul_'rtype_code`);
39
40void
41matmul_'rtype_code` ('rtype` * const restrict retarray,
42	gfc_array_l1 * const restrict a, gfc_array_l1 * const restrict b)
43{
44  const GFC_LOGICAL_1 * restrict abase;
45  const GFC_LOGICAL_1 * restrict bbase;
46  'rtype_name` * restrict dest;
47  index_type rxstride;
48  index_type rystride;
49  index_type xcount;
50  index_type ycount;
51  index_type xstride;
52  index_type ystride;
53  index_type x;
54  index_type y;
55  int a_kind;
56  int b_kind;
57
58  const GFC_LOGICAL_1 * restrict pa;
59  const GFC_LOGICAL_1 * restrict 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->base_addr == NULL)
69    {
70      if (GFC_DESCRIPTOR_RANK (a) == 1)
71        {
72	  GFC_DIMENSION_SET(retarray->dim[0], 0,
73	                    GFC_DESCRIPTOR_EXTENT(b,1) - 1, 1);
74        }
75      else if (GFC_DESCRIPTOR_RANK (b) == 1)
76        {
77	  GFC_DIMENSION_SET(retarray->dim[0], 0,
78	                    GFC_DESCRIPTOR_EXTENT(a,0) - 1, 1);
79        }
80      else
81        {
82	  GFC_DIMENSION_SET(retarray->dim[0], 0,
83	                    GFC_DESCRIPTOR_EXTENT(a,0) - 1, 1);
84
85          GFC_DIMENSION_SET(retarray->dim[1], 0,
86	                    GFC_DESCRIPTOR_EXTENT(b,1) - 1,
87			    GFC_DESCRIPTOR_EXTENT(retarray,0));
88        }
89
90      retarray->base_addr
91	= xmallocarray (size0 ((array_t *) retarray), sizeof ('rtype_name`));
92      retarray->offset = 0;
93    }
94    else if (unlikely (compile_options.bounds_check))
95      {
96	index_type ret_extent, arg_extent;
97
98	if (GFC_DESCRIPTOR_RANK (a) == 1)
99	  {
100	    arg_extent = GFC_DESCRIPTOR_EXTENT(b,1);
101	    ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
102	    if (arg_extent != ret_extent)
103	      runtime_error ("Incorrect extent in return array in"
104			     " MATMUL intrinsic: is %ld, should be %ld",
105			     (long int) ret_extent, (long int) arg_extent);
106	  }
107	else if (GFC_DESCRIPTOR_RANK (b) == 1)
108	  {
109	    arg_extent = GFC_DESCRIPTOR_EXTENT(a,0);
110	    ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
111	    if (arg_extent != ret_extent)
112	      runtime_error ("Incorrect extent in return array in"
113			     " MATMUL intrinsic: is %ld, should be %ld",
114			     (long int) ret_extent, (long int) arg_extent);
115	  }
116	else
117	  {
118	    arg_extent = GFC_DESCRIPTOR_EXTENT(a,0);
119	    ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
120	    if (arg_extent != ret_extent)
121	      runtime_error ("Incorrect extent in return array in"
122			     " MATMUL intrinsic for dimension 1:"
123			     " is %ld, should be %ld",
124			     (long int) ret_extent, (long int) arg_extent);
125
126	    arg_extent = GFC_DESCRIPTOR_EXTENT(b,1);
127	    ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,1);
128	    if (arg_extent != ret_extent)
129	      runtime_error ("Incorrect extent in return array in"
130			     " MATMUL intrinsic for dimension 2:"
131			     " is %ld, should be %ld",
132			     (long int) ret_extent, (long int) arg_extent);
133	  }
134      }
135
136  abase = a->base_addr;
137  a_kind = GFC_DESCRIPTOR_SIZE (a);
138
139  if (a_kind == 1 || a_kind == 2 || a_kind == 4 || a_kind == 8
140#ifdef HAVE_GFC_LOGICAL_16
141     || a_kind == 16
142#endif
143     )
144    abase = GFOR_POINTER_TO_L1 (abase, a_kind);
145  else
146    internal_error (NULL, "Funny sized logical array");
147
148  bbase = b->base_addr;
149  b_kind = GFC_DESCRIPTOR_SIZE (b);
150
151  if (b_kind == 1 || b_kind == 2 || b_kind == 4 || b_kind == 8
152#ifdef HAVE_GFC_LOGICAL_16
153     || b_kind == 16
154#endif
155     )
156    bbase = GFOR_POINTER_TO_L1 (bbase, b_kind);
157  else
158    internal_error (NULL, "Funny sized logical array");
159
160  dest = retarray->base_addr;
161'
162sinclude(`matmul_asm_'rtype_code`.m4')dnl
163`
164  if (GFC_DESCRIPTOR_RANK (retarray) == 1)
165    {
166      rxstride = GFC_DESCRIPTOR_STRIDE(retarray,0);
167      rystride = rxstride;
168    }
169  else
170    {
171      rxstride = GFC_DESCRIPTOR_STRIDE(retarray,0);
172      rystride = GFC_DESCRIPTOR_STRIDE(retarray,1);
173    }
174
175  /* If we have rank 1 parameters, zero the absent stride, and set the size to
176     one.  */
177  if (GFC_DESCRIPTOR_RANK (a) == 1)
178    {
179      astride = GFC_DESCRIPTOR_STRIDE_BYTES(a,0);
180      count = GFC_DESCRIPTOR_EXTENT(a,0);
181      xstride = 0;
182      rxstride = 0;
183      xcount = 1;
184    }
185  else
186    {
187      astride = GFC_DESCRIPTOR_STRIDE_BYTES(a,1);
188      count = GFC_DESCRIPTOR_EXTENT(a,1);
189      xstride = GFC_DESCRIPTOR_STRIDE_BYTES(a,0);
190      xcount = GFC_DESCRIPTOR_EXTENT(a,0);
191    }
192  if (GFC_DESCRIPTOR_RANK (b) == 1)
193    {
194      bstride = GFC_DESCRIPTOR_STRIDE_BYTES(b,0);
195      assert(count == GFC_DESCRIPTOR_EXTENT(b,0));
196      ystride = 0;
197      rystride = 0;
198      ycount = 1;
199    }
200  else
201    {
202      bstride = GFC_DESCRIPTOR_STRIDE_BYTES(b,0);
203      assert(count == GFC_DESCRIPTOR_EXTENT(b,0));
204      ystride = GFC_DESCRIPTOR_STRIDE_BYTES(b,1);
205      ycount = GFC_DESCRIPTOR_EXTENT(b,1);
206    }
207
208  for (y = 0; y < ycount; y++)
209    {
210      for (x = 0; x < xcount; x++)
211        {
212          /* Do the summation for this element.  For real and integer types
213             this is the same as DOT_PRODUCT.  For complex types we use do
214             a*b, not conjg(a)*b.  */
215          pa = abase;
216          pb = bbase;
217          *dest = 0;
218
219          for (n = 0; n < count; n++)
220            {
221              if (*pa && *pb)
222                {
223                  *dest = 1;
224                  break;
225                }
226              pa += astride;
227              pb += bstride;
228            }
229
230          dest += rxstride;
231          abase += xstride;
232        }
233      abase -= xstride * xcount;
234      bbase += ystride;
235      dest += rystride - (rxstride * xcount);
236    }
237}
238
239#endif
240'
241