OSDN Git Service

b7ea96978008eec24eaee8c07067f3e28e54519e
[pf3gnuchains/gcc-fork.git] / 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 }