fflas-ffpack-common 2.2.2-5 / usr / include / fflas-ffpack / fflas / fflas_fgemm.inl

/* -*- mode: C++; tab-width: 4; indent-tabs-mode: t; c-basic-offset: 4 -*- */
// vim:sts=4:sw=4:ts=4:noet:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s

/* fflas/fflas_fgemm.inl
 * Copyright (C) 2005 Clement Pernet
 * Copyright (C) 2014 the FFLAS-FFPACK group
 *
 * Written by Clement Pernet  < Clement.Pernet@imag.fr >
 *            Brice Boyer (briceboyer) <boyer.brice@gmail.com>
 *
 *
 * ========LICENCE========
 * This file is part of the library FFLAS-FFPACK.
 *
 * FFLAS-FFPACK is free software: you can redistribute it and/or modify
 * it under the terms of the  GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 * This library 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
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 * ========LICENCE========
 *.
 */

#ifndef __FFLASFFPACK_fgemm_INL
#define __FFLASFFPACK_fgemm_INL

#include <givaro/modular.h>
#include <givaro/modular-balanced.h>
#include "fflas-ffpack/utils/debug.h"

namespace FFLAS { namespace Protected{

		template <typename FloatElement, class Field, class FieldMode>
		inline typename Field::Element_ptr
		fgemm_convert (const Field& F,
					   const FFLAS_TRANSPOSE ta,
					   const FFLAS_TRANSPOSE tb,
					   const size_t m, const size_t n, const size_t k,
					   const typename Field::Element alpha,
					   typename Field::ConstElement_ptr A,const size_t lda,
					   typename Field::ConstElement_ptr B,const size_t ldb,
					   const typename Field::Element beta,
					   typename Field::Element_ptr C, const size_t ldc,
					   MMHelper<Field, MMHelperAlgo::Winograd, FieldMode> & H)
		{
			// CP: lda, ldb, ldc can be zero (if m,n or k is 0) and since  this may have not 
			// been checked by the caller at this point.
			// FFLASFFPACK_check(lda);
			// FFLASFFPACK_check(ldb);
			// FFLASFFPACK_check(ldc);

			Givaro::ModularBalanced<FloatElement> G((FloatElement) F.characteristic());
			FloatElement tmp,alphaf, betaf;
			// This conversion is quite tricky, but convert and init are required
			// in sequence e.g. for when F is a ModularBalanced field and alpha == -1
			F.convert (tmp, beta);
			G.init(betaf, tmp);
			F.convert (tmp, alpha);
			G.init(alphaf, tmp);

			FloatElement* Af = FFLAS::fflas_new(G, m, k);
			FloatElement* Bf = FFLAS::fflas_new(G, k, n);
			FloatElement* Cf = FFLAS::fflas_new(G, m, n);

			size_t ma, ka, kb, nb; //mb, na
			if (ta == FflasTrans) { ma = k; ka = m; }
			else { ma = m; ka = k; }
			if (tb == FflasTrans) { kb = n; nb = k; }
			else {  kb = k; nb = n; }
			size_t ldaf = ka, ldbf = nb, ldcf= n;

			fconvert(F, ma, ka, Af, ka, A, lda);
			freduce(G, ma, ka, Af, ka);
			fconvert(F, kb, nb, Bf, nb, B, ldb);
			freduce(G, kb, nb, Bf, nb);

			if (!F.isZero(beta)){
				fconvert(F, m, n, Cf, n, C, ldc);
				freduce (G, m, n, Cf, n);
			}
			MMHelper<Givaro::ModularBalanced<FloatElement>, 
					 MMHelperAlgo::Winograd> 
				HG(G,H.recLevel, ParSeqHelper::Sequential());
			fgemm (G, ta, tb, m, n, k, alphaf, Af, ldaf, Bf, ldbf, betaf, Cf, ldcf, HG);

			finit (F, m, n, Cf, n, C, ldc);

			fflas_delete (Af);
			fflas_delete (Bf);
			fflas_delete (Cf);
			return C;
		}
	}//Protected
}//FFLAS

