[thirdparty/gcc.git] / libgfortran / generated / matmul_i8.c

/* Implementation of the MATMUL intrinsic
   Copyright (C) 2002-2016 Free Software Foundation, Inc.
   Contributed by Paul Brook <paul@nowt.org>

This file is part of the GNU Fortran runtime library (libgfortran).

Libgfortran is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public
License as published by the Free Software Foundation; either
version 3 of the License, or (at your option) any later version.

Libgfortran is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

Under Section 7 of GPL version 3, you are granted additional
permissions described in the GCC Runtime Library Exception, version
3.1, as published by the Free Software Foundation.

You should have received a copy of the GNU General Public License and
a copy of the GCC Runtime Library Exception along with this program;
see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
<http://www.gnu.org/licenses/>.  */

#include "libgfortran.h"
#include <stdlib.h>
#include <string.h>
#include <assert.h>


#if defined (HAVE_GFC_INTEGER_8)

/* Prototype for the BLAS ?gemm subroutine, a pointer to which can be
   passed to us by the front-end, in which case we'll call it for large
   matrices.  */

typedef void (*blas_call)(const char *, const char *, const int *, const int *,
                          const int *, const GFC_INTEGER_8 *, const GFC_INTEGER_8 *,
                          const int *, const GFC_INTEGER_8 *, const int *,
                          const GFC_INTEGER_8 *, GFC_INTEGER_8 *, const int *,
                          int, int);

/* The order of loops is different in the case of plain matrix
   multiplication C=MATMUL(A,B), and in the frequent special case where
   the argument A is the temporary result of a TRANSPOSE intrinsic:
   C=MATMUL(TRANSPOSE(A),B).  Transposed temporaries are detected by
   looking at their strides.

   The equivalent Fortran pseudo-code is:

   DIMENSION A(M,COUNT), B(COUNT,N), C(M,N)
   IF (.NOT.IS_TRANSPOSED(A)) THEN
     C = 0
     DO J=1,N
       DO K=1,COUNT
         DO I=1,M
           C(I,J) = C(I,J)+A(I,K)*B(K,J)
   ELSE
     DO J=1,N
       DO I=1,M
         S = 0
         DO K=1,COUNT
           S = S+A(I,K)*B(K,J)
         C(I,J) = S
   ENDIF
*/

/* If try_blas is set to a nonzero value, then the matmul function will
   see if there is a way to perform the matrix multiplication by a call
   to the BLAS gemm function.  */

extern void matmul_i8 (gfc_array_i8 * const restrict retarray, 
	gfc_array_i8 * const restrict a, gfc_array_i8 * const restrict b, int try_blas,
	int blas_limit, blas_call gemm);
export_proto(matmul_i8);

