matmul_l8.c 6.8 KB

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