fflas-ffpack-common 2.2.2-4 / usr / include / fflas-ffpack / fflas / fflas_helpers.inl

/* -*- mode: C++; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 8 -*- */
// vim:sts=8:sw=8:ts=8:noet:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
/*
 * Copyright (C) 2014 the FFLAS-FFPACK group
 *
 * Written by Clement Pernet <Clement.Pernet@imag.fr>
 *
 *
 * ========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========
 *.
 */

/** @file fflas/fflas_mmhelper.h
 * @brief  Matrix-Matrix Helper class
 */

#ifndef __FFLASFFPACK_fflas_fflas_mmhelper_INL
#define __FFLASFFPACK_fflas_fflas_mmhelper_INL

#include "fflas-ffpack/field/field-traits.h"
#include "fflas-ffpack/paladin/parallel.h"
#include "fflas-ffpack/utils/flimits.h"

#include <algorithm> // std::max

namespace FFLAS{ namespace Protected{
	/** \brief Computes the number of recursive levels to perform.
	 *
	 * \param m the common dimension in the product AxB
	 */
	template<class Field>
	int WinogradSteps (const Field & F, const size_t & m);

}//Protected
}//FFLAS

namespace FFLAS {
	
	namespace Protected{
		template <class DFE> inline size_t min_types(DFE& k) {return static_cast<size_t>(k);}
#if __FFLASFFPACK_SIZEOF_LONG == 4
		template <> inline size_t min_types(double& k) {return static_cast<size_t>(std::min(k,double(std::numeric_limits<size_t>::max())));}
		template <> inline size_t min_types(int64_t& k) {return static_cast<size_t>(std::min(k,int64_t(std::numeric_limits<size_t>::max())));}
#endif
		template <> inline size_t min_types(RecInt::rint<6>& k) {return static_cast<size_t>(uint64_t(std::min(k,RecInt::rint<6>(uint64_t(std::numeric_limits<size_t>::max())))));}
		template <> inline size_t min_types(RecInt::rint<7>& k) {return static_cast<size_t>(uint64_t(std::min(k,RecInt::rint<7>(uint64_t(std::numeric_limits<size_t>::max())))));}
		template <> inline size_t min_types(RecInt::rint<8>& k) {return static_cast<size_t>(uint64_t(std::min(k,RecInt::rint<8>(uint64_t(std::numeric_limits<size_t>::max())))));}
		template <> inline size_t min_types(RecInt::rint<9>& k) {return static_cast<size_t>(uint64_t(std::min(k,RecInt::rint<9>(uint64_t(std::numeric_limits<size_t>::max())))));}
		template <> inline size_t min_types(RecInt::rint<10>& k) {return static_cast<size_t>(uint64_t(std::min(k,RecInt::rint<10>(uint64_t(std::numeric_limits<size_t>::max())))));}
		template <> inline size_t min_types(Givaro::Integer& k) {return static_cast<size_t>(uint64_t(std::min(k,Givaro::Integer(uint64_t(std::numeric_limits<size_t>::max())))));}

		template <class T>
		inline bool unfit(T x){return false;}
		template <>
		inline bool unfit(int64_t x){return (x>limits<int32_t>::max());}
		template <size_t K>
		inline bool unfit(RecInt::rint<K> x){return (x > RecInt::rint<K>(limits<RecInt::rint<K-1>>::max()));}
		template <>
		inline bool unfit(RecInt::rint<6> x){return (x > limits<int32_t>::max());}
	}

	namespace MMHelperAlgo{
		struct Auto{};
		struct Classic{};
		struct Winograd{};
		struct WinogradPar{};
		struct Bini{};
	}

	template<class ModeT, class ParSeq>
	struct AlgoChooser{typedef MMHelperAlgo::Winograd value;};
	template<class ParSeq>
	struct AlgoChooser<ModeCategories::ConvertTo<ElementCategories::RNSElementTag>, ParSeq>{typedef MMHelperAlgo::Classic value;};

	template<class Field,
		 typename AlgoTrait = MMHelperAlgo::Auto,
		 typename ModeTrait = typename ModeTraits<Field>::value,
		 typename ParSeqTrait = ParSeqHelper::Sequential >
	struct MMHelper;
	/*! FGEMM  Helper for Default and ConvertTo modes of operation
	 */
	template<class Field,
		 typename AlgoTrait,
		 typename ParSeqTrait >
	struct MMHelper<Field, AlgoTrait, ModeCategories::DefaultTag, ParSeqTrait>
	{
		typedef MMHelper<Field,AlgoTrait, ModeCategories::DefaultTag,ParSeqTrait> Self_t;
		int recLevel ;
		ParSeqTrait parseq;