namespace FFLAS{ namespace Protected{
		template <class Field, class Element, class AlgoT, class ParSeqTrait>
		inline bool NeedPreAddReduction (Element& Outmin, Element& Outmax,
										 Element& Op1min, Element& Op1max,
										 Element& Op2min, Element& Op2max,
										 MMHelper<Field, AlgoT, ModeCategories::LazyTag, ParSeqTrait >& WH)
		{
			Outmin = Op1min + Op2min;
			Outmax = Op1max + Op2max;
			if (WH.MaxStorableValue - Op1max < Op2max ||
				WH.MaxStorableValue + Op1min < -Op2min){
				// Reducing both Op1 and Op2
				Op1min = Op2min = WH.FieldMin;
				Op1max = Op2max = WH.FieldMax;
				Outmin = 2*WH.FieldMin;
				Outmax = 2*WH.FieldMax;
				return true;
			} else return false;
		}

		template <class Field, class Element, class AlgoT, class ModeT, class ParSeqTrait>
		inline bool NeedPreAddReduction (Element& Outmin, Element& Outmax,
										 Element& Op1min, Element& Op1max,
										 Element& Op2min, Element& Op2max,
										 MMHelper<Field, AlgoT, ModeT, ParSeqTrait >& WH)
		{
			Outmin = WH.FieldMin;
			Outmax = WH.FieldMax;
			return false;
		}

		template <class Field, class Element, class AlgoT, class ParSeqTrait>
		inline bool NeedPreSubReduction (Element& Outmin, Element& Outmax,
										 Element& Op1min, Element& Op1max,
										 Element& Op2min, Element& Op2max,
										 MMHelper<Field, AlgoT, ModeCategories::LazyTag, ParSeqTrait >& WH)
		{
			Outmin = Op1min - Op2max;
			Outmax = Op1max - Op2min;
			if (WH.MaxStorableValue - Op1max < -Op2min || 
				WH.MaxStorableValue - Op2max < -Op1min){
				// Reducing both Op1 and Op2
				Op1min = Op2min = WH.FieldMin;
				Op1max = Op2max = WH.FieldMax;
				Outmin = WH.FieldMin-WH.FieldMax;
				Outmax = -Outmin;
				return true;
			} else return false;
		}

		template <class Field, class Element, class AlgoT, class ModeT, class ParSeqTrait>
		inline bool NeedPreSubReduction (Element& Outmin, Element& Outmax,
										 Element& Op1min, Element& Op1max,
										 Element& Op2min, Element& Op2max,
										 MMHelper<Field, AlgoT, ModeT, ParSeqTrait >& WH)
		{
			// Necessary? -> CP: Yes, for generic Mode of op
			Outmin = WH.FieldMin;
			Outmax = WH.FieldMax;
			return false;
		}

		//Probable bug here due to overflow of int64_t
		template<class Field, class Element, class AlgoT, class ParSeqTrait>
		inline bool NeedDoublePreAddReduction (Element& Outmin, Element& Outmax,
											   Element& Op1min, Element& Op1max,
											   Element& Op2min, Element& Op2max, Element beta,
											   MMHelper<Field, AlgoT, ModeCategories::LazyTag, ParSeqTrait >& WH)
		{
			// Testing if P5 need to be reduced
			Outmin =  std::min(beta*Op2min,beta*Op2max);
			Outmax =  std::max(beta*Op2min,beta*Op2max);
			if (Op1max > WH.MaxStorableValue-Outmax || 
				-Op1min > WH.MaxStorableValue+Outmin){
				Outmin += WH.FieldMin;
				Outmax += WH.FieldMax;
				return true;
			} else{
				Outmin += Op1min;
				Outmax += Op1max;
				return false;
			}
		}

