//
//  Copyright (c) 2018-2019, Cem Bassoy, cem.bassoy@gmail.com
//
//  Distributed under the Boost Software License, Version 1.0. (See
//  accompanying file LICENSE_1_0.txt or copy at
//  http://www.boost.org/LICENSE_1_0.txt)
//
//  The authors gratefully acknowledge the support of
//  Fraunhofer IOSB, Ettlingen, Germany
//


#ifndef BOOST_UBLAS_TENSOR_MULTIPLICATION
#define BOOST_UBLAS_TENSOR_MULTIPLICATION

#include <cassert>

namespace boost {
namespace numeric {
namespace ublas {
namespace detail {
namespace recursive {


/** @brief Computes the tensor-times-tensor product for q contraction modes
 *
 * Implements C[i1,...,ir,j1,...,js] = sum( A[i1,...,ir+q] * B[j1,...,js+q]  )
 *
 * nc[x]         = na[phia[x]  ] for 1 <= x <= r
 * nc[r+x]       = nb[phib[x]  ] for 1 <= x <= s
 * na[phia[r+x]] = nb[phib[s+x]] for 1 <= x <= q
 *
 * @note is used in function ttt
 *
 * @param k  zero-based recursion level starting with 0
 * @param r  number of non-contraction indices of A
 * @param s  number of non-contraction indices of B
 * @param q  number of contraction indices with q > 0
 * @param phia pointer to the permutation tuple of length q+r for A
 * @param phib pointer to the permutation tuple of length q+s for B
 * @param c  pointer to the output tensor C with rank(A)=r+s
 * @param nc pointer to the extents of tensor C
 * @param wc pointer to the strides of tensor C
 * @param a  pointer to the first input tensor with rank(A)=r+q
 * @param na pointer to the extents of the first input tensor A
 * @param wa pointer to the strides of the first input tensor A
 * @param b  pointer to the second input tensor B with rank(B)=s+q
 * @param nb pointer to the extents of the second input tensor B
 * @param wb pointer to the strides of the second input tensor B
*/

template <class PointerOut, class PointerIn1, class PointerIn2, class SizeType>
void ttt(SizeType const k,
         SizeType const r, SizeType const s, SizeType const q,
         SizeType const*const phia, SizeType const*const phib,
         PointerOut c, SizeType const*const nc, SizeType const*const wc,
         PointerIn1 a, SizeType const*const na, SizeType const*const wa,
         PointerIn2 b, SizeType const*const nb, SizeType const*const wb)
{
	if(k < r)
	{
		assert(nc[k] == na[phia[k]-1]);
		for(size_t ic = 0u; ic < nc[k]; a += wa[phia[k]-1], c += wc[k], ++ic)
			ttt(k+1, r, s, q,  phia,phib,  c, nc, wc,   a, na, wa,   b, nb, wb);
	}
	else if(k < r+s)
	{
		assert(nc[k] == nb[phib[k-r]-1]);
		for(size_t ic = 0u; ic < nc[k]; b += wb[phib[k-r]-1], c += wc[k], ++ic)
			ttt(k+1, r, s, q,  phia, phib,    c, nc, wc,   a, na, wa,   b, nb, wb);
	}
	else if(k < r+s+q-1)
	{
		assert(na[phia[k-s]-1] == nb[phib[k-r]-1]);
		for(size_t ia = 0u; ia < na[phia[k-s]-1]; a += wa[phia[k-s]-1], b += wb[phib[k-r]-1], ++ia)
			ttt(k+1, r, s, q,  phia, phib,  c, nc, wc,   a, na, wa,   b, nb, wb);
	}
	else
	{
		assert(na[phia[k-s]-1] == nb[phib[k-r]-1]);
		for(size_t ia = 0u; ia < na[phia[k-s]-1]; a += wa[phia[k-s]-1], b += wb[phib[k-r]-1], ++ia)
			*c += *a * *b;
	}
}




/** @brief Computes the tensor-times-tensor product for q contraction modes
 *
 * Implements C[i1,...,ir,j1,...,js] = sum( A[i1,...,ir+q] * B[j1,...,js+q]  )
 *
 * @note no permutation tuple is used
 *
 * nc[x]   = na[x  ] for 1 <= x <= r
 * nc[r+x] = nb[x  ] for 1 <= x <= s
 * na[r+x] = nb[s+x] for 1 <= x <= q
 *
 * @note is used in function ttt
 *
 * @param k  zero-based recursion level starting with 0
 * @param r  number of non-contraction indices of A
 * @param s  number of non-contraction indices of B
 * @param q  number of contraction indices with q > 0
 * @param c  pointer to the output tensor C with rank(A)=r+s
 * @param nc pointer to the extents of tensor C
 * @param wc pointer to the strides of tensor C
 * @param a  pointer to the first input tensor with rank(A)=r+q
 * @param na pointer to the extents of the first input tensor A
 * @param wa pointer to the strides of the first input tensor A
 * @param b  pointer to the second input tensor B with rank(B)=s+q
 * @param nb pointer to the extents of the second input tensor B
 * @param wb pointer to the strides of the second input tensor B
*/

template <class PointerOut, class PointerIn1, class PointerIn2, class SizeType>
void ttt(SizeType const k,
         SizeType const r, SizeType const s, SizeType const q,
         PointerOut c, SizeType const*const nc, SizeType const*const wc,
         PointerIn1 a, SizeType const*const na, SizeType const*const wa,
         PointerIn2 b, SizeType const*const nb, SizeType const*const wb)
{
	if(k < r)
	{
		assert(nc[k] == na[k]);
		for(size_t ic = 0u; ic < nc[k]; a += wa[k], c += wc[k], ++ic)
			ttt(k+1, r, s, q,  c, nc, wc,   a, na, wa,   b, nb, wb);
	}
	else if(k < r+s)
	{
		assert(nc[k] == nb[k-r]);
		for(size_t ic = 0u; ic < nc[k]; b += wb[k-r], c += wc[k], ++ic)
			ttt(k+1, r, s, q,  c, nc, wc,   a, na, wa,   b, nb, wb);
	}
	else if(k < r+s+q-1)
	{
		assert(na[k-s] == nb[k-r]);
		for(size_t ia = 0u; ia < na[k-s]; a += wa[k-s], b += wb[k-r], ++ia)
			ttt(k+1, r, s, q,  c, nc, wc,   a, na, wa,   b, nb, wb);
	}
	else
	{
		assert(na[k-s] == nb[k-r]);
		for(size_t ia = 0u; ia < na[k-s]; a += wa[k-s], b += wb[k-r], ++ia)
			*c += *a * *b;
	}
}


/** @brief Computes the tensor-times-matrix product for the contraction mode m > 0
 *
 * Implements C[i1,i2,...,im-1,j,im+1,...,ip] = sum(A[i1,i2,...,im,...,ip] * B[j,im])
 *
 * @note is used in function ttm
 *
 * @param m  zero-based contraction mode with 0<m<p
 * @param r  zero-based recursion level starting with p-1
 * @param c  pointer to the output tensor
 * @param nc pointer to the extents of tensor c
 * @param wc pointer to the strides of tensor c
 * @param a  pointer to the first input tensor
 * @param na pointer to the extents of input tensor a
 * @param wa pointer to the strides of input tensor a
 * @param b  pointer to the second input tensor
 * @param nb pointer to the extents of input tensor b
 * @param wb pointer to the strides of input tensor b
*/

template <class PointerOut, class PointerIn1, class PointerIn2, class SizeType>
void ttm(SizeType const m,  SizeType const r,
         PointerOut c, SizeType const*const nc, SizeType const*const wc,
         PointerIn1 a, SizeType const*const na, SizeType const*const wa,
         PointerIn2 b, SizeType const*const nb, SizeType const*const wb)
{

	if(r == m) {
		ttm(m, r-1, c, nc, wc,    a, na, wa,    b, nb, wb);
	}
	else if(r == 0){
		for(auto i0 = 0ul; i0 < nc[0]; c += wc[0], a += wa[0], ++i0) {
			auto cm = c;
			auto b0 = b;
			for(auto i0 = 0ul; i0 < nc[m]; cm += wc[m], b0 += wb[0], ++i0){
				auto am = a;
				auto b1 = b0;
				for(auto i1 = 0ul; i1 < nb[1]; am += wa[m], b1 += wb[1], ++i1)
					*cm += *am * *b1;
			}
		}
	}

	else{
		for(auto i = 0ul; i < na[r]; c += wc[r], a += wa[r], ++i)
			ttm(m, r-1, c, nc, wc,    a, na, wa,    b, nb, wb);
	}
}

/** @brief Computes the tensor-times-matrix product for the contraction mode m = 0
 *
 * Implements C[j,i2,...,ip] = sum(A[i1,i2,...,ip] * B[j,i1])
 *
 * @note is used in function ttm
 *
 * @param m  zero-based contraction mode with 0<m<p
 * @param r  zero-based recursion level starting with p-1
 * @param c  pointer to the output tensor
 * @param nc pointer to the extents of tensor c
 * @param wc pointer to the strides of tensor c
 * @param a  pointer to the first input tensor
 * @param na pointer to the extents of input tensor a
 * @param wa pointer to the strides of input tensor a
 * @param b  pointer to the second input tensor
 * @param nb pointer to the extents of input tensor b
 * @param wb pointer to the strides of input tensor b
*/
template <class PointerOut, class PointerIn1, class PointerIn2, class SizeType>
void ttm0( SizeType const r,
           PointerOut c, SizeType const*const nc, SizeType const*const wc,
           PointerIn1 a, SizeType const*const na, SizeType const*const wa,
           PointerIn2 b, SizeType const*const nb, SizeType const*const wb)
{

	if(r > 1){
		for(auto i = 0ul; i < na[r]; c += wc[r], a += wa[r], ++i)
			ttm0(r-1, c, nc, wc,    a, na, wa,    b, nb, wb);
	}
	else{
		for(auto i1 = 0ul; i1 < nc[1]; c += wc[1], a += wa[1], ++i1) {
			auto cm = c;
			auto b0 = b;
			// r == m == 0
			for(auto i0 = 0ul; i0 < nc[0]; cm += wc[0], b0 += wb[0], ++i0){

				auto am = a;
				auto b1 = b0;
				for(auto i1 = 0u; i1 < nb[1]; am += wa[0], b1 += wb[1], ++i1){

					*cm += *am * *b1;
				}
			}
		}
	}
}


//////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////


/** @brief Computes the tensor-times-vector product for the contraction mode m > 0
 *
 * Implements C[i1,i2,...,im-1,im+1,...,ip] = sum(A[i1,i2,...,im,...,ip] * b[im])
 *
 * @note is used in function ttv
 *
 * @param m  zero-based contraction mode with 0<m<p
 * @param r  zero-based recursion level starting with p-1 for tensor A
 * @param q  zero-based recursion level starting with p-1 for tensor C
 * @param c  pointer to the output tensor
 * @param nc pointer to the extents of tensor c
 * @param wc pointer to the strides of tensor c
 * @param a  pointer to the first input tensor
 * @param na pointer to the extents of input tensor a
 * @param wa pointer to the strides of input tensor a
 * @param b  pointer to the second input tensor
*/

template <class PointerOut, class PointerIn1, class PointerIn2, class SizeType>
void ttv( SizeType const m, SizeType const r, SizeType const q,
          PointerOut c, SizeType const*const nc, SizeType const*const wc,
          PointerIn1 a, SizeType const*const na, SizeType const*const wa,
          PointerIn2 b)
{

	if(r == m) {
		ttv(m, r-1, q, c, nc, wc,    a, na, wa,    b);
	}
	else if(r == 0){
		for(auto i0 = 0u; i0 < na[0]; c += wc[0], a += wa[0], ++i0) {
			auto c1 = c; auto a1 = a; auto b1 = b;
			for(auto im = 0u; im < na[m]; a1 += wa[m], ++b1, ++im)
				*c1 += *a1 * *b1;
		}
	}
	else{
		for(auto i = 0u; i < na[r]; c += wc[q], a += wa[r], ++i)
			ttv(m, r-1, q-1, c, nc, wc,    a, na, wa,    b);
	}
}


/** @brief Computes the tensor-times-vector product for the contraction mode m = 0
 *
 * Implements C[i2,...,ip] = sum(A[i1,...,ip] * b[i1])
 *
 * @note is used in function ttv
 *
 * @param m  zero-based contraction mode with m=0
 * @param r  zero-based recursion level starting with p-1
 * @param c  pointer to the output tensor
 * @param nc pointer to the extents of tensor c
 * @param wc pointer to the strides of tensor c
 * @param a  pointer to the first input tensor
 * @param na pointer to the extents of input tensor a
 * @param wa pointer to the strides of input tensor a
 * @param b  pointer to the second input tensor
*/
template <class PointerOut, class PointerIn1, class PointerIn2, class SizeType>
void ttv0(SizeType const r,
          PointerOut c, SizeType const*const nc, SizeType const*const wc,
          PointerIn1 a, SizeType const*const na, SizeType const*const wa,
          PointerIn2 b)
{

	if(r > 1){
		for(auto i = 0u; i < na[r]; c += wc[r-1], a += wa[r], ++i)
			ttv0(r-1, c, nc, wc,    a, na, wa,    b);
	}
	else{
		for(auto i1 = 0u; i1 < na[1]; c += wc[0], a += wa[1], ++i1)
		{
			auto c1 = c; auto a1 = a; auto b1 = b;
			for(auto i0 = 0u; i0 < na[0]; a1 += wa[0], ++b1, ++i0)
				*c1 += *a1 * *b1;
		}
	}
}


/** @brief Computes the matrix-times-vector product
 *
 * Implements C[i1] = sum(A[i1,i2] * b[i2]) or C[i2] = sum(A[i1,i2] * b[i1])
 *
 * @note is used in function ttv
 *
 * @param[in]  m  zero-based contraction mode with m=0 or m=1
 * @param[out] c  pointer to the output tensor C
 * @param[in]  nc pointer to the extents of tensor C
 * @param[in]  wc pointer to the strides of tensor C
 * @param[in]  a  pointer to the first input tensor A
 * @param[in]  na pointer to the extents of input tensor A
 * @param[in]  wa pointer to the strides of input tensor A
 * @param[in]  b  pointer to the second input tensor B
*/
template <class PointerOut, class PointerIn1, class PointerIn2, class SizeType>
void mtv(SizeType const m,
         PointerOut c, SizeType const*const   , SizeType const*const wc,
         PointerIn1 a, SizeType const*const na, SizeType const*const wa,
         PointerIn2 b)
{
	// decides whether matrix multiplied with vector or vector multiplied with matrix
	const auto o = (m == 0) ? 1 : 0;

	for(auto io = 0u; io < na[o]; c += wc[o], a += wa[o], ++io) {
		auto c1 = c; auto a1 = a; auto b1 = b;
		for(auto im = 0u; im < na[m]; a1 += wa[m], ++b1, ++im)
			*c1 += *a1 * *b1;
	}
}


/** @brief Computes the matrix-times-matrix product
 *
 * Implements C[i1,i3] = sum(A[i1,i2] * B[i2,i3])
 *
 * @note is used in function ttm
 *
 * @param[out] c  pointer to the output tensor C
 * @param[in]  nc pointer to the extents of tensor C
 * @param[in]  wc pointer to the strides of tensor C
 * @param[in]  a  pointer to the first input tensor A
 * @param[in]  na pointer to the extents of input tensor A
 * @param[in]  wa pointer to the strides of input tensor A
 * @param[in]  b  pointer to the second input tensor B
 * @param[in]  nb pointer to the extents of input tensor B
 * @param[in]  wb pointer to the strides of input tensor B
*/
template <class PointerOut, class PointerIn1, class PointerIn2, class SizeType>
void mtm(PointerOut c, SizeType const*const nc, SizeType const*const wc,
         PointerIn1 a, SizeType const*const na, SizeType const*const wa,
         PointerIn2 b, SizeType const*const nb, SizeType const*const wb)
{

	// C(i,j) = A(i,k) * B(k,j)

	assert(nc[0] == na[0]);
	assert(nc[1] == nb[1]);
	assert(na[1] == nb[0]);

	auto cj = c; auto bj = b;
	for(auto j = 0u; j < nc[1]; cj += wc[1], bj += wb[1], ++j) {

		auto bk = bj; auto ak = a;
		for(auto k = 0u; k < na[1]; ak += wa[1], bk += wb[0], ++k) {

			auto ci = cj; auto ai = ak;
			for(auto i = 0u; i < na[0]; ai += wa[0], ci += wc[0], ++i){
				*ci += *ai * *bk;
			}

		}

	}
}



/** @brief Computes the inner product of two tensors
 *
 * Implements c = sum(A[i1,i2,...,ip] * B[i1,i2,...,ip])
 *
 * @note is used in function inner
 *
 * @param r   zero-based recursion level starting with p-1
 * @param n   pointer to the extents of input or output tensor
 * @param a   pointer to the first input tensor
 * @param wa  pointer to the strides of input tensor a
 * @param b   pointer to the second input tensor
 * @param wb  pointer to the strides of tensor b
 * @param v   previously computed value (start with v = 0).
 * @return    inner product of two tensors.
*/
template <class PointerIn1, class PointerIn2, class value_t, class SizeType>
value_t inner(SizeType const r, SizeType const*const n,
              PointerIn1  a, SizeType const*const wa,
              PointerIn2  b, SizeType const*const wb,
              value_t v)
{
	if(r == 0)
		for(auto i0 = 0u; i0 < n[0]; a += wa[0], b += wb[0], ++i0)
			v += *a * *b;
	else
		for(auto ir = 0u; ir < n[r]; a += wa[r], b += wb[r], ++ir)
			v = inner(r-1, n,   a, wa,    b, wb, v);
	return v;
}


template <class PointerOut, class PointerIn1, class PointerIn2, class SizeType>
void outer_2x2(SizeType const pa,
               PointerOut c, SizeType const*const   , SizeType const*const wc,
               PointerIn1 a, SizeType const*const na, SizeType const*const wa,
               PointerIn2 b, SizeType const*const nb, SizeType const*const wb)
{
	//	assert(rc == 3);
	//	assert(ra == 1);
	//	assert(rb == 1);

	for(auto ib1 = 0u; ib1 < nb[1]; b += wb[1], c += wc[pa+1], ++ib1) {
		auto c2 = c;
		auto b0 = b;
		for(auto ib0 = 0u; ib0 < nb[0]; b0 += wb[0], c2 += wc[pa], ++ib0) {
			const auto b = *b0;
			auto c1 = c2;
			auto a1 = a;
			for(auto ia1 = 0u; ia1 < na[1]; a1 += wa[1], c1 += wc[1], ++ia1) {
				auto a0 = a1;
				auto c0 = c1;
				for(SizeType ia0 = 0u; ia0 < na[0]; a0 += wa[0], c0 += wc[0], ++ia0)
					*c0 = *a0 * b;
			}
		}
	}
}

/** @brief Computes the outer product of two tensors
 *
 * Implements C[i1,...,ip,j1,...,jq] = A[i1,i2,...,ip] * B[j1,j2,...,jq]
 *
 * @note called by outer
 *
 *
 * @param[in]  pa number of dimensions (rank) of the first input tensor A with pa > 0
 *
 * @param[in]  rc recursion level for C that starts with pc-1
 * @param[out] c  pointer to the output tensor
 * @param[in]  nc pointer to the extents of output tensor c
 * @param[in]  wc pointer to the strides of output tensor c
 *
 * @param[in]  ra recursion level for A that starts with pa-1
 * @param[in]  a  pointer to the first input tensor
 * @param[in]  na pointer to the extents of the first input tensor a
 * @param[in]  wa pointer to the strides of the first input tensor a
 *
 * @param[in]  rb recursion level for B that starts with pb-1
 * @param[in]  b  pointer to the second input tensor
 * @param[in]  nb pointer to the extents of the second input tensor b
 * @param[in]  wb pointer to the strides of the second input tensor b
*/
template<class PointerOut, class PointerIn1, class PointerIn2, class SizeType>
void outer(SizeType const pa,
           SizeType const rc, PointerOut c, SizeType const*const nc, SizeType const*const wc,
           SizeType const ra, PointerIn1 a, SizeType const*const na, SizeType const*const wa,
           SizeType const rb, PointerIn2 b, SizeType const*const nb, SizeType const*const wb)
{
	if(rb > 1)
		for(auto ib = 0u; ib < nb[rb]; b += wb[rb], c += wc[rc], ++ib)
			outer(pa, rc-1, c, nc, wc,    ra, a, na, wa,    rb-1, b, nb, wb);
	else if(ra > 1)
		for(auto ia = 0u; ia < na[ra]; a += wa[ra], c += wc[ra], ++ia)
			outer(pa, rc-1, c, nc, wc,   ra-1, a, na, wa,   rb, b, nb, wb);
	else
		outer_2x2(pa, c, nc, wc,   a, na, wa,    b, nb, wb); //assert(ra==1 && rb==1 && rc==3);
}




/** @brief Computes the outer product with permutation tuples
 *
 * Implements C[i1,...,ir,j1,...,js] = sum( A[i1,...,ir] * B[j1,...,js]  )
 *
 * nc[x]         = na[phia[x]] for 1 <= x <= r
 * nc[r+x]       = nb[phib[x]] for 1 <= x <= s
 *
 * @note maybe called by ttt function
 *
 * @param k  zero-based recursion level starting with 0
 * @param r  number of non-contraction indices of A
 * @param s  number of non-contraction indices of B
 * @param phia pointer to the permutation tuple of length r for A
 * @param phib pointer to the permutation tuple of length s for B
 * @param c  pointer to the output tensor C with rank(A)=r+s
 * @param nc pointer to the extents of tensor C
 * @param wc pointer to the strides of tensor C
 * @param a  pointer to the first input tensor with rank(A)=r
 * @param na pointer to the extents of the first input tensor A
 * @param wa pointer to the strides of the first input tensor A
 * @param b  pointer to the second input tensor B with rank(B)=s
 * @param nb pointer to the extents of the second input tensor B
 * @param wb pointer to the strides of the second input tensor B
*/

template <class PointerOut, class PointerIn1, class PointerIn2, class SizeType>
void outer(SizeType const k,
           SizeType const r, SizeType const s,
           SizeType const*const phia, SizeType const*const phib,
           PointerOut c, SizeType const*const nc, SizeType const*const wc,
           PointerIn1 a, SizeType const*const na, SizeType const*const wa,
           PointerIn2 b, SizeType const*const nb, SizeType const*const wb)
{
	if(k < r)
	{
		assert(nc[k] == na[phia[k]-1]);
		for(size_t ic = 0u; ic < nc[k]; a += wa[phia[k]-1], c += wc[k], ++ic)
			outer(k+1, r, s,   phia,phib,  c, nc, wc,   a, na, wa,   b, nb, wb);
	}
	else if(k < r+s-1)
	{
		assert(nc[k] == nb[phib[k-r]-1]);
		for(size_t ic = 0u; ic < nc[k]; b += wb[phib[k-r]-1], c += wc[k], ++ic)
			outer(k+1, r, s, phia, phib,    c, nc, wc,   a, na, wa,   b, nb, wb);
	}
	else
	{
		assert(nc[k] == nb[phib[k-r]-1]);
		for(size_t ic = 0u; ic < nc[k]; b += wb[phib[k-r]-1], c += wc[k], ++ic)
			*c = *a * *b;
	}
}


} // namespace recursive
} // namespace detail
} // namespace ublas
} // namespace numeric
} // namespace boost




