OSDN Git Service

2009-06-21 Thomas Koenig <tkoenig@gcc.gnu.org>
[pf3gnuchains/gcc-fork.git] / libgfortran / generated / matmul_l8.c
1 /* Implementation of the MATMUL intrinsic
2    Copyright 2002, 2005, 2006, 2007, 2009 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 3 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 General Public License for more details.
16
17 Under Section 7 of GPL version 3, you are granted additional
18 permissions described in the GCC Runtime Library Exception, version
19 3.1, as published by the Free Software Foundation.
20
21 You should have received a copy of the GNU General Public License and
22 a copy of the GCC Runtime Library Exception along with this program;
23 see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
24 <http://www.gnu.org/licenses/>.  */
25
26 #include "libgfortran.h"
27 #include <stdlib.h>
28 #include <assert.h>
29
30
31 #if defined (HAVE_GFC_LOGICAL_8)
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
36 extern void matmul_l8 (gfc_array_l8 * const restrict, 
37         gfc_array_l1 * const restrict, gfc_array_l1 * const restrict);
38 export_proto(matmul_l8);
39
40 void
41 matmul_l8 (gfc_array_l8 * 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   GFC_LOGICAL_8 * 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->data == 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->data
91         = internal_malloc_size (sizeof (GFC_LOGICAL_8) * size0 ((array_t *) retarray));
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->data;
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->data;
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->data;
161
162
163   if (GFC_DESCRIPTOR_RANK (retarray) == 1)
164     {
165       rxstride = GFC_DESCRIPTOR_STRIDE(retarray,0);
166       rystride = rxstride;
167     }
168   else
169     {
170       rxstride = GFC_DESCRIPTOR_STRIDE(retarray,0);
171       rystride = GFC_DESCRIPTOR_STRIDE(retarray,1);
172     }
173
174   /* If we have rank 1 parameters, zero the absent stride, and set the size to
175      one.  */
176   if (GFC_DESCRIPTOR_RANK (a) == 1)
177     {
178       astride = GFC_DESCRIPTOR_STRIDE_BYTES(a,0);
179       count = GFC_DESCRIPTOR_EXTENT(a,0);
180       xstride = 0;
181       rxstride = 0;
182       xcount = 1;
183     }
184   else
185     {
186       astride = GFC_DESCRIPTOR_STRIDE_BYTES(a,1);
187       count = GFC_DESCRIPTOR_EXTENT(a,1);
188       xstride = GFC_DESCRIPTOR_STRIDE_BYTES(a,0);
189       xcount = GFC_DESCRIPTOR_EXTENT(a,0);
190     }
191   if (GFC_DESCRIPTOR_RANK (b) == 1)
192     {
193       bstride = GFC_DESCRIPTOR_STRIDE_BYTES(b,0);
194       assert(count == GFC_DESCRIPTOR_EXTENT(b,0));
195       ystride = 0;
196       rystride = 0;
197       ycount = 1;
198     }
199   else
200     {
201       bstride = GFC_DESCRIPTOR_STRIDE_BYTES(b,0);
202       assert(count == GFC_DESCRIPTOR_EXTENT(b,0));
203       ystride = GFC_DESCRIPTOR_STRIDE_BYTES(b,1);
204       ycount = GFC_DESCRIPTOR_EXTENT(b,1);
205     }
206
207   for (y = 0; y < ycount; y++)
208     {
209       for (x = 0; x < xcount; x++)
210         {
211           /* Do the summation for this element.  For real and integer types
212              this is the same as DOT_PRODUCT.  For complex types we use do
213              a*b, not conjg(a)*b.  */
214           pa = abase;
215           pb = bbase;
216           *dest = 0;
217
218           for (n = 0; n < count; n++)
219             {
220               if (*pa && *pb)
221                 {
222                   *dest = 1;
223                   break;
224                 }
225               pa += astride;
226               pb += bstride;
227             }
228
229           dest += rxstride;
230           abase += xstride;
231         }
232       abase -= xstride * xcount;
233       bbase += ystride;
234       dest += rystride - (rxstride * xcount);
235     }
236 }
237
238 #endif
239