		template<class Field, class Element, class AlgoT, class ModeT, class ParSeqTrait>
		inline bool NeedDoublePreAddReduction (Element& Outmin, Element& Outmax,
											   Element& Op1min, Element& Op1max,
											   Element& Op2min, Element& Op2max, Element beta,
											   MMHelper<Field, AlgoT, ModeT, ParSeqTrait>& WH)
		{
			Outmin = WH.FieldMin;
			Outmax = WH.FieldMax;
			return false;
		}

		template <class Field, class AlgoT, class ParSeqTrait>
		inline void ScalAndReduce (const Field& F, const size_t N,
								   const typename Field::Element alpha,
								   typename Field::Element_ptr X, const size_t incX,
								   const MMHelper<Field, AlgoT, ModeCategories::LazyTag, ParSeqTrait >& H)
		{
			if (!F.isOne(alpha) && !F.isMOne(alpha)){
				typename MMHelper<Field, AlgoT, ModeCategories::LazyTag, ParSeqTrait >::DFElt al; 
				F.convert(al, alpha);
				if (al < 0) al = -al;
				if (std::max(-H.Outmin, H.Outmax) > H.MaxStorableValue/al){
					freduce (F, N, X, incX);
					fscalin (F, N, alpha, X, incX);
				} else {
					fscalin (H.delayedField, N, alpha, X, incX);
					freduce (F, N, X, incX);
				}
			} else
				freduce (F, N, X, incX);
		}

		template <class Field, class AlgoT, class ParSeqTrait>
		inline void ScalAndReduce (const Field& F, const size_t M, const size_t N,
								   const typename Field::Element alpha,
								   typename Field::Element_ptr A, const size_t lda,
								   const MMHelper<Field, AlgoT, ModeCategories::LazyTag, ParSeqTrait >& H)
		{
			if (!F.isOne(alpha) && !F.isMOne(alpha)){
				typename MMHelper<Field, AlgoT, ModeCategories::LazyTag, ParSeqTrait >::DFElt al; 
				F.convert(al, alpha);
				if (al<0) al = -al;
				if (std::max(-H.Outmin, H.Outmax) > H.MaxStorableValue/al){
					freduce (F, M, N, A, lda);
					fscalin (F, M, N, alpha, A, lda);
				} else {
					fscalin (H.delayedField, M, N, alpha, (typename MMHelper<Field, AlgoT, ModeCategories::LazyTag, ParSeqTrait >::DFElt*)A, lda);
					freduce (F, M, N, A, lda);
				}
			} else
				freduce (F, M, N, A, lda);
		}

	} // Protected
} // FFLAS

namespace FFLAS {

	template<class Field>
	inline  typename Field::Element_ptr
	fgemm (const Field& F,
	       const FFLAS_TRANSPOSE ta,
	       const FFLAS_TRANSPOSE tb,
	       const size_t m, const size_t n, const size_t k,
	       const typename Field::Element alpha,
	       typename Field::ConstElement_ptr A, const size_t lda,
	       typename Field::ConstElement_ptr B, const size_t ldb,
	       const typename Field::Element beta,
	       typename Field::Element_ptr C, const size_t ldc,
	       MMHelper<Field, MMHelperAlgo::Winograd, ModeCategories::ConvertTo<ElementCategories::MachineFloatTag>, ParSeqHelper::Sequential> & H)
	{
		if (F.cardinality() < DOUBLE_TO_FLOAT_CROSSOVER)
			return Protected::fgemm_convert<float,Field>(F,ta,tb,m,n,k,alpha,A,lda,B,ldb,beta,C,ldc,H);
		else if (16*F.cardinality() < Givaro::ModularBalanced<double>::maxCardinality())
			return Protected::fgemm_convert<double,Field>(F,ta,tb,m,n,k,alpha,A,lda,B,ldb,beta,C,ldc,H);
		// else if (Protected::AreEqual<typename Field::Element,int64_t>::value) {
		// 	    // Stays over int64_t
		// 	MMHelper<Field, MMHelperAlgo::Winograd, ModeCategories::DelayedTag, ParSeqHelper::Sequential> HG(H);
		// 	H.Outmin=HG.Outmin;
		// 	H.Outmax=HG.Outmax;
		// 	return fgemm(F,ta,tb,m,n,k,alpha,A,lda,B,ldb,beta,C,ldc,HG);
			
		//	}
		else {
			// Fall back case: used 
			FFPACK::failure()(__func__,__LINE__,"Invalid ConvertTo Mode for this field");	
		}
		return C;
	}
}// FFLAS


