具有完全自定义类型的特征稀疏矩阵乘法

如何解决具有完全自定义类型的特征稀疏矩阵乘法

我已经说服 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 设置本征的类型特征才能做到这一点，并设置这些模板 struct Eigen::ScalarBinaryOpTraits >{ typedef 产品返回类型； };

我曾想我可以使用“..sum.. -> Accumulator，为 Prod+Prod->Acc 提供重载来暗示存在一个单独的产品，但这不起作用。

用例=混合精度与小产品，但许多产品需要更大的累加器，并且做非算术的事情（信息在图形间流动）恰好适合稀疏矩阵mul所需的存储和遍历模式优化。

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

如果我屈服于 NIH 并且只是推出我自己的，我的设计将如下 - 只需使用 decltype 从运算符重载推断单独的 sum/product 类型的存在（提供具有 summation 重载产生的自定义产品类型）所需的累加器）。但它会在内部使用 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*VecElem\n");return Prod{};}
auto operator+(const Prod& a,const Prod& b){printf("Prod+Prod\n");return Prod{};}
auto operator+=(Prod& a,const Prod& b){printf("Prod+=Prod\n");return a;}
auto operator+=(Acc& a,const Prod& b){printf("Acc+=Prod\n");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,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,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,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 %d\n",v.row(),v.col());
        }
    }

    for (decltype(tmp)::InnerIterator v(tmp); v;++v){
        printf("%d\n",v.index());
    }
}

int main(int argc,const char**){
    eigen_experiment();
    return 0;
}

解决方法

感谢@chtz，只需将两个运算符特征设置为“累加器”（所需的输出），而实际的运算符重载返回临时“产品”类型，它就可以工作。我没想过要尝试它，但它似乎只是使用 trait 来分配输出，当然它不会影响实际的操作符。调试打印似乎确认它在做我想做的事情。