		MMHelper(){}
		MMHelper(const Field& F, size_t m, size_t k, size_t n, ParSeqTrait _PS) : recLevel(-1), parseq(_PS) {}
		MMHelper(const Field& F, int w, ParSeqTrait _PS=ParSeqTrait()) : recLevel(w), parseq(_PS) {}

		// copy constructor from other Field and Algo Traits
		template<class F2, typename AlgoT2, typename FT2, typename PS2>
		MMHelper(MMHelper<F2, AlgoT2, FT2, PS2>& WH) : recLevel(WH.recLevel), parseq(WH.parseq) {}

		friend std::ostream& operator<<(std::ostream& out, const Self_t& M)
		{
			return out <<"Helper: "
				   <<typeid(AlgoTrait).name()<<' '
				   <<typeid(ModeCategories::DefaultTag).name()<< ' '
				   << M.parseq <<std::endl
				   <<"  recLevel = "<<M.recLevel<<std::endl;
		}
	}; 
	template<class Field,
		 typename AlgoTrait,
		 typename Dest,
		 typename ParSeqTrait>
	struct MMHelper<Field, AlgoTrait, ModeCategories::ConvertTo<Dest>, ParSeqTrait>
	{
		typedef MMHelper<Field,AlgoTrait, ModeCategories::ConvertTo<Dest>,ParSeqTrait> Self_t;
		int recLevel ;
		ParSeqTrait parseq;

		MMHelper(){}
		MMHelper(const Field& F, size_t m, size_t k, size_t n, ParSeqTrait _PS) : recLevel(-1), parseq(_PS) {}
		MMHelper(const Field& F, int w, ParSeqTrait _PS=ParSeqTrait()) : recLevel(w), parseq(_PS) {}

		// copy constructor from other Field and Algo Traits
		template<class F2, typename AlgoT2, typename FT2, typename PS2>
		MMHelper(MMHelper<F2, AlgoT2, FT2, PS2>& WH) : recLevel(WH.recLevel), parseq(WH.parseq) {}

		friend std::ostream& operator<<(std::ostream& out, const Self_t& M)
		{
			return out <<"Helper: "
				   <<typeid(AlgoTrait).name()<<' '
				   <<typeid(ModeCategories::ConvertTo<Dest>).name()<< ' '
				   << M.parseq <<std::endl
				   <<"  recLevel = "<<M.recLevel<<std::endl;
		}
	}; 
	    // MMHelper for Delayed and Lazy Modes of operation
	template<class Field,
		 typename AlgoTrait,
		 typename ModeTrait,
		 typename ParSeqTrait>
	struct MMHelper {
		typedef MMHelper<Field,AlgoTrait,ModeTrait,ParSeqTrait> Self_t;
		typedef typename associatedDelayedField<const Field>::type DelayedField_t;
		typedef typename associatedDelayedField<const Field>::field DelayedField;
		typedef typename DelayedField::Element DFElt;
		int recLevel ;
		DFElt FieldMin, FieldMax, Amin, Amax, Bmin, Bmax, Cmin, Cmax, Outmin, Outmax;
		DFElt MaxStorableValue;
	
		const DelayedField_t delayedField;
		ParSeqTrait parseq;
		void initC(){Cmin = FieldMin; Cmax = FieldMax;}
		void initA(){Amin = FieldMin; Amax = FieldMax;}
		void initB(){Bmin = FieldMin; Bmax = FieldMax;}
		void initOut(){Outmin = FieldMin; Outmax = FieldMax;}


		size_t MaxDelayedDim(DFElt beta)
		{
			if (MaxStorableValue < DFElt(0))
				//Infinte precision delayed field
				return std::numeric_limits<size_t>::max();

			DFElt absbeta;
			delayedField.init(absbeta,beta);
			if (beta < 0) absbeta = -beta;
			// This cast is needed when Cmin base type is int8/16_t,
			// getting -Cmin returns a int, not the same base type.
			DFElt diff = MaxStorableValue - absbeta
				* std::max(static_cast<const DFElt&>(-Cmin), Cmax);
			DFElt AB = std::max(static_cast<const DFElt&>(-Amin), Amax)
				* std::max(static_cast<const DFElt&>(-Bmin), Bmax);
			if ((diff < DFElt(0u))||(AB<DFElt(0u))) return 0;


			DFElt kmax = diff/AB;
			return FFLAS::Protected::min_types<DFElt>(kmax);
			// if (kmax > std::numeric_limits<size_t>::max())
			// 	return std::numeric_limits<size_t>::max();
			// else
			// 	return kmax;
		}
		bool Aunfit(){ return Protected::unfit(std::max(static_cast<const DFElt&>(-Amin),Amax));}
		bool Bunfit(){ return Protected::unfit(std::max(static_cast<const DFElt&>(-Bmin),Bmax));}
		void setOutBounds(const size_t k, const DFElt alpha, const DFElt beta)
		{
			if (beta<0){
				Outmin = beta*Cmax;
				Outmax = beta*Cmin;
			} else {
				Outmin = beta*Cmin;
				Outmax = beta*Cmax;
			}
			if (alpha >0){
				Outmin += DFElt(k)*alpha*std::min(Amin*Bmax, Amax*Bmin);
				Outmax += DFElt(k)*alpha*std::max(Amin*Bmin, Amax*Bmax);
			}else{
				Outmin += DFElt(k)*alpha*std::max(Amin*Bmin, Amax*Bmax);
				Outmax += DFElt(k)*alpha*std::min(Amin*Bmax, Amax*Bmin);
			}
		}