// fgemm
namespace FFLAS {

	template<typename Field>
	inline typename Field::Element_ptr
	fgemm( const Field& F,
	       const FFLAS_TRANSPOSE ta,
	       const FFLAS_TRANSPOSE tb,
	       const size_t m,
	       const size_t n,
	       const size_t k,
	       const typename Field::Element alpha,
	       typename Field::ConstElement_ptr A, const size_t lda,
	       typename Field::ConstElement_ptr B, const size_t ldb,
	       const typename Field::Element beta,
	       typename Field::Element_ptr C, const size_t ldc,
	       const ParSeqHelper::Sequential seq)
	{
		MMHelper<Field, MMHelperAlgo::Auto, typename FFLAS::ModeTraits<Field>::value, ParSeqHelper::Sequential > HW (F, m, k, n, seq);
		return 	fgemm (F, ta, tb, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc, HW);
	}

	template<typename Field,class Cut,class Param>
	inline typename Field::Element_ptr
	fgemm( const Field& F,
	       const FFLAS_TRANSPOSE ta,
	       const FFLAS_TRANSPOSE tb,
	       const size_t m,
	       const size_t n,
	       const size_t k,
	       const typename Field::Element alpha,
	       typename Field::ConstElement_ptr A, const size_t lda,
	       typename Field::ConstElement_ptr B, const size_t ldb,
	       const typename Field::Element beta,
	       typename Field::Element_ptr C, const size_t ldc,
	       const ParSeqHelper::Parallel<Cut,Param> par)
	{

		MMHelper<Field, MMHelperAlgo::Auto, typename FFLAS::ModeTraits<Field>::value, ParSeqHelper::Parallel<Cut,Param> > HW (F, m, k, n, par);
		return 	fgemm (F, ta, tb, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc, HW);
	}

	template<typename Field>
	inline typename Field::Element_ptr
	fgemm( const Field& F,
	       const FFLAS_TRANSPOSE ta,
	       const FFLAS_TRANSPOSE tb,
	       const size_t m,
	       const size_t n,
	       const size_t k,
	       const typename Field::Element alpha,
	       typename Field::ConstElement_ptr A, const size_t lda,
	       typename Field::ConstElement_ptr B, const size_t ldb,
	       const typename Field::Element beta,
	       typename Field::Element_ptr C, const size_t ldc)
	{
		if (!m || !n) {return C;}
		
		if (!k || F.isZero (alpha)){
			fscalin(F, m, n, beta, C, ldc);
		 	return C;
		}
		Checker_fgemm<Field> checker(F,m,n,k,beta,C,ldc);
		fgemm(F,ta,tb,m,n,k,alpha,A,lda,B,ldb,beta,C,ldc,FFLAS::ParSeqHelper::Sequential());
		checker.check(ta,tb,alpha,A,lda,B,ldb,C);
		return C;
	}

	template<typename Field, class ModeT, class ParSeq>
	inline typename Field::Element_ptr
	fgemm( const Field& F,
	       const FFLAS_TRANSPOSE ta,
	       const FFLAS_TRANSPOSE tb,
	       const size_t m,
	       const size_t n,
	       const size_t k,
	       const typename Field::Element alpha,
	       typename Field::ConstElement_ptr A, const size_t lda,
	       typename Field::ConstElement_ptr B, const size_t ldb,
	       const typename Field::Element beta,
	       typename Field::Element_ptr C, const size_t ldc,
	       MMHelper<Field, MMHelperAlgo::Auto, ModeT, ParSeq> & H)
	{
		MMHelper<Field, typename AlgoChooser<ModeT, ParSeq>::value, ModeT, ParSeq> HW (H);
		return fgemm(F,ta,tb,m,n,k,alpha,A,lda,B,ldb,beta,C,ldc,HW);
	}

