//
|
// 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
|