void
matmul_i8 (gfc_array_i8 * const restrict retarray, 
	gfc_array_i8 * const restrict a, gfc_array_i8 * const restrict b, int try_blas,
	int blas_limit, blas_call gemm)
{
  const GFC_INTEGER_8 * restrict abase;
  const GFC_INTEGER_8 * restrict bbase;
  GFC_INTEGER_8 * restrict dest;

  index_type rxstride, rystride, axstride, aystride, bxstride, bystride;
  index_type x, y, n, count, xcount, ycount;

  assert (GFC_DESCRIPTOR_RANK (a) == 2
          || GFC_DESCRIPTOR_RANK (b) == 2);

/* C[xcount,ycount] = A[xcount, count] * B[count,ycount]

   Either A or B (but not both) can be rank 1:

   o One-dimensional argument A is implicitly treated as a row matrix
     dimensioned [1,count], so xcount=1.

   o One-dimensional argument B is implicitly treated as a column matrix
     dimensioned [count, 1], so ycount=1.
  */

  if (retarray->base_addr == NULL)
    {
      if (GFC_DESCRIPTOR_RANK (a) == 1)
        {
	  GFC_DIMENSION_SET(retarray->dim[0], 0,
	                    GFC_DESCRIPTOR_EXTENT(b,1) - 1, 1);
        }
      else if (GFC_DESCRIPTOR_RANK (b) == 1)
        {
	  GFC_DIMENSION_SET(retarray->dim[0], 0,
	                    GFC_DESCRIPTOR_EXTENT(a,0) - 1, 1);
        }
      else
        {
	  GFC_DIMENSION_SET(retarray->dim[0], 0,
	                    GFC_DESCRIPTOR_EXTENT(a,0) - 1, 1);

          GFC_DIMENSION_SET(retarray->dim[1], 0,
	                    GFC_DESCRIPTOR_EXTENT(b,1) - 1,
			    GFC_DESCRIPTOR_EXTENT(retarray,0));
        }

      retarray->base_addr
	= xmallocarray (size0 ((array_t *) retarray), sizeof (GFC_INTEGER_8));
      retarray->offset = 0;
    }
    else if (unlikely (compile_options.bounds_check))
      {
	index_type ret_extent, arg_extent;

	if (GFC_DESCRIPTOR_RANK (a) == 1)
	  {
	    arg_extent = GFC_DESCRIPTOR_EXTENT(b,1);
	    ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
	    if (arg_extent != ret_extent)
	      runtime_error ("Incorrect extent in return array in"
			     " MATMUL intrinsic: is %ld, should be %ld",
			     (long int) ret_extent, (long int) arg_extent);
	  }
	else if (GFC_DESCRIPTOR_RANK (b) == 1)
	  {
	    arg_extent = GFC_DESCRIPTOR_EXTENT(a,0);
	    ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
	    if (arg_extent != ret_extent)
	      runtime_error ("Incorrect extent in return array in"
			     " MATMUL intrinsic: is %ld, should be %ld",
			     (long int) ret_extent, (long int) arg_extent);	    
	  }
	else
	  {
	    arg_extent = GFC_DESCRIPTOR_EXTENT(a,0);
	    ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
	    if (arg_extent != ret_extent)
	      runtime_error ("Incorrect extent in return array in"
			     " MATMUL intrinsic for dimension 1:"
			     " is %ld, should be %ld",
			     (long int) ret_extent, (long int) arg_extent);

	    arg_extent = GFC_DESCRIPTOR_EXTENT(b,1);
	    ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,1);
	    if (arg_extent != ret_extent)
	      runtime_error ("Incorrect extent in return array in"
			     " MATMUL intrinsic for dimension 2:"
			     " is %ld, should be %ld",
			     (long int) ret_extent, (long int) arg_extent);
	  }
      }


  if (GFC_DESCRIPTOR_RANK (retarray) == 1)
    {
      /* One-dimensional result may be addressed in the code below
	 either as a row or a column matrix. We want both cases to
	 work. */
      rxstride = rystride = GFC_DESCRIPTOR_STRIDE(retarray,0);
    }
  else
    {
      rxstride = GFC_DESCRIPTOR_STRIDE(retarray,0);
      rystride = GFC_DESCRIPTOR_STRIDE(retarray,1);
    }


  if (GFC_DESCRIPTOR_RANK (a) == 1)
    {
      /* Treat it as a a row matrix A[1,count]. */
      axstride = GFC_DESCRIPTOR_STRIDE(a,0);
      aystride = 1;

      xcount = 1;
      count = GFC_DESCRIPTOR_EXTENT(a,0);
    }
  else
    {
      axstride = GFC_DESCRIPTOR_STRIDE(a,0);
      aystride = GFC_DESCRIPTOR_STRIDE(a,1);

      count = GFC_DESCRIPTOR_EXTENT(a,1);
      xcount = GFC_DESCRIPTOR_EXTENT(a,0);
    }

  if (count != GFC_DESCRIPTOR_EXTENT(b,0))
    {
      if (count > 0 || GFC_DESCRIPTOR_EXTENT(b,0) > 0)
	runtime_error ("dimension of array B incorrect in MATMUL intrinsic");
    }

  if (GFC_DESCRIPTOR_RANK (b) == 1)
    {
      /* Treat it as a column matrix B[count,1] */
      bxstride = GFC_DESCRIPTOR_STRIDE(b,0);

      /* bystride should never be used for 1-dimensional b.
	 in case it is we want it to cause a segfault, rather than
	 an incorrect result. */
      bystride = 0xDEADBEEF;
      ycount = 1;
    }
  else
    {
      bxstride = GFC_DESCRIPTOR_STRIDE(b,0);
      bystride = GFC_DESCRIPTOR_STRIDE(b,1);
      ycount = GFC_DESCRIPTOR_EXTENT(b,1);
    }

  abase = a->base_addr;
  bbase = b->base_addr;
  dest = retarray->base_addr;


  /* Now that everything is set up, we're performing the multiplication
     itself.  */