//////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////


#include <stdexcept>

namespace boost {
namespace numeric {
namespace ublas {

/** @brief Computes the tensor-times-vector product
 *
 * Implements
 *   C[i1,i2,...,im-1,im+1,...,ip] = sum(A[i1,i2,...,im,...,ip] * b[im]) for m>1 and
 *   C[i2,...,ip]                  = sum(A[i1,...,ip]           * b[i1]) for m=1
 *
 * @note calls detail::ttv, detail::ttv0 or detail::mtv
 *
 * @param[in]  m  contraction mode with 0 < m <= p
 * @param[in]  p  number of dimensions (rank) of the first input tensor with p > 0
 * @param[out] c  pointer to the output tensor with rank p-1
 * @param[in]  nc pointer to the extents of tensor c
 * @param[in]  wc pointer to the strides of tensor c
 * @param[in]  a  pointer to the first input tensor
 * @param[in]  na pointer to the extents of input tensor a
 * @param[in]  wa pointer to the strides of input tensor a
 * @param[in]  b  pointer to the second input tensor
 * @param[in]  nb pointer to the extents of input tensor b
 * @param[in]  wb pointer to the strides of input tensor b
*/
template <class PointerOut, class PointerIn1, class PointerIn2, class SizeType>
void ttv(SizeType const m, SizeType const p,
         PointerOut c,       SizeType const*const nc, SizeType const*const wc,
         const PointerIn1 a, SizeType const*const na, SizeType const*const wa,
         const PointerIn2 b, SizeType const*const nb, SizeType const*const wb)
{
	static_assert( std::is_pointer<PointerOut>::value & std::is_pointer<PointerIn1>::value & std::is_pointer<PointerIn2>::value,
	               "Static error in boost::numeric::ublas::ttv: Argument types for pointers are not pointer types.");

	if( m == 0)
		throw std::length_error("Error in boost::numeric::ublas::ttv: Contraction mode must be greater than zero.");

	if( p < m )
		throw std::length_error("Error in boost::numeric::ublas::ttv: Rank must be greater equal the modus.");

	if( p == 0)
		throw std::length_error("Error in boost::numeric::ublas::ttv: Rank must be greater than zero.");

	if(c == nullptr || a == nullptr || b == nullptr)
		throw std::length_error("Error in boost::numeric::ublas::ttv: Pointers shall not be null pointers.");

	for(auto i = 0u; i < m-1; ++i)
		if(na[i] != nc[i])
			throw std::length_error("Error in boost::numeric::ublas::ttv: Extents (except of dimension mode) of A and C must be equal.");

	for(auto i = m; i < p; ++i)
		if(na[i] != nc[i-1])
			throw std::length_error("Error in boost::numeric::ublas::ttv: Extents (except of dimension mode) of A and C must be equal.");

	const auto max = std::max(nb[0], nb[1]);
	if(  na[m-1] != max)
		throw std::length_error("Error in boost::numeric::ublas::ttv: Extent of dimension mode of A and b must be equal.");


	if((m != 1) && (p > 2))
		detail::recursive::ttv(m-1, p-1, p-2, c, nc, wc,    a, na, wa,   b);
	else if ((m == 1) && (p > 2))
		detail::recursive::ttv0(p-1, c, nc, wc,  a, na, wa,   b);
	else if( p == 2 )
		detail::recursive::mtv(m-1, c, nc, wc,  a, na, wa,   b);
	else /*if( p == 1 )*/{
		auto v = std::remove_pointer_t<std::remove_cv_t<PointerOut>>{};
		*c = detail::recursive::inner(SizeType(0), na, a, wa, b, wb, v);
	}

}



/** @brief Computes the tensor-times-matrix product
 *
 * Implements
 *   C[i1,i2,...,im-1,j,im+1,...,ip] = sum(A[i1,i2,...,im,...,ip] * B[j,im]) for m>1 and
 *   C[j,i2,...,ip]                  = sum(A[i1,i2,...,ip]        * B[j,i1]) for m=1
 *
 * @note calls detail::ttm or detail::ttm0
 *
 * @param[in]  m  contraction mode with 0 < m <= p
 * @param[in]  p  number of dimensions (rank) of the first input tensor with p > 0
 * @param[out] c  pointer to the output tensor with rank p-1
 * @param[in]  nc pointer to the extents of tensor c
 * @param[in]  wc pointer to the strides of tensor c
 * @param[in]  a  pointer to the first input tensor
 * @param[in]  na pointer to the extents of input tensor a
 * @param[in]  wa pointer to the strides of input tensor a
 * @param[in]  b  pointer to the second input tensor
 * @param[in]  nb pointer to the extents of input tensor b
 * @param[in]  wb pointer to the strides of input tensor b
*/

template <class PointerIn1, class PointerIn2, class PointerOut, class SizeType>
void ttm(SizeType const m, SizeType const p,
         PointerOut c,       SizeType const*const nc, SizeType const*const wc,
         const PointerIn1 a, SizeType const*const na, SizeType const*const wa,
         const PointerIn2 b, SizeType const*const nb, SizeType const*const wb)
{

	static_assert( std::is_pointer<PointerOut>::value & std::is_pointer<PointerIn1>::value & std::is_pointer<PointerIn2>::value,
	               "Static error in boost::numeric::ublas::ttm: Argument types for pointers are not pointer types.");

	if( m == 0 )
		throw std::length_error("Error in boost::numeric::ublas::ttm: Contraction mode must be greater than zero.");

	if( p < m )
		throw std::length_error("Error in boost::numeric::ublas::ttm: Rank must be greater equal than the specified mode.");

	if( p == 0)
		throw std::length_error("Error in boost::numeric::ublas::ttm:Rank must be greater than zero.");

	if(c == nullptr || a == nullptr || b == nullptr)
		throw std::length_error("Error in boost::numeric::ublas::ttm: Pointers shall not be null pointers.");

	for(auto i = 0u; i < m-1; ++i)
		if(na[i] != nc[i])
			throw std::length_error("Error in boost::numeric::ublas::ttm: Extents (except of dimension mode) of A and C must be equal.");

	for(auto i = m; i < p; ++i)
		if(na[i] != nc[i])
			throw std::length_error("Error in boost::numeric::ublas::ttm: Extents (except of dimension mode) of A and C must be equal.");

	if(na[m-1] != nb[1])
		throw std::length_error("Error in boost::numeric::ublas::ttm: 2nd Extent of B and M-th Extent of A must be the equal.");

	if(nc[m-1] != nb[0])
		throw std::length_error("Error in boost::numeric::ublas::ttm: 1nd Extent of B and M-th Extent of C must be the equal.");

	if ( m != 1 )
		detail::recursive::ttm (m-1, p-1, c, nc, wc,    a, na, wa,   b, nb, wb);
	else /*if (m == 1 && p >  2)*/
		detail::recursive::ttm0(     p-1, c, nc, wc,    a, na, wa,   b, nb, wb);

}


/** @brief Computes the tensor-times-tensor product
 *
 * Implements C[i1,...,ir,j1,...,js] = sum( A[i1,...,ir+q] * B[j1,...,js+q]  )
 *
 * @note calls detail::recursive::ttt or ttm or ttv or inner or outer
 *
 * nc[x]         = na[phia[x]  ] for 1 <= x <= r
 * nc[r+x]       = nb[phib[x]  ] for 1 <= x <= s
 * na[phia[r+x]] = nb[phib[s+x]] for 1 <= x <= q
 *
 * @param[in]  pa number of dimensions (rank) of the first input tensor a with pa > 0
 * @param[in]  pb number of dimensions (rank) of the second input tensor b with pb > 0
 * @param[in]	 q  number of contraction dimensions with pa >= q and pb >= q and q >= 0
 * @param[in]	 phia pointer to a permutation tuple for the first input tensor a
 * @param[in]	 phib pointer to a permutation tuple for the second input tensor b
 * @param[out] c  pointer to the output tensor with rank p-1
 * @param[in]  nc pointer to the extents of tensor c
 * @param[in]  wc pointer to the strides of tensor c
 * @param[in]  a  pointer to the first input tensor
 * @param[in]  na pointer to the extents of input tensor a
 * @param[in]  wa pointer to the strides of input tensor a
 * @param[in]  b  pointer to the second input tensor
 * @param[in]  nb pointer to the extents of input tensor b
 * @param[in]  wb pointer to the strides of input tensor b
*/

template <class PointerIn1, class PointerIn2, class PointerOut, class SizeType>
void ttt(SizeType const pa, SizeType const pb, SizeType const q,
         SizeType const*const phia, SizeType const*const phib,
         PointerOut c, SizeType const*const nc, SizeType const*const wc,
         PointerIn1 a, SizeType const*const na, SizeType const*const wa,
         PointerIn2 b, SizeType const*const nb, SizeType const*const wb)
{
	static_assert( std::is_pointer<PointerOut>::value & std::is_pointer<PointerIn1>::value & std::is_pointer<PointerIn2>::value,
	               "Static error in boost::numeric::ublas::ttm: Argument types for pointers are not pointer types.");

	if( pa == 0 || pb == 0)
		throw std::length_error("Error in boost::numeric::ublas::ttt: tensor order must be greater zero.");

	if( q > pa && q > pb)
		throw std::length_error("Error in boost::numeric::ublas::ttt: number of contraction must be smaller than or equal to the tensor order.");


	SizeType const r = pa - q;
	SizeType const s = pb - q;

	if(c == nullptr || a == nullptr || b == nullptr)
		throw std::length_error("Error in boost::numeric::ublas::ttm: Pointers shall not be null pointers.");

	for(auto i = 0ul; i < r; ++i)
		if( na[phia[i]-1] != nc[i] )
			throw std::length_error("Error in boost::numeric::ublas::ttt: dimensions of lhs and res tensor not correct.");

	for(auto i = 0ul; i < s; ++i)
		if( nb[phib[i]-1] != nc[r+i] )
			throw std::length_error("Error in boost::numeric::ublas::ttt: dimensions of rhs and res not correct.");

	for(auto i = 0ul; i < q; ++i)
		if( nb[phib[s+i]-1] != na[phia[r+i]-1] )
			throw std::length_error("Error in boost::numeric::ublas::ttt: dimensions of lhs and rhs not correct.");


	if(q == 0ul)
		detail::recursive::outer(SizeType{0},r,s,  phia,phib, c,nc,wc, a,na,wa, b,nb,wb);
	else
		detail::recursive::ttt(SizeType{0},r,s,q,  phia,phib, c,nc,wc, a,na,wa, b,nb,wb);
}



/** @brief Computes the tensor-times-tensor product
 *
 * Implements C[i1,...,ir,j1,...,js] = sum( A[i1,...,ir+q] * B[j1,...,js+q]  )
 *
 * @note calls detail::recursive::ttt or ttm or ttv or inner or outer
 *
 * nc[x]   = na[x  ] for 1 <= x <= r
 * nc[r+x] = nb[x  ] for 1 <= x <= s
 * na[r+x] = nb[s+x] for 1 <= x <= q
 *
 * @param[in]  pa number of dimensions (rank) of the first input tensor a with pa > 0
 * @param[in]  pb number of dimensions (rank) of the second input tensor b with pb > 0
 * @param[in]	 q  number of contraction dimensions with pa >= q and pb >= q and q >= 0
 * @param[out] c  pointer to the output tensor with rank p-1
 * @param[in]  nc pointer to the extents of tensor c
 * @param[in]  wc pointer to the strides of tensor c
 * @param[in]  a  pointer to the first input tensor
 * @param[in]  na pointer to the extents of input tensor a
 * @param[in]  wa pointer to the strides of input tensor a
 * @param[in]  b  pointer to the second input tensor
 * @param[in]  nb pointer to the extents of input tensor b
 * @param[in]  wb pointer to the strides of input tensor b
*/

template <class PointerIn1, class PointerIn2, class PointerOut, class SizeType>
void ttt(SizeType const pa, SizeType const pb, SizeType const q,
         PointerOut c, SizeType const*const nc, SizeType const*const wc,
         PointerIn1 a, SizeType const*const na, SizeType const*const wa,
         PointerIn2 b, SizeType const*const nb, SizeType const*const wb)
{
	static_assert( std::is_pointer<PointerOut>::value & std::is_pointer<PointerIn1>::value & std::is_pointer<PointerIn2>::value,
	               "Static error in boost::numeric::ublas::ttm: Argument types for pointers are not pointer types.");

	if( pa == 0 || pb == 0)
		throw std::length_error("Error in boost::numeric::ublas::ttt: tensor order must be greater zero.");

	if( q > pa && q > pb)
		throw std::length_error("Error in boost::numeric::ublas::ttt: number of contraction must be smaller than or equal to the tensor order.");


	SizeType const r  = pa - q;
	SizeType const s  = pb - q;
	SizeType const pc = r+s;

	if(c == nullptr || a == nullptr || b == nullptr)
		throw std::length_error("Error in boost::numeric::ublas::ttm: Pointers shall not be null pointers.");

	for(auto i = 0ul; i < r; ++i)
		if( na[i] != nc[i] )
			throw std::length_error("Error in boost::numeric::ublas::ttt: dimensions of lhs and res tensor not correct.");

	for(auto i = 0ul; i < s; ++i)
		if( nb[i] != nc[r+i] )
			throw std::length_error("Error in boost::numeric::ublas::ttt: dimensions of rhs and res not correct.");

	for(auto i = 0ul; i < q; ++i)
		if( nb[s+i] != na[r+i] )
			throw std::length_error("Error in boost::numeric::ublas::ttt: dimensions of lhs and rhs not correct.");

	using value_type = std::decay_t<decltype(*c)>;



	if(q == 0ul)
		detail::recursive::outer(pa, pc-1, c,nc,wc, pa-1, a,na,wa, pb-1, b,nb,wb);
	else if(r == 0ul && s == 0ul)
		*c = detail::recursive::inner(q-1, na, a,wa,  b,wb, value_type(0) );
	else
		detail::recursive::ttt(SizeType{0},r,s,q, c,nc,wc, a,na,wa, b,nb,wb);
}


/** @brief Computes the inner product of two tensors
 *
 * Implements c = sum(A[i1,i2,...,ip] * B[i1,i2,...,ip])
 *
 * @note calls detail::inner
 *
 * @param[in] p  number of dimensions (rank) of the first input tensor with p > 0
 * @param[in] n  pointer to the extents of input or output tensor
 * @param[in] a  pointer to the first input tensor
 * @param[in] wa pointer to the strides of input tensor a
 * @param[in] b  pointer to the second input tensor
 * @param[in] wb pointer to the strides of input tensor b
 * @param[in] v  inital value
 *
 * @return   inner product of two tensors.
*/
template <class PointerIn1, class PointerIn2, class value_t, class SizeType>
auto inner(const SizeType    p, SizeType const*const n,
           const PointerIn1  a, SizeType const*const wa,
           const PointerIn2  b, SizeType const*const wb,
           value_t v)
{
	static_assert( std::is_pointer<PointerIn1>::value && std::is_pointer<PointerIn2>::value,
	               "Static error in boost::numeric::ublas::inner: Argument types for pointers must be pointer types.");
	if(p<2)
		throw std::length_error("Error in boost::numeric::ublas::inner: Rank must be greater than zero.");
	if(a == nullptr || b == nullptr)
		throw std::length_error("Error in boost::numeric::ublas::inner: Pointers shall not be null pointers.");

	return detail::recursive::inner(p-1, n, a, wa, b, wb, v);

}


/** @brief Computes the outer product of two tensors
 *
 * Implements C[i1,...,ip,j1,...,jq] = A[i1,i2,...,ip] * B[j1,j2,...,jq]
 *
 * @note calls detail::outer
 *
 * @param[out] c  pointer to the output tensor
 * @param[in]  pc number of dimensions (rank) of the output tensor c with pc > 0
 * @param[in]  nc pointer to the extents of output tensor c
 * @param[in]  wc pointer to the strides of output tensor c
 * @param[in]  a  pointer to the first input tensor
 * @param[in]  pa number of dimensions (rank) of the first input tensor a with pa > 0
 * @param[in]  na pointer to the extents of the first input tensor a
 * @param[in]  wa pointer to the strides of the first input tensor a
 * @param[in]  b  pointer to the second input tensor
 * @param[in]  pb number of dimensions (rank) of the second input tensor b with pb > 0
 * @param[in]  nb pointer to the extents of the second input tensor b
 * @param[in]  wb pointer to the strides of the second input tensor b
*/
template <class PointerOut, class PointerIn1, class PointerIn2, class SizeType>
void outer(PointerOut c,       SizeType const pc, SizeType const*const nc, SizeType const*const wc,
           const PointerIn1 a, SizeType const pa, SizeType const*const na, SizeType const*const wa,
           const PointerIn2 b, SizeType const pb, SizeType const*const nb, SizeType const*const wb)
{
	static_assert( std::is_pointer<PointerIn1>::value & std::is_pointer<PointerIn2>::value & std::is_pointer<PointerOut>::value,
	               "Static error in boost::numeric::ublas::outer: argument types for pointers must be pointer types.");
	if(pa < 2u || pb < 2u)
		throw std::length_error("Error in boost::numeric::ublas::outer: number of extents of lhs and rhs tensor must be equal or greater than two.");
	if((pa + pb) != pc)
		throw std::length_error("Error in boost::numeric::ublas::outer: number of extents of lhs plus rhs tensor must be equal to the number of extents of C.");
	if(a == nullptr || b == nullptr || c == nullptr)
		throw std::length_error("Error in boost::numeric::ublas::outer: pointers shall not be null pointers.");

	detail::recursive::outer(pa, pc-1, c, nc, wc,   pa-1, a, na, wa,   pb-1, b, nb, wb);

}




}
}
}

#endif