实现 Eigen 的块或 Boost 的项目方法

Question

我正在努力提高我的 C++ 技能——特别是模板的使用。我正在创建一个过滤数学库，它通过使用模板兼容不同的向量/矩阵类，包括 Boost 和 Eigen。到目前为止，我对该库非常满意，但在扩展功能时遇到了问题。

Eigen 和 Boost 的不同之处在于，后者没有乘法运算符（请参阅 1），因此这为我的实现提供了一个选择。我决定确保模板库中使用的矩阵/向量方法和运算符使用为 Eigen 定义的那些，以便至少有一个库可以很好地与数学库配合使用。在我的库中，我会使用 * 来乘以不能用于 boost 的矩阵。对于 boost，我为矩阵和向量创建了一个包装类，以允许通过自定义运算符使用 *。矩阵包装器示例如下：

class matrixBoost{
        private:
                boost::numeric::ublas::matrix<double> value;

        public:
                matrixBoost(){
                }

                matrixBoost(boost::numeric::ublas::matrix<double> value):value(value){}

                static matrixBoost Identity(int dimension,int dimension2){
                        return matrixBoost(boost::numeric::ublas::identity_matrix<double>(dimension,dimension2));
                }

                matrixBoost transpose(){
                        matrixBoost t(boost::numeric::ublas::trans(value));
                        return t;
                }

                boost::numeric::ublas::matrix<double> getSystemValue() const{
                        return value;
                }
};

//example operator - one of many
matrixBoost operator*(const matrixBoost &lhs, const matrixBoost &rhs){
        matrixBoost res(boost::numeric::ublas::prod(lhs.getSystemValue(),rhs.getSystemValue()));
        return res;
}

我现在正在尝试将 Eigen 的 block 添加到 boost 的包装类中，该类允许在 Eigen 的库中进行以下行为：

Eigen::MatrixXd m(2,2);m<<1,2,3,4;
Eigen::MatrixXd n = m.block(0, 0, 1, 2) + m.block(0,0,1,2);
m.block(0,0,1,2) = n;

现在与m 相等

2 4
3 4

我的第一个问题是有人可以链接或展示如何编码块函数的示例（甚至不考虑包装类）。我试过用谷歌搜索，但要么谷歌变得更糟，要么我没有使用正确的关键词——结果挤满了 operator[] 结果。我很难理解。我尝试在 Eigen 库源中搜索，并想象代码必须位于 here 但模板语言对我来说有点难以轻松解析。

我知道boost也有类似here所示的概念：

project(A, r1, r2) = C;    // assign the submatrix of A specified by the two index ranges r1 and r2
C = project(A, r1, r2);    // the submatrix of A specified by the two index ranges r1 and r2

我不确定是否可以向我的包装类添加块函数。现在我有以下内容：

matrixBoost block(int index1,int index2, int length1, int length2){
    boost::numeric::ublas::range r1(i1,l1);
    boost::numeric::ublas::range r2(i2,l2);
    return matrixBoost(boost::numeric::ublas::project(value,r1,r2));
}

这适用于方法位于等号右侧的情况。但是，对于如何继续允许在左侧调用该方法，我遇到了困难。也许我应该利用boost的项目方法。任何帮助将非常感激！

谢谢。

Answer 1

于是我找到了解决boost类包装问题的方法，但是解决方法有点乱。上述问题中推导出的block 函数必须返回一个引用。如果它在等号的右侧，那么我们需要返回包裹在matrixBoost 中的子矩阵，以用于可能在表达式中的其他 */+ 乘法。但是，如果在等号的左侧调用它呢？

我的解决方案是使用布尔标志 (subMatrixCopy)，将其设置为 true，并通过使用 std::swap 将完整矩阵 value 备份到 valBack 并返回子矩阵。此修改将允许在右侧进行正确的表达式评估。对于左侧，一旦调用了= 运算符，布尔标志就会确保将右侧正确分配给valBack 中备份矩阵值的指定块。

#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/qvm/mat_operations.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/numeric/ublas/triangular.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <lapacke.h>

class matrixBoost{
   private:
        boost::numeric::ublas::matrix<double> value;
        boost::numeric::ublas::range r1;
        boost::numeric::ublas::range r2;
        boost::numeric::ublas::matrix<double>valBack;
        bool subMatrixCopy = false;

        void setSystemValue(boost::numeric::ublas::matrix<double> v){
            value = v;
        }

   public:
        matrixBoost(){}

        matrixBoost(int rowCount, int colCount){
            value = boost::numeric::ublas::zero_matrix<double>(rowCount,colCount);
        }

        matrixBoost(boost::numeric::ublas::matrix<double> value):value(value){}

        static matrixBoost Identity(int dimension,int dimension2){
           return matrixBoost(boost::numeric::ublas::identity_matrix<double>(dimension,dimension2));
         }

         matrixBoost transpose(){
              matrixBoost t(boost::numeric::ublas::trans(value));
              return t;
          }

          boost::numeric::ublas::matrix<double> getSystemValue() const{
              return value;
          }

         matrixBoost & block(int i1, int i2, int l1, int l2){
             if(subMatrixCopy){
                 std::swap(valBack,value);
              }
              subMatrixCopy = true;
              r1 = boost::numeric::ublas::range(i1,i1+l1);
              r2 = boost::numeric::ublas::range(i2,i2+l2);
              valBack = boost::numeric::ublas::project(value,r1,r2);
              std::swap(valBack,value);
              return *this;
        }

        matrixBoost &operator=( const matrixBoost other){
            if(subMatrixCopy){
                subMatrixCopy = false;
                boost::numeric::ublas::project(valBack,r1,r2) = other.getSystemValue();
                std::swap(valBack,value);
            }else{
                value = other.getSystemValue();
            }
            return *this;//better be no chaining of equal signs
        }
};

这是我用来测试的实际代码示例。有用！但是，我相信该解决方案会阻止执行某些操作，例如链接等号。

int main(){
    boost::numeric::ublas::matrix<double> value(3,3);
    value(0,0) = 1;
    value(1,1) = 2;
    value(2,2) = 3;
    matrixBoost valueWrap(value);

    boost::numeric::ublas::vector<double> v(3);
    v(0) = 1;
    v(1) = 1;
    v(2) = 1;
    cout<<vectorBoost::solve(valueWrap,vectorBoost(v))<<endl;

    cout<<valueWrap<<endl;

    matrixBoost newMatrix = valueWrap.block(1,1,2,2);
    cout<<newMatrix<<endl;
    cout<<valueWrap.block(1,1,2,2)<<endl;

    valueWrap.block(1,1,2,2) = 
        valueWrap.block(1,1,2,2)*newMatrix + newMatrix;
    cout<<valueWrap<<endl;
}