#define POW3(x) (((float) (x)) * ((float) (x)) * ((float) (x)))

  if (try_blas && rxstride == 1 && (axstride == 1 || aystride == 1)
      && (bxstride == 1 || bystride == 1)
      && (((float) xcount) * ((float) ycount) * ((float) count)
          > POW3(blas_limit)))
  {
    const int m = xcount, n = ycount, k = count, ldc = rystride;
    const GFC_INTEGER_8 one = 1, zero = 0;
    const int lda = (axstride == 1) ? aystride : axstride,
              ldb = (bxstride == 1) ? bystride : bxstride;

    if (lda > 0 && ldb > 0 && ldc > 0 && m > 1 && n > 1 && k > 1)
      {
        assert (gemm != NULL);
        gemm (axstride == 1 ? "N" : "T", bxstride == 1 ? "N" : "T", &m, &n, &k,
              &one, abase, &lda, bbase, &ldb, &zero, dest, &ldc, 1, 1);
        return;
      }
  }

  if (rxstride == 1 && axstride == 1 && bxstride == 1)
    {
      const GFC_INTEGER_8 * restrict bbase_y;
      GFC_INTEGER_8 * restrict dest_y;
      const GFC_INTEGER_8 * restrict abase_n;
      GFC_INTEGER_8 bbase_yn;

      if (rystride == xcount)
	memset (dest, 0, (sizeof (GFC_INTEGER_8) * xcount * ycount));
      else
	{
	  for (y = 0; y < ycount; y++)
	    for (x = 0; x < xcount; x++)
	      dest[x + y*rystride] = (GFC_INTEGER_8)0;
	}

      for (y = 0; y < ycount; y++)
	{
	  bbase_y = bbase + y*bystride;
	  dest_y = dest + y*rystride;
	  for (n = 0; n < count; n++)
	    {
	      abase_n = abase + n*aystride;
	      bbase_yn = bbase_y[n];
	      for (x = 0; x < xcount; x++)
		{
		  dest_y[x] += abase_n[x] * bbase_yn;
		}
	    }
	}
    }
  else if (rxstride == 1 && aystride == 1 && bxstride == 1)
    {
      if (GFC_DESCRIPTOR_RANK (a) != 1)
	{
	  const GFC_INTEGER_8 *restrict abase_x;
	  const GFC_INTEGER_8 *restrict bbase_y;
	  GFC_INTEGER_8 *restrict dest_y;
	  GFC_INTEGER_8 s;

	  for (y = 0; y < ycount; y++)
	    {
	      bbase_y = &bbase[y*bystride];
	      dest_y = &dest[y*rystride];
	      for (x = 0; x < xcount; x++)
		{
		  abase_x = &abase[x*axstride];
		  s = (GFC_INTEGER_8) 0;
		  for (n = 0; n < count; n++)
		    s += abase_x[n] * bbase_y[n];
		  dest_y[x] = s;
		}
	    }
	}
      else
	{
	  const GFC_INTEGER_8 *restrict bbase_y;
	  GFC_INTEGER_8 s;

	  for (y = 0; y < ycount; y++)
	    {
	      bbase_y = &bbase[y*bystride];
	      s = (GFC_INTEGER_8) 0;
	      for (n = 0; n < count; n++)
		s += abase[n*axstride] * bbase_y[n];
	      dest[y*rystride] = s;
	    }
	}
    }
  else if (axstride < aystride)
    {
      for (y = 0; y < ycount; y++)
	for (x = 0; x < xcount; x++)
	  dest[x*rxstride + y*rystride] = (GFC_INTEGER_8)0;

      for (y = 0; y < ycount; y++)
	for (n = 0; n < count; n++)
	  for (x = 0; x < xcount; x++)
	    /* dest[x,y] += a[x,n] * b[n,y] */
	    dest[x*rxstride + y*rystride] += abase[x*axstride + n*aystride] * bbase[n*bxstride + y*bystride];
    }
  else if (GFC_DESCRIPTOR_RANK (a) == 1)
    {
      const GFC_INTEGER_8 *restrict bbase_y;
      GFC_INTEGER_8 s;

      for (y = 0; y < ycount; y++)
	{
	  bbase_y = &bbase[y*bystride];
	  s = (GFC_INTEGER_8) 0;
	  for (n = 0; n < count; n++)
	    s += abase[n*axstride] * bbase_y[n*bxstride];
	  dest[y*rxstride] = s;
	}
    }
  else
    {
      const GFC_INTEGER_8 *restrict abase_x;
      const GFC_INTEGER_8 *restrict bbase_y;
      GFC_INTEGER_8 *restrict dest_y;
      GFC_INTEGER_8 s;

      for (y = 0; y < ycount; y++)
	{
	  bbase_y = &bbase[y*bystride];
	  dest_y = &dest[y*rystride];
	  for (x = 0; x < xcount; x++)
	    {
	      abase_x = &abase[x*axstride];
	      s = (GFC_INTEGER_8) 0;
	      for (n = 0; n < count; n++)
		s += abase_x[n*aystride] * bbase_y[n*bxstride];
	      dest_y[x*rxstride] = s;
	    }
	}
    }
}

