1 /* Implementation of the MATMUL intrinsic
2 Copyright 2002, 2005, 2006, 2007 Free Software Foundation, Inc.
3 Contributed by Paul Brook <paul@nowt.org>
5 This file is part of the GNU Fortran 95 runtime library (libgfortran).
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.
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
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.
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. */
31 #include "libgfortran.h"
37 #if defined (HAVE_GFC_INTEGER_8)
39 /* Prototype for the BLAS ?gemm subroutine, a pointer to which can be
40 passed to us by the front-end, in which case we'll call it for large
43 typedef void (*blas_call
)(const char *, const char *, const int *, const int *,
44 const int *, const GFC_INTEGER_8
*, const GFC_INTEGER_8
*,
45 const int *, const GFC_INTEGER_8
*, const int *,
46 const GFC_INTEGER_8
*, GFC_INTEGER_8
*, const int *,
49 /* The order of loops is different in the case of plain matrix
50 multiplication C=MATMUL(A,B), and in the frequent special case where
51 the argument A is the temporary result of a TRANSPOSE intrinsic:
52 C=MATMUL(TRANSPOSE(A),B). Transposed temporaries are detected by
53 looking at their strides.
55 The equivalent Fortran pseudo-code is:
57 DIMENSION A(M,COUNT), B(COUNT,N), C(M,N)
58 IF (.NOT.IS_TRANSPOSED(A)) THEN
63 C(I,J) = C(I,J)+A(I,K)*B(K,J)
74 /* If try_blas is set to a nonzero value, then the matmul function will
75 see if there is a way to perform the matrix multiplication by a call
76 to the BLAS gemm function. */
78 extern void matmul_i8 (gfc_array_i8
* const restrict retarray
,
79 gfc_array_i8
* const restrict a
, gfc_array_i8
* const restrict b
, int try_blas
,
80 int blas_limit
, blas_call gemm
);
81 export_proto(matmul_i8
);
84 matmul_i8 (gfc_array_i8
* const restrict retarray
,
85 gfc_array_i8
* const restrict a
, gfc_array_i8
* const restrict b
, int try_blas
,
86 int blas_limit
, blas_call gemm
)
88 const GFC_INTEGER_8
* restrict abase
;
89 const GFC_INTEGER_8
* restrict bbase
;
90 GFC_INTEGER_8
* restrict dest
;
92 index_type rxstride
, rystride
, axstride
, aystride
, bxstride
, bystride
;
93 index_type x
, y
, n
, count
, xcount
, ycount
;
95 assert (GFC_DESCRIPTOR_RANK (a
) == 2
96 || GFC_DESCRIPTOR_RANK (b
) == 2);
98 /* C[xcount,ycount] = A[xcount, count] * B[count,ycount]
100 Either A or B (but not both) can be rank 1:
102 o One-dimensional argument A is implicitly treated as a row matrix
103 dimensioned [1,count], so xcount=1.
105 o One-dimensional argument B is implicitly treated as a column matrix
106 dimensioned [count, 1], so ycount=1.
109 if (retarray
->data
== NULL
)
111 if (GFC_DESCRIPTOR_RANK (a
) == 1)
113 retarray
->dim
[0].lbound
= 0;
114 retarray
->dim
[0].ubound
= b
->dim
[1].ubound
- b
->dim
[1].lbound
;
115 retarray
->dim
[0].stride
= 1;
117 else if (GFC_DESCRIPTOR_RANK (b
) == 1)
119 retarray
->dim
[0].lbound
= 0;
120 retarray
->dim
[0].ubound
= a
->dim
[0].ubound
- a
->dim
[0].lbound
;
121 retarray
->dim
[0].stride
= 1;
125 retarray
->dim
[0].lbound
= 0;
126 retarray
->dim
[0].ubound
= a
->dim
[0].ubound
- a
->dim
[0].lbound
;
127 retarray
->dim
[0].stride
= 1;
129 retarray
->dim
[1].lbound
= 0;
130 retarray
->dim
[1].ubound
= b
->dim
[1].ubound
- b
->dim
[1].lbound
;
131 retarray
->dim
[1].stride
= retarray
->dim
[0].ubound
+1;
135 = internal_malloc_size (sizeof (GFC_INTEGER_8
) * size0 ((array_t
*) retarray
));
136 retarray
->offset
= 0;
140 if (GFC_DESCRIPTOR_RANK (retarray
) == 1)
142 /* One-dimensional result may be addressed in the code below
143 either as a row or a column matrix. We want both cases to
145 rxstride
= rystride
= retarray
->dim
[0].stride
;
149 rxstride
= retarray
->dim
[0].stride
;
150 rystride
= retarray
->dim
[1].stride
;
154 if (GFC_DESCRIPTOR_RANK (a
) == 1)
156 /* Treat it as a a row matrix A[1,count]. */
157 axstride
= a
->dim
[0].stride
;
161 count
= a
->dim
[0].ubound
+ 1 - a
->dim
[0].lbound
;
165 axstride
= a
->dim
[0].stride
;
166 aystride
= a
->dim
[1].stride
;
168 count
= a
->dim
[1].ubound
+ 1 - a
->dim
[1].lbound
;
169 xcount
= a
->dim
[0].ubound
+ 1 - a
->dim
[0].lbound
;
172 if (count
!= b
->dim
[0].ubound
+ 1 - b
->dim
[0].lbound
)
174 if (count
> 0 || b
->dim
[0].ubound
+ 1 - b
->dim
[0].lbound
> 0)
175 runtime_error ("dimension of array B incorrect in MATMUL intrinsic");
178 if (GFC_DESCRIPTOR_RANK (b
) == 1)
180 /* Treat it as a column matrix B[count,1] */
181 bxstride
= b
->dim
[0].stride
;
183 /* bystride should never be used for 1-dimensional b.
184 in case it is we want it to cause a segfault, rather than
185 an incorrect result. */
186 bystride
= 0xDEADBEEF;
191 bxstride
= b
->dim
[0].stride
;
192 bystride
= b
->dim
[1].stride
;
193 ycount
= b
->dim
[1].ubound
+ 1 - b
->dim
[1].lbound
;
198 dest
= retarray
->data
;
201 /* Now that everything is set up, we're performing the multiplication
204 #define POW3(x) (((float) (x)) * ((float) (x)) * ((float) (x)))
206 if (try_blas
&& rxstride
== 1 && (axstride
== 1 || aystride
== 1)
207 && (bxstride
== 1 || bystride
== 1)
208 && (((float) xcount
) * ((float) ycount
) * ((float) count
)
211 const int m
= xcount
, n
= ycount
, k
= count
, ldc
= rystride
;
212 const GFC_INTEGER_8 one
= 1, zero
= 0;
213 const int lda
= (axstride
== 1) ? aystride
: axstride
,
214 ldb
= (bxstride
== 1) ? bystride
: bxstride
;
216 if (lda
> 0 && ldb
> 0 && ldc
> 0 && m
> 1 && n
> 1 && k
> 1)
218 assert (gemm
!= NULL
);
219 gemm (axstride
== 1 ? "N" : "T", bxstride
== 1 ? "N" : "T", &m
, &n
, &k
,
220 &one
, abase
, &lda
, bbase
, &ldb
, &zero
, dest
, &ldc
, 1, 1);
225 if (rxstride
== 1 && axstride
== 1 && bxstride
== 1)
227 const GFC_INTEGER_8
* restrict bbase_y
;
228 GFC_INTEGER_8
* restrict dest_y
;
229 const GFC_INTEGER_8
* restrict abase_n
;
230 GFC_INTEGER_8 bbase_yn
;
232 if (rystride
== xcount
)
233 memset (dest
, 0, (sizeof (GFC_INTEGER_8
) * xcount
* ycount
));
236 for (y
= 0; y
< ycount
; y
++)
237 for (x
= 0; x
< xcount
; x
++)
238 dest
[x
+ y
*rystride
] = (GFC_INTEGER_8
)0;
241 for (y
= 0; y
< ycount
; y
++)
243 bbase_y
= bbase
+ y
*bystride
;
244 dest_y
= dest
+ y
*rystride
;
245 for (n
= 0; n
< count
; n
++)
247 abase_n
= abase
+ n
*aystride
;
248 bbase_yn
= bbase_y
[n
];
249 for (x
= 0; x
< xcount
; x
++)
251 dest_y
[x
] += abase_n
[x
] * bbase_yn
;
256 else if (rxstride
== 1 && aystride
== 1 && bxstride
== 1)
258 if (GFC_DESCRIPTOR_RANK (a
) != 1)
260 const GFC_INTEGER_8
*restrict abase_x
;
261 const GFC_INTEGER_8
*restrict bbase_y
;
262 GFC_INTEGER_8
*restrict dest_y
;
265 for (y
= 0; y
< ycount
; y
++)
267 bbase_y
= &bbase
[y
*bystride
];
268 dest_y
= &dest
[y
*rystride
];
269 for (x
= 0; x
< xcount
; x
++)
271 abase_x
= &abase
[x
*axstride
];
272 s
= (GFC_INTEGER_8
) 0;
273 for (n
= 0; n
< count
; n
++)
274 s
+= abase_x
[n
] * bbase_y
[n
];
281 const GFC_INTEGER_8
*restrict bbase_y
;
284 for (y
= 0; y
< ycount
; y
++)
286 bbase_y
= &bbase
[y
*bystride
];
287 s
= (GFC_INTEGER_8
) 0;
288 for (n
= 0; n
< count
; n
++)
289 s
+= abase
[n
*axstride
] * bbase_y
[n
];
290 dest
[y
*rystride
] = s
;
294 else if (axstride
< aystride
)
296 for (y
= 0; y
< ycount
; y
++)
297 for (x
= 0; x
< xcount
; x
++)
298 dest
[x
*rxstride
+ y
*rystride
] = (GFC_INTEGER_8
)0;
300 for (y
= 0; y
< ycount
; y
++)
301 for (n
= 0; n
< count
; n
++)
302 for (x
= 0; x
< xcount
; x
++)
303 /* dest[x,y] += a[x,n] * b[n,y] */
304 dest
[x
*rxstride
+ y
*rystride
] += abase
[x
*axstride
+ n
*aystride
] * bbase
[n
*bxstride
+ y
*bystride
];
306 else if (GFC_DESCRIPTOR_RANK (a
) == 1)
308 const GFC_INTEGER_8
*restrict bbase_y
;
311 for (y
= 0; y
< ycount
; y
++)
313 bbase_y
= &bbase
[y
*bystride
];
314 s
= (GFC_INTEGER_8
) 0;
315 for (n
= 0; n
< count
; n
++)
316 s
+= abase
[n
*axstride
] * bbase_y
[n
*bxstride
];
317 dest
[y
*rxstride
] = s
;
322 const GFC_INTEGER_8
*restrict abase_x
;
323 const GFC_INTEGER_8
*restrict bbase_y
;
324 GFC_INTEGER_8
*restrict dest_y
;
327 for (y
= 0; y
< ycount
; y
++)
329 bbase_y
= &bbase
[y
*bystride
];
330 dest_y
= &dest
[y
*rystride
];
331 for (x
= 0; x
< xcount
; x
++)
333 abase_x
= &abase
[x
*axstride
];
334 s
= (GFC_INTEGER_8
) 0;
335 for (n
= 0; n
< count
; n
++)
336 s
+= abase_x
[n
*aystride
] * bbase_y
[n
*bxstride
];
337 dest_y
[x
*rxstride
] = s
;