1 /* Implementation of the MATMUL intrinsic
2 Copyright 2002, 2005 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. */
35 #include "libgfortran.h"
37 #if defined (HAVE_GFC_INTEGER_8)
39 /* The order of loops is different in the case of plain matrix
40 multiplication C=MATMUL(A,B), and in the frequent special case where
41 the argument A is the temporary result of a TRANSPOSE intrinsic:
42 C=MATMUL(TRANSPOSE(A),B). Transposed temporaries are detected by
43 looking at their strides.
45 The equivalent Fortran pseudo-code is:
47 DIMENSION A(M,COUNT), B(COUNT,N), C(M,N)
48 IF (.NOT.IS_TRANSPOSED(A)) THEN
53 C(I,J) = C(I,J)+A(I,K)*B(K,J)
64 extern void matmul_i8 (gfc_array_i8
* const restrict retarray
,
65 gfc_array_i8
* const restrict a
, gfc_array_i8
* const restrict b
);
66 export_proto(matmul_i8
);
69 matmul_i8 (gfc_array_i8
* const restrict retarray
,
70 gfc_array_i8
* const restrict a
, gfc_array_i8
* const restrict b
)
72 const GFC_INTEGER_8
* restrict abase
;
73 const GFC_INTEGER_8
* restrict bbase
;
74 GFC_INTEGER_8
* restrict dest
;
76 index_type rxstride
, rystride
, axstride
, aystride
, bxstride
, bystride
;
77 index_type x
, y
, n
, count
, xcount
, ycount
;
79 assert (GFC_DESCRIPTOR_RANK (a
) == 2
80 || GFC_DESCRIPTOR_RANK (b
) == 2);
82 /* C[xcount,ycount] = A[xcount, count] * B[count,ycount]
84 Either A or B (but not both) can be rank 1:
86 o One-dimensional argument A is implicitly treated as a row matrix
87 dimensioned [1,count], so xcount=1.
89 o One-dimensional argument B is implicitly treated as a column matrix
90 dimensioned [count, 1], so ycount=1.
93 if (retarray
->data
== NULL
)
95 if (GFC_DESCRIPTOR_RANK (a
) == 1)
97 retarray
->dim
[0].lbound
= 0;
98 retarray
->dim
[0].ubound
= b
->dim
[1].ubound
- b
->dim
[1].lbound
;
99 retarray
->dim
[0].stride
= 1;
101 else if (GFC_DESCRIPTOR_RANK (b
) == 1)
103 retarray
->dim
[0].lbound
= 0;
104 retarray
->dim
[0].ubound
= a
->dim
[0].ubound
- a
->dim
[0].lbound
;
105 retarray
->dim
[0].stride
= 1;
109 retarray
->dim
[0].lbound
= 0;
110 retarray
->dim
[0].ubound
= a
->dim
[0].ubound
- a
->dim
[0].lbound
;
111 retarray
->dim
[0].stride
= 1;
113 retarray
->dim
[1].lbound
= 0;
114 retarray
->dim
[1].ubound
= b
->dim
[1].ubound
- b
->dim
[1].lbound
;
115 retarray
->dim
[1].stride
= retarray
->dim
[0].ubound
+1;
119 = internal_malloc_size (sizeof (GFC_INTEGER_8
) * size0 ((array_t
*) retarray
));
120 retarray
->offset
= 0;
123 if (retarray
->dim
[0].stride
== 0)
124 retarray
->dim
[0].stride
= 1;
126 /* This prevents constifying the input arguments. */
127 if (a
->dim
[0].stride
== 0)
128 a
->dim
[0].stride
= 1;
129 if (b
->dim
[0].stride
== 0)
130 b
->dim
[0].stride
= 1;
133 if (GFC_DESCRIPTOR_RANK (retarray
) == 1)
135 /* One-dimensional result may be addressed in the code below
136 either as a row or a column matrix. We want both cases to
138 rxstride
= rystride
= retarray
->dim
[0].stride
;
142 rxstride
= retarray
->dim
[0].stride
;
143 rystride
= retarray
->dim
[1].stride
;
147 if (GFC_DESCRIPTOR_RANK (a
) == 1)
149 /* Treat it as a a row matrix A[1,count]. */
150 axstride
= a
->dim
[0].stride
;
154 count
= a
->dim
[0].ubound
+ 1 - a
->dim
[0].lbound
;
158 axstride
= a
->dim
[0].stride
;
159 aystride
= a
->dim
[1].stride
;
161 count
= a
->dim
[1].ubound
+ 1 - a
->dim
[1].lbound
;
162 xcount
= a
->dim
[0].ubound
+ 1 - a
->dim
[0].lbound
;
165 assert(count
== b
->dim
[0].ubound
+ 1 - b
->dim
[0].lbound
);
167 if (GFC_DESCRIPTOR_RANK (b
) == 1)
169 /* Treat it as a column matrix B[count,1] */
170 bxstride
= b
->dim
[0].stride
;
172 /* bystride should never be used for 1-dimensional b.
173 in case it is we want it to cause a segfault, rather than
174 an incorrect result. */
175 bystride
= 0xDEADBEEF;
180 bxstride
= b
->dim
[0].stride
;
181 bystride
= b
->dim
[1].stride
;
182 ycount
= b
->dim
[1].ubound
+ 1 - b
->dim
[1].lbound
;
187 dest
= retarray
->data
;
189 if (rxstride
== 1 && axstride
== 1 && bxstride
== 1)
191 const GFC_INTEGER_8
* restrict bbase_y
;
192 GFC_INTEGER_8
* restrict dest_y
;
193 const GFC_INTEGER_8
* restrict abase_n
;
194 GFC_INTEGER_8 bbase_yn
;
196 if (rystride
== ycount
)
197 memset (dest
, 0, (sizeof (GFC_INTEGER_8
) * size0((array_t
*) retarray
)));
200 for (y
= 0; y
< ycount
; y
++)
201 for (x
= 0; x
< xcount
; x
++)
202 dest
[x
+ y
*rystride
] = (GFC_INTEGER_8
)0;
205 for (y
= 0; y
< ycount
; y
++)
207 bbase_y
= bbase
+ y
*bystride
;
208 dest_y
= dest
+ y
*rystride
;
209 for (n
= 0; n
< count
; n
++)
211 abase_n
= abase
+ n
*aystride
;
212 bbase_yn
= bbase_y
[n
];
213 for (x
= 0; x
< xcount
; x
++)
215 dest_y
[x
] += abase_n
[x
] * bbase_yn
;
220 else if (rxstride
== 1 && aystride
== 1 && bxstride
== 1)
222 const GFC_INTEGER_8
*restrict abase_x
;
223 const GFC_INTEGER_8
*restrict bbase_y
;
224 GFC_INTEGER_8
*restrict dest_y
;
227 for (y
= 0; y
< ycount
; y
++)
229 bbase_y
= &bbase
[y
*bystride
];
230 dest_y
= &dest
[y
*rystride
];
231 for (x
= 0; x
< xcount
; x
++)
233 abase_x
= &abase
[x
*axstride
];
234 s
= (GFC_INTEGER_8
) 0;
235 for (n
= 0; n
< count
; n
++)
236 s
+= abase_x
[n
] * bbase_y
[n
];
241 else if (axstride
< aystride
)
243 for (y
= 0; y
< ycount
; y
++)
244 for (x
= 0; x
< xcount
; x
++)
245 dest
[x
*rxstride
+ y
*rystride
] = (GFC_INTEGER_8
)0;
247 for (y
= 0; y
< ycount
; y
++)
248 for (n
= 0; n
< count
; n
++)
249 for (x
= 0; x
< xcount
; x
++)
250 /* dest[x,y] += a[x,n] * b[n,y] */
251 dest
[x
*rxstride
+ y
*rystride
] += abase
[x
*axstride
+ n
*aystride
] * bbase
[n
*bxstride
+ y
*bystride
];
255 const GFC_INTEGER_8
*restrict abase_x
;
256 const GFC_INTEGER_8
*restrict bbase_y
;
257 GFC_INTEGER_8
*restrict dest_y
;
260 for (y
= 0; y
< ycount
; y
++)
262 bbase_y
= &bbase
[y
*bystride
];
263 dest_y
= &dest
[y
*rystride
];
264 for (x
= 0; x
< xcount
; x
++)
266 abase_x
= &abase
[x
*axstride
];
267 s
= (GFC_INTEGER_8
) 0;
268 for (n
= 0; n
< count
; n
++)
269 s
+= abase_x
[n
*aystride
] * bbase_y
[n
*bxstride
];
270 dest_y
[x
*rxstride
] = s
;