#endif
Commit	Line	Data
4ee9c684	1	/* Implementation of the MATMUL intrinsic
f1717362	2	Copyright (C) 2002-2016 Free Software Foundation, Inc.
4ee9c684	3	Contributed by Paul Brook <paul@nowt.org>
4ee9c684	4
553877d9	5	This file is part of the GNU Fortran runtime library (libgfortran).
4ee9c684	6
4ee9c684	7	Libgfortran is free software; you can redistribute it and/or
b417ea8c	8	modify it under the terms of the GNU General Public
4ee9c684	9	License as published by the Free Software Foundation; either
6bc9506f	10	version 3 of the License, or (at your option) any later version.
4ee9c684	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
b417ea8c	15	GNU General Public License for more details.
4ee9c684	16
6bc9506f	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/>. */
4ee9c684	25
41f2d5e8	26	#include "libgfortran.h"
4ee9c684	27	#include <stdlib.h>
4e867f90	28	#include <string.h>
4ee9c684	29	#include <assert.h>
41f2d5e8	30
4ee9c684	31
920e54ef	32	#if defined (HAVE_GFC_INTEGER_8)
920e54ef	33
4e8e57b0	34	/* Prototype for the BLAS ?gemm subroutine, a pointer to which can be
	35	passed to us by the front-end, in which case we'll call it for large
	36	matrices. */
	37
	38	typedef void (blas_call)(const char , const char , const int , const int *,
	39	const int , const GFC_INTEGER_8 , const GFC_INTEGER_8 *,
	40	const int , const GFC_INTEGER_8 , const int *,
	41	const GFC_INTEGER_8 , GFC_INTEGER_8 , const int *,
	42	int, int);
	43
