特征稀疏矩阵与完全自定义类型相乘

Question

我已经说服Eigen做这个了

class A{...};  class B{...}; class Product{...};
auto operator*(const A&,const B&){return Product{..}} 
..operator+..
//"SparseMatrix<A>*SparseVector<B> -> SparseVector<Product>"

SparseMatrix<A> mymat_A; //...
SparseVector<B> myvec_B; //...
SparseVector<Product> result= mymat_A*myvec_B;

然而，我真正想要的是不同的组件产品和输出累加器类型，例如

class A{}class B{} class Product{} class Accumulator;
auto operator*(const A&,const B&)->Product{..}
auto operator+=(Accumulator& lhs,const Product& p){lhs+=p;return lhs;}
SparseVector<Accumulator> result = mymat_A*myvec_B

这可能吗？我需要在这里设置A、B、Prod的特征类型特征，才能走到这一步，并设置这些模板<gt；struct Eigen:：ScalarBinaryOpTraits<A、 B，特征：：internal:：scalar_product_op<A、 B>gt；{typedef产品退货类型；}；

我曾以为我可以用"；。。总和<产品，产品&gt-&gt；收集器，为Prod+Prod->Acc，但这不起作用。

用例=精度与小产品混合，但许多产品需要更大的累加器，并进行非算术操作(信息在图中流动(，刚好适合稀疏矩阵mul所需的存储和遍历模式优化。

如果这不可能，也许有人可以推荐一个可以做到这一点的替代稀疏矩阵库。

如果我屈服于NIH，只推出我自己的产品，我的设计如下-只需使用decltype从运算符重载中推断出单独的求和/乘积类型的存在(提供具有求和重载的自定义乘积类型，从而产生所需的累加器(。但它将在内部使用Acc+=Prod(在声明decltype(Prod+Prod(累加器初始化为"0"之后(

template<typename A,typename B>
SparseVector<decltype(A()*B()+A()*B())>  operator*(const SparseMatrix<A>& ,const SparseVector<B>&){...}
// supply appropriate operator+, +=...

当前代码(如果Eigen在路径中，则使用clang编译的单个源文件(

#pragma once
#include <Eigen/Sparse>
// setup for MatElem*VecElem->Prod
// goal is to add Prod+Prod->Acc,  Acc+=Prod,  Acc(Prod) Acc=0
class MatElem {public:MatElem(){} MatElem(int x){}};
class VecElem {public:VecElem(){} VecElem(int x){}};
class Prod {public:Prod(){} Prod(int x){}};
class Acc {public:Acc(){} Acc(const Prod&){printf("Acc(Prod)");} Acc(Prod&){printf("Acc(Prod)");} Acc(int x){}};
auto operator*(const MatElem& a,const VecElem& b){printf("MatElem*VecElemn");return Prod{};}
auto operator+(const Prod& a,const Prod& b){printf("Prod+Prodn");return Prod{};}
auto operator+=(Prod& a,const Prod& b){printf("Prod+=Prodn");return a;}
auto operator+=(Acc& a,const Prod& b){printf("Acc+=Prodn");return a;}
template<>
class Eigen::NumTraits<MatElem> {
public:
typedef MatElem Real;
typedef MatElem NonInteger; 
typedef MatElem Literal;
typedef MatElem Nested;
enum {
IsInteger=0,
IsSigned=1,
RequireInitialization=1,
IsComplex=0,
ReadCost=1,
AddCost=1,
MulCost=1

};
auto static epsilon(){return MatElem();}
auto static dummy_precision(){return MatElem();}
auto static highest(){return MatElem();}
auto static lowest(){return MatElem();}
auto static digist10(){return 5;}
};
template<>
class Eigen::NumTraits<VecElem> {
public:
typedef VecElem Real;
typedef VecElem NonInteger; 
typedef VecElem Literal;
typedef VecElem Nested;
enum {
IsInteger=0,
IsSigned=1,
RequireInitialization=1,
IsComplex=0,
ReadCost=1,
AddCost=1,
MulCost=1

};
auto static epsilon(){return VecElem{};}
auto static dummy_precision(){return VecElem{};}
auto static highest(){return VecElem{};}
auto static lowest(){return VecElem{};}
auto static digist10(){return 5;}
};
template<>
class Eigen::NumTraits<Prod> {
public:
typedef Prod Real;
typedef Prod NonInteger;    
typedef Prod Literal;
typedef Prod Nested;
enum {
IsInteger=0,
IsSigned=1,
RequireInitialization=1,
IsComplex=0,
ReadCost=1,
AddCost=1,
MulCost=1

};
auto static epsilon(){return Prod{};}
auto static dummy_precision(){return Prod{};}
auto static highest(){return Prod{};}
auto static lowest(){return Prod{};}
auto static digist10(){return 5;}
};
template<>
class Eigen::NumTraits<Acc> {
public:
typedef Acc Real;
typedef Acc NonInteger; 
typedef Acc Literal;
typedef Acc Nested;
enum {
IsInteger=0,
IsSigned=1,
RequireInitialization=1,
IsComplex=0,
ReadCost=1,
AddCost=1,
MulCost=1

};
auto static epsilon(){return Acc{};}
auto static dummy_precision(){return Acc{};}
auto static highest(){return Acc{};}
auto static lowest(){return Acc{};}
auto static digist10(){return 5;}
};
template<>
struct Eigen::ScalarBinaryOpTraits<MatElem,VecElem,Eigen::internal::scalar_product_op<MatElem, VecElem> >{
typedef Prod ReturnType;
};
template<>
struct Eigen::ScalarBinaryOpTraits<Prod,Prod,Eigen::internal::scalar_sum_op<Prod, Prod> >{
typedef Prod ReturnType;
};

void eigen_experiment() {
Eigen::SparseMatrix<MatElem> mymat(3,3);
mymat.insert(0,0)=MatElem{};
mymat.insert(0,1)=MatElem{};
mymat.insert(1,0)=MatElem{};
mymat.insert(1,1)=MatElem{};
Eigen::SparseVector<VecElem> myvec(3);
myvec.insert(0)=VecElem{};
myvec.insert(1)=VecElem{};
// Can't seem to do this with "Acc", even if supplying appropriate OpTraits etc above.
Eigen::SparseVector<Prod> tmp=mymat*myvec;

for (int k=0; k<mymat.outerSize(); ++k){
for (decltype(mymat)::InnerIterator v(mymat,k); v;++v){
printf("%d %dn",v.row(),v.col());
}
}
for (decltype(tmp)::InnerIterator v(tmp); v;++v){
printf("%dn",v.index());
}
}
int main(int argc,const char**){
eigen_experiment();
return 0;
}

Answer 1

感谢@chtz，只是将两个运算符特性都设置为"累加器"；(所需输出(，而实际操作员过载则返回临时"输出"；产品"；类型，它只是工作。我没有想过尝试它，但它似乎只是使用特性来分配输出，当然它不会影响实际的运算符。调试打印似乎证实了它在做我想要的事情。