	template<class Field>
	inline typename Field::Element_ptr
	fgemm( const Field& F,
	       const FFLAS_TRANSPOSE ta,
	       const FFLAS_TRANSPOSE tb,
	       const size_t m,
	       const size_t n,
	       const size_t k,
	       const typename Field::Element alpha,
	       typename Field::ConstElement_ptr A, const size_t lda,
	       typename Field::ConstElement_ptr B, const size_t ldb,
	       const typename Field::Element beta,
	       typename Field::Element_ptr C, const size_t ldc,
	       MMHelper<Field, MMHelperAlgo::Winograd, ModeCategories::DelayedTag, ParSeqHelper::Sequential> & H)
	{		
		if (!m || !n) {return C;}

		if (!k || F.isZero (alpha)){
			fscalin(F, m, n, beta, C, ldc);
			return C;
		}
#ifndef NDEBUG
		/*  check if alpha is invertible.
		 *  XXX do it in F.isInvertible(Element&) ?
		 *  XXX do it in return status of F.inv(Element&,Element&)
		 */
		typename Field::Element e ;
		F.assign(e,beta);
		F.divin(e,alpha);
		F.mulin(e,alpha);
		FFLASFFPACK_check(F.areEqual(e,beta));
#endif

#if 0
		// detect fgemv
		if (n == 1 and ...) {}
		// detect fger
		if (k==1 and ...) {}
#endif
		if (Protected::AreEqual<Field, Givaro::Modular<double> >::value ||
		    Protected::AreEqual<Field, Givaro::ModularBalanced<double> >::value){
			//Givaro::Modular<double> need to switch to float if p too small
			if (F.characteristic() < DOUBLE_TO_FLOAT_CROSSOVER)
				return Protected::fgemm_convert<float,Field>(F,ta,tb,m,n,k,alpha,A,lda,B,ldb,beta,C,ldc,H);
		}
		if (Protected::AreEqual<Field, Givaro::Modular<int64_t> >::value ||
		    Protected::AreEqual<Field, Givaro::ModularBalanced<int64_t> >::value)
			if (16*F.cardinality() < Givaro::ModularBalanced<double>::maxCardinality())
				return Protected::fgemm_convert<double,Field>(F,ta,tb,m,n,k,alpha,A,lda,B,ldb,beta,C,ldc,H);
		
		typename Field::Element alpha_,beta_;
		if ( !F.isOne(alpha) && !F.isMOne(alpha)){
			F.assign (alpha_, F.one);
			F.div (beta_, beta, alpha);
		} else {
			F.assign (alpha_,alpha);
			F.assign (beta_,beta);
		}
		MMHelper<Field, MMHelperAlgo::Winograd, ModeCategories::LazyTag>  HD(H);
		// std::cerr<<"\n Delayed -> Lazy alpha_ = "<<alpha_<<std::endl;
		// std::cerr<<" A = "<<*A<<"\n B = "<<*B<<"\n C = "<<*C<<"\n alpha, beta ="<<alpha<<" "<<beta<<std::endl;
		fgemm (F, ta, tb, m, n, k, alpha_, A, lda, B, ldb, beta_, C, ldc, HD);
		// std::cerr<<"Sortie de fgemm Lazy C = "<<*C<<std::endl;
		Protected::ScalAndReduce (F, m, n, alpha, C, ldc, HD);
		// std::cerr<<"Sortie de ScalAndReduce C = "<<*C<<std::endl;

		H.initOut();

		return C;
	}
} // FFLAS

// #include "fflas_fgemm/matmul_algos.inl"
#include "fflas_fgemm/fgemm_classical.inl"
#include "fflas_fgemm/fgemm_winograd.inl"
// #include "fflas_fgemm/gemm_bini.inl"