018ef8b8	44	/* The order of loops is different in the case of plain matrix
	45	multiplication C=MATMUL(A,B), and in the frequent special case where
	46	the argument A is the temporary result of a TRANSPOSE intrinsic:
	47	C=MATMUL(TRANSPOSE(A),B). Transposed temporaries are detected by
	48	looking at their strides.
	49
	50	The equivalent Fortran pseudo-code is:
4e867f90	51
4e867f90	52	DIMENSION A(M,COUNT), B(COUNT,N), C(M,N)
018ef8b8	53	IF (.NOT.IS_TRANSPOSED(A)) THEN
	54	C = 0
	55	DO J=1,N
	56	DO K=1,COUNT
	57	DO I=1,M
	58	C(I,J) = C(I,J)+A(I,K)*B(K,J)
	59	ELSE
	60	DO J=1,N
4e867f90	61	DO I=1,M
018ef8b8	62	S = 0
018ef8b8	63	DO K=1,COUNT
4e8e57b0	64	S = S+A(I,K)*B(K,J)
018ef8b8	65	C(I,J) = S
018ef8b8	66	ENDIF
4e867f90	67	*/
4e867f90	68
4e8e57b0	69	/* If try_blas is set to a nonzero value, then the matmul function will
	70	see if there is a way to perform the matrix multiplication by a call
	71	to the BLAS gemm function. */
	72
02013060	73	extern void matmul_i8 (gfc_array_i8 * const restrict retarray,
4e8e57b0	74	gfc_array_i8 * const restrict a, gfc_array_i8 * const restrict b, int try_blas,
4e8e57b0	75	int blas_limit, blas_call gemm);
820b4fbd	76	export_proto(matmul_i8);
7b6cb5bd	77
4ee9c684	78	void
02013060	79	matmul_i8 (gfc_array_i8 * const restrict retarray,
4e8e57b0	80	gfc_array_i8 * const restrict a, gfc_array_i8 * const restrict b, int try_blas,
4e8e57b0	81	int blas_limit, blas_call gemm)
4ee9c684	82	{
02013060	83	const GFC_INTEGER_8 * restrict abase;
	84	const GFC_INTEGER_8 * restrict bbase;
	85	GFC_INTEGER_8 * restrict dest;
4e867f90	86
	87	index_type rxstride, rystride, axstride, aystride, bxstride, bystride;
	88	index_type x, y, n, count, xcount, ycount;
4ee9c684	89
	90	assert (GFC_DESCRIPTOR_RANK (a) == 2
	91	\|\| GFC_DESCRIPTOR_RANK (b) == 2);
4a43abf4	92
4e867f90	93	/* C[xcount,ycount] = A[xcount, count] * B[count,ycount]
	94
	95	Either A or B (but not both) can be rank 1:
	96
	97	o One-dimensional argument A is implicitly treated as a row matrix
	98	dimensioned [1,count], so xcount=1.
	99
	100	o One-dimensional argument B is implicitly treated as a column matrix
	101	dimensioned [count, 1], so ycount=1.
	102	*/
	103
553877d9	104	if (retarray->base_addr == NULL)
4a43abf4	105	{
	106	if (GFC_DESCRIPTOR_RANK (a) == 1)
	107	{
827aef63	108	GFC_DIMENSION_SET(retarray->dim[0], 0,
827aef63	109	GFC_DESCRIPTOR_EXTENT(b,1) - 1, 1);
4a43abf4	110	}
	111	else if (GFC_DESCRIPTOR_RANK (b) == 1)
	112	{
827aef63	113	GFC_DIMENSION_SET(retarray->dim[0], 0,
827aef63	114	GFC_DESCRIPTOR_EXTENT(a,0) - 1, 1);
4a43abf4	115	}
	116	else
	117	{
827aef63	118	GFC_DIMENSION_SET(retarray->dim[0], 0,
827aef63	119	GFC_DESCRIPTOR_EXTENT(a,0) - 1, 1);
c0d02c82	120
827aef63	121	GFC_DIMENSION_SET(retarray->dim[1], 0,
	122	GFC_DESCRIPTOR_EXTENT(b,1) - 1,
	123	GFC_DESCRIPTOR_EXTENT(retarray,0));
4a43abf4	124	}
c0d02c82	125
553877d9	126	retarray->base_addr
af1e9051	127	= xmallocarray (size0 ((array_t *) retarray), sizeof (GFC_INTEGER_8));
93830de1	128	retarray->offset = 0;
4a43abf4	129	}
c7fb575f	130	else if (unlikely (compile_options.bounds_check))
8dbe504b	131	{
	132	index_type ret_extent, arg_extent;
	133
	134	if (GFC_DESCRIPTOR_RANK (a) == 1)
	135	{
827aef63	136	arg_extent = GFC_DESCRIPTOR_EXTENT(b,1);
827aef63	137	ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
8dbe504b	138	if (arg_extent != ret_extent)
	139	runtime_error ("Incorrect extent in return array in"
	140	" MATMUL intrinsic: is %ld, should be %ld",
	141	(long int) ret_extent, (long int) arg_extent);
	142	}
	143	else if (GFC_DESCRIPTOR_RANK (b) == 1)
	144	{
827aef63	145	arg_extent = GFC_DESCRIPTOR_EXTENT(a,0);
827aef63	146	ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
8dbe504b	147	if (arg_extent != ret_extent)
	148	runtime_error ("Incorrect extent in return array in"
	149	" MATMUL intrinsic: is %ld, should be %ld",
	150	(long int) ret_extent, (long int) arg_extent);
	151	}
	152	else
	153	{
827aef63	154	arg_extent = GFC_DESCRIPTOR_EXTENT(a,0);
827aef63	155	ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
8dbe504b	156	if (arg_extent != ret_extent)
	157	runtime_error ("Incorrect extent in return array in"
	158	" MATMUL intrinsic for dimension 1:"
	159	" is %ld, should be %ld",
	160	(long int) ret_extent, (long int) arg_extent);
	161
827aef63	162	arg_extent = GFC_DESCRIPTOR_EXTENT(b,1);
827aef63	163	ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,1);
8dbe504b	164	if (arg_extent != ret_extent)
	165	runtime_error ("Incorrect extent in return array in"
	166	" MATMUL intrinsic for dimension 2:"
	167	" is %ld, should be %ld",
	168	(long int) ret_extent, (long int) arg_extent);
	169	}
	170	}