		bool checkA(const Field& F, const FFLAS::FFLAS_TRANSPOSE ta, const size_t M, const size_t N,
			    typename Field::ConstElement_ptr A, const size_t lda )
		{
#ifdef DEBUG
			for (size_t i=0; i<M;++i)
				for (size_t j=0; j<N;++j){
					const typename Field::Element x = (ta == FFLAS::FflasNoTrans)? A[i*lda+j] : A[i+j*lda];
					if (x > Amax || x < Amin){
						std::cerr<<"Error in "<<Amin<<" < = A["<<i<<", "<<j<<"] ="<<x<<" <= "<<Amax<<std::endl;
						return false;
					}
				}
#endif
			return true;
		}

		bool checkB(const Field& F, const FFLAS::FFLAS_TRANSPOSE tb, const size_t M, const size_t N,
			    typename Field::ConstElement_ptr B, const size_t ldb)
		{
#ifdef DEBUG
			for (size_t i=0; i<M;++i)
				for (size_t j=0; j<N;++j){
					const typename Field::Element x = (tb == FFLAS::FflasNoTrans)? B[i*ldb+j] : B[i+j*ldb];
					if (x > Bmax || x < Bmin){
						std::cerr<<"Error in "<<Bmin<<" < = B["<<i<<", "<<j<<"] ="<<B[i*ldb+j]<<" <= "<<Bmax<<std::endl;
						return false;
					}
				}
#endif
			return true;
		}

		bool checkOut(const Field& F, const size_t M, const size_t N,
			      typename Field::ConstElement_ptr A, const size_t lda ){
#ifdef DEBUG
			for (size_t i=0; i<M;++i)
				for (size_t j=0; j<N;++j)
					if ((A[i*lda+j]>Outmax) || (A[i*lda+j]<Outmin)){
						std::cerr<<"Error in "<<Outmin<<" <= Out["<<i<<", "<<j<<"] = "<<A[i*lda+j]<<" <= "<<Outmax<<std::endl;
						return false;
					}
#endif
			return true;
		}

		MMHelper(){}
		//TODO: delayedField constructor has a >0 characteristic even when it is a Double/FloatDomain
		// correct but semantically not satisfactory
		MMHelper(const Field& F, size_t m, size_t k, size_t n, ParSeqTrait _PS) :
			recLevel(-1),
			FieldMin((DFElt)F.minElement()), FieldMax((DFElt)F.maxElement()),
			Amin(FieldMin), Amax(FieldMax),
			Bmin(FieldMin), Bmax(FieldMax),
			Cmin(FieldMin), Cmax(FieldMax),
			Outmin(0), Outmax(0),
			MaxStorableValue ((DFElt)(limits<typename DelayedField::Element>::max())),
			delayedField(F),
			// delayedField((typename Field::Element)F.characteristic()),
			parseq(_PS)
		{
		}

		MMHelper(const Field& F, int w, ParSeqTrait _PS=ParSeqTrait()) :
			recLevel(w), 
			FieldMin((DFElt)F.minElement()), FieldMax((DFElt)F.maxElement()),
			Amin(FieldMin), Amax(FieldMax),
			Bmin(FieldMin), Bmax(FieldMax),
			Cmin(FieldMin), Cmax(FieldMax),
			Outmin(0), Outmax(0),
			MaxStorableValue ((DFElt)(limits<typename DelayedField::Element>::max())),
			delayedField(F),
			parseq(_PS)
	        {
		}

		// copy constructor from other Field and Algo Traits
		template<class F2, typename AlgoT2, typename FT2, typename PS2>
		MMHelper(MMHelper<F2, AlgoT2, FT2, PS2>& WH) :
			recLevel(WH.recLevel),
			FieldMin(WH.FieldMin), FieldMax(WH.FieldMax),
			Amin(WH.Amin), Amax(WH.Amax),
			Bmin(WH.Bmin), Bmax(WH.Bmax),
			Cmin(WH.Cmin), Cmax(WH.Cmax),
			Outmin(WH.Outmin), Outmax(WH.Outmax),
			MaxStorableValue(WH.MaxStorableValue),
			delayedField(WH.delayedField),
			parseq(WH.parseq)
		{
		}