// fsquare
namespace FFLAS {
	template < class Field >
	inline typename Field::Element_ptr
	fsquare (const Field& F,
			 const FFLAS_TRANSPOSE ta,
			 const size_t n, const typename Field::Element alpha,
			 typename Field::ConstElement_ptr A, const size_t lda,
			 const typename Field::Element beta,
			 typename Field::Element_ptr C, const size_t ldc)
	{

		double alphad, betad;
		F.convert (alphad, alpha);
		if (F.isMOne (beta))
			betad = -1.0;
		else
			F.convert (betad, beta);

		//! @bug why double ?
		// Double  matrices initialisation
		Givaro::DoubleDomain::Element_ptr Ad = fflas_new (Givaro::DoubleDomain(),n,n);
		Givaro::DoubleDomain::Element_ptr Cd = fflas_new (Givaro::DoubleDomain(),n,n);
		// Conversion finite Field = >  double
		fconvert (F, n, n, Ad, n, A, lda);
		if (!F.isZero(beta)) fconvert(F, n, n, Cd, n, C, ldc);

		// Call to the blas Multiplication
		FFLASFFPACK_check(n);
		cblas_dgemm (CblasRowMajor, (CBLAS_TRANSPOSE)ta,
					 (CBLAS_TRANSPOSE)ta, (int)n, (int)n, (int)n,
					 (Givaro::DoubleDomain::Element) alphad, Ad, (int)n, Ad, (int)n,
					 (Givaro::DoubleDomain::Element) betad, Cd, (int)n);
		// Conversion double = >  Finite Field
		fflas_delete (Ad);
		finit (F,n,n, Cd, n, C, ldc);
		fflas_delete (Cd);
		return C;
	}

	namespace Protected {

		// F is Modular(Balanced)<float/double>
		template < class Field >
		inline typename Field::Element_ptr
		fsquareCommon (const Field& F,
					   const FFLAS_TRANSPOSE ta,
					   const size_t n, const typename Field::Element alpha,
					   typename Field::ConstElement_ptr A, const size_t lda,
					   const typename Field::Element beta,
					   typename Field::Element_ptr C, const size_t ldc)
		{
			if (C==A) {
				typename Field::Element_ptr Ad = fflas_new (F, n, n);
				fassign(F,n,n,A,lda,Ad,n);
				fgemm (F, ta, ta, n, n, n, alpha, Ad, n, Ad, n, beta, C, ldc);
				fflas_delete (Ad);
			}
			else
				fgemm (F, ta, ta, n, n, n, alpha, A, lda, A, lda, beta, C, ldc);
			return C;

		}

	} // Protected

	template <>
	inline double* fsquare (const  Givaro::ModularBalanced<double> & F,
							const FFLAS_TRANSPOSE ta,
							const size_t n, const double alpha,
							const double* A, const size_t lda,
							const double beta,
							double* C, const size_t ldc)
	{
		return Protected::fsquareCommon(F,ta,n,alpha,A,lda,beta,C,ldc);
	}

	template <>
	inline float * fsquare (const  Givaro::ModularBalanced<float> & F,
							const FFLAS_TRANSPOSE ta,
							const size_t n, const float alpha,
							const float* A, const size_t lda,
							const float beta,
							float* C, const size_t ldc)
	{
		return Protected::fsquareCommon(F,ta,n,alpha,A,lda,beta,C,ldc);
	}

	template <>
	inline double* fsquare (const  Givaro::Modular<double> & F,
							const FFLAS_TRANSPOSE ta,
							const size_t n, const double alpha,
							const double* A, const size_t lda,
							const double beta,
							double* C, const size_t ldc)
	{
		return Protected::fsquareCommon(F,ta,n,alpha,A,lda,beta,C,ldc);
	}

	template <>
	inline float * fsquare (const  Givaro::Modular<float> & F,
							const FFLAS_TRANSPOSE ta,
							const size_t n, const float alpha,
							const float* A, const size_t lda,
							const float beta,
							float* C, const size_t ldc)
	{
		return Protected::fsquareCommon(F,ta,n,alpha,A,lda,beta,C,ldc);
	}

} // FFLAS

#endif // __FFLASFFPACK_fgemm_INL