4a43abf4	171
4ee9c684	172
	173	if (GFC_DESCRIPTOR_RANK (retarray) == 1)
	174	{
4e867f90	175	/* One-dimensional result may be addressed in the code below
	176	either as a row or a column matrix. We want both cases to
	177	work. */
827aef63	178	rxstride = rystride = GFC_DESCRIPTOR_STRIDE(retarray,0);
4ee9c684	179	}
	180	else
	181	{
827aef63	182	rxstride = GFC_DESCRIPTOR_STRIDE(retarray,0);
827aef63	183	rystride = GFC_DESCRIPTOR_STRIDE(retarray,1);
4ee9c684	184	}
4ee9c684	185
4e867f90	186
4ee9c684	187	if (GFC_DESCRIPTOR_RANK (a) == 1)
4ee9c684	188	{
4e867f90	189	/* Treat it as a a row matrix A[1,count]. */
827aef63	190	axstride = GFC_DESCRIPTOR_STRIDE(a,0);
4e867f90	191	aystride = 1;
4e867f90	192
4ee9c684	193	xcount = 1;
827aef63	194	count = GFC_DESCRIPTOR_EXTENT(a,0);
4ee9c684	195	}
	196	else
	197	{
827aef63	198	axstride = GFC_DESCRIPTOR_STRIDE(a,0);
827aef63	199	aystride = GFC_DESCRIPTOR_STRIDE(a,1);
4e867f90	200
827aef63	201	count = GFC_DESCRIPTOR_EXTENT(a,1);
827aef63	202	xcount = GFC_DESCRIPTOR_EXTENT(a,0);
4ee9c684	203	}
4e867f90	204
827aef63	205	if (count != GFC_DESCRIPTOR_EXTENT(b,0))
710284a1	206	{
827aef63	207	if (count > 0 \|\| GFC_DESCRIPTOR_EXTENT(b,0) > 0)
710284a1	208	runtime_error ("dimension of array B incorrect in MATMUL intrinsic");
710284a1	209	}
4e867f90	210
4ee9c684	211	if (GFC_DESCRIPTOR_RANK (b) == 1)
4ee9c684	212	{
4e867f90	213	/* Treat it as a column matrix B[count,1] */
827aef63	214	bxstride = GFC_DESCRIPTOR_STRIDE(b,0);
4e867f90	215
	216	/* bystride should never be used for 1-dimensional b.
	217	in case it is we want it to cause a segfault, rather than
	218	an incorrect result. */
c0d02c82	219	bystride = 0xDEADBEEF;
4ee9c684	220	ycount = 1;
	221	}
	222	else
	223	{
827aef63	224	bxstride = GFC_DESCRIPTOR_STRIDE(b,0);
	225	bystride = GFC_DESCRIPTOR_STRIDE(b,1);
	226	ycount = GFC_DESCRIPTOR_EXTENT(b,1);
4ee9c684	227	}
4ee9c684	228
553877d9	229	abase = a->base_addr;
	230	bbase = b->base_addr;
	231	dest = retarray->base_addr;