		MMHelper(const Field& F, int w,
			 DFElt _Amin, DFElt _Amax,
			 DFElt _Bmin, DFElt _Bmax,
			 DFElt _Cmin, DFElt _Cmax,
             ParSeqTrait _PS=ParSeqTrait()):
			recLevel(w), FieldMin((DFElt)F.minElement()), FieldMax((DFElt)F.maxElement()),
			Amin(_Amin), Amax(_Amax),
			Bmin(_Bmin), Bmax(_Bmax),
			Cmin(_Cmin), Cmax(_Cmax),
			Outmin(0),Outmax(0),
			MaxStorableValue(limits<typename DelayedField::Element>::max()),
			delayedField(F),
            parseq(_PS)
		{
		}

		friend std::ostream& operator<<(std::ostream& out, const Self_t& M)
		{
			return out <<"Helper: "
				   <<typeid(AlgoTrait).name()<<' '
				   <<typeid(ModeTrait).name()<< ' '
				   << M.parseq <<std::endl
				   <<" DelayedField = "<<typeid(DelayedField).name()<<std::endl
				   <<"  recLevel = "<<M.recLevel<<std::endl
				   <<"  FieldMin = "<<M.FieldMin<<" FieldMax = "<<M.FieldMax<<std::endl
				   <<"  MaxStorableValue = "<< M.MaxStorableValue <<std::endl
				   <<"  Amin = "<<M.Amin<<" Amax = "<<M.Amax<<std::endl
				   <<"  Bmin = "<<M.Bmin<<" Bmax = "<<M.Bmax<<std::endl
				   <<"  Cmin = "<<M.Cmin<<" Cmax = "<<M.Cmax<<std::endl
				   <<"  Outmin = "<<M.Outmin<<" Outmax = "<<M.Outmax<<std::endl;
		}
	}; // MMHelper


	    // to be used in the future, when Winograd's algorithm will be made generic wrt the ModeTrait
	// template <class Field, class AlgoT, class ParSeqH>
	// void copyOutBounds(const MMHelper<Field,AlgoT,ModeCategories::DelayedTag, ParSeqH> &Source,
	// 		   MMHelper<Field,AlgoT,ModeCategories::DelayedTag, ParSeqH> & Dest){
	// 	Dest.Outmax = Source.Outmax;
	// 	Dest.Outmin = Source.Outmin;
	// }
	// template <class Field, class AlgoT, class ParSeqH>
	// void copyOutBounds(const MMHelper<Field,AlgoT,ModeCategories::LazyTag, ParSeqH> &Source,
	// 		   MMHelper<Field,AlgoT,ModeCategories::LazyTag, ParSeqH> & Dest){
	// 	Dest.Outmax = Source.Outmax;
	// 	Dest.Outmin = Source.Outmin;
	// }
	// template <class MMH1, class MMH2>
	// void copyOutBounds(const MMH1 &Source, MMH2 & Dest){}
	/*! StructureHelper for ftrsm
	*/
	namespace StructureHelper {
		struct Recursive{};
		struct Iterative{};
		struct Hybrid{};
	}

	/*! TRSM Helper
	*/
	template<typename RecIterTrait = StructureHelper::Recursive, typename ParSeqTrait = ParSeqHelper::Sequential>
	struct TRSMHelper {
		ParSeqTrait parseq;
		template<class Cut,class Param>
		TRSMHelper(ParSeqHelper::Parallel<Cut,Param> _PS):parseq(_PS){}
		TRSMHelper(ParSeqHelper::Sequential _PS):parseq(_PS){}
		template<typename RIT, typename PST>
		TRSMHelper(TRSMHelper<RIT,PST>& _TH):parseq(_TH.parseq){}

		template<class Dom, class Algo=FFLAS::MMHelperAlgo::Winograd, class ModeT=typename FFLAS::ModeTraits<Dom>::value>
		FFLAS::MMHelper<Dom, Algo, ModeT, ParSeqTrait> pMMH (Dom& D, size_t m, size_t k, size_t n, ParSeqTrait p) const {
			return FFLAS::MMHelper<Dom, Algo, ModeT, ParSeqTrait>(D,m,k,n,p);
		}

		template<class Dom, class Algo=FFLAS::MMHelperAlgo::Winograd, class ModeT=typename FFLAS::ModeTraits<Dom>::value>
		FFLAS::MMHelper<Dom, Algo, ModeT, ParSeqTrait> pMMH (Dom& D, size_t m, size_t k, size_t n) const {
			return pMMH(D,m,k,n,this->parseq);
		}

	};




} // FFLAS
#endif