4e867f90	232
4e8e57b0	233
	234	/* Now that everything is set up, we're performing the multiplication
	235	itself. */
	236
	237	#define POW3(x) (((float) (x)) * ((float) (x)) * ((float) (x)))
	238
	239	if (try_blas && rxstride == 1 && (axstride == 1 \|\| aystride == 1)
	240	&& (bxstride == 1 \|\| bystride == 1)
	241	&& (((float) xcount) * ((float) ycount) * ((float) count)
	242	> POW3(blas_limit)))
	243	{
	244	const int m = xcount, n = ycount, k = count, ldc = rystride;
	245	const GFC_INTEGER_8 one = 1, zero = 0;
	246	const int lda = (axstride == 1) ? aystride : axstride,
	247	ldb = (bxstride == 1) ? bystride : bxstride;
	248
	249	if (lda > 0 && ldb > 0 && ldc > 0 && m > 1 && n > 1 && k > 1)
	250	{
	251	assert (gemm != NULL);
	252	gemm (axstride == 1 ? "N" : "T", bxstride == 1 ? "N" : "T", &m, &n, &k,
	253	&one, abase, &lda, bbase, &ldb, &zero, dest, &ldc, 1, 1);
	254	return;
	255	}
	256	}
	257
4e867f90	258	if (rxstride == 1 && axstride == 1 && bxstride == 1)
4ee9c684	259	{
02013060	260	const GFC_INTEGER_8 * restrict bbase_y;
	261	GFC_INTEGER_8 * restrict dest_y;
	262	const GFC_INTEGER_8 * restrict abase_n;
4e867f90	263	GFC_INTEGER_8 bbase_yn;
4e867f90	264
407f9409	265	if (rystride == xcount)
407f9409	266	memset (dest, 0, (sizeof (GFC_INTEGER_8) * xcount * ycount));
e967d28c	267	else
	268	{
	269	for (y = 0; y < ycount; y++)
	270	for (x = 0; x < xcount; x++)
	271	dest[x + y*rystride] = (GFC_INTEGER_8)0;
	272	}
4e867f90	273
	274	for (y = 0; y < ycount; y++)
	275	{
	276	bbase_y = bbase + y*bystride;
	277	dest_y = dest + y*rystride;
	278	for (n = 0; n < count; n++)
	279	{
	280	abase_n = abase + n*aystride;
	281	bbase_yn = bbase_y[n];
	282	for (x = 0; x < xcount; x++)
	283	{
	284	dest_y[x] += abase_n[x] * bbase_yn;
	285	}
	286	}
	287	}
	288	}
018ef8b8	289	else if (rxstride == 1 && aystride == 1 && bxstride == 1)
018ef8b8	290	{
dd5b9961	291	if (GFC_DESCRIPTOR_RANK (a) != 1)
	292	{
	293	const GFC_INTEGER_8 *restrict abase_x;
	294	const GFC_INTEGER_8 *restrict bbase_y;
	295	GFC_INTEGER_8 *restrict dest_y;
	296	GFC_INTEGER_8 s;
018ef8b8	297
dd5b9961	298	for (y = 0; y < ycount; y++)
	299	{
	300	bbase_y = &bbase[y*bystride];
	301	dest_y = &dest[y*rystride];
	302	for (x = 0; x < xcount; x++)
	303	{
	304	abase_x = &abase[x*axstride];
	305	s = (GFC_INTEGER_8) 0;
	306	for (n = 0; n < count; n++)
	307	s += abase_x[n] * bbase_y[n];
	308	dest_y[x] = s;
	309	}
	310	}
	311	}
	312	else
018ef8b8	313	{
dd5b9961	314	const GFC_INTEGER_8 *restrict bbase_y;
	315	GFC_INTEGER_8 s;
	316
	317	for (y = 0; y < ycount; y++)
018ef8b8	318	{
dd5b9961	319	bbase_y = &bbase[y*bystride];
018ef8b8	320	s = (GFC_INTEGER_8) 0;
018ef8b8	321	for (n = 0; n < count; n++)
dd5b9961	322	s += abase[naxstride] bbase_y[n];
dd5b9961	323	dest[y*rystride] = s;
018ef8b8	324	}
	325	}
	326	}
	327	else if (axstride < aystride)
4e867f90	328	{
	329	for (y = 0; y < ycount; y++)
	330	for (x = 0; x < xcount; x++)
	331	dest[xrxstride + yrystride] = (GFC_INTEGER_8)0;
	332
	333	for (y = 0; y < ycount; y++)
	334	for (n = 0; n < count; n++)
	335	for (x = 0; x < xcount; x++)
	336	/* dest[x,y] += a[x,n] * b[n,y] */
	337	dest[xrxstride + yrystride] += abase[xaxstride + naystride] * bbase[nbxstride + ybystride];
4ee9c684	338	}
609b676a	339	else if (GFC_DESCRIPTOR_RANK (a) == 1)
	340	{
	341	const GFC_INTEGER_8 *restrict bbase_y;
	342	GFC_INTEGER_8 s;
	343
	344	for (y = 0; y < ycount; y++)
	345	{
	346	bbase_y = &bbase[y*bystride];
	347	s = (GFC_INTEGER_8) 0;
	348	for (n = 0; n < count; n++)
	349	s += abase[naxstride] bbase_y[n*bxstride];
	350	dest[y*rxstride] = s;
	351	}
	352	}
018ef8b8	353	else
	354	{
	355	const GFC_INTEGER_8 *restrict abase_x;
	356	const GFC_INTEGER_8 *restrict bbase_y;
	357	GFC_INTEGER_8 *restrict dest_y;
	358	GFC_INTEGER_8 s;
	359
	360	for (y = 0; y < ycount; y++)
	361	{
	362	bbase_y = &bbase[y*bystride];
	363	dest_y = &dest[y*rystride];
	364	for (x = 0; x < xcount; x++)
	365	{
	366	abase_x = &abase[x*axstride];
	367	s = (GFC_INTEGER_8) 0;
	368	for (n = 0; n < count; n++)
	369	s += abase_x[naystride] bbase_y[n*bxstride];
	370	dest_y[x*rxstride] = s;
	371	}
	372	}
	373	}
4ee9c684	374	}
920e54ef	375
920e54ef	376	#endif