实现Eigen的block或者Boost的project方法

Question

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

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

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 的包装器 class 中，它允许在 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

我的第一个问题是有人可以 link 或展示如何编写块函数的示例（甚至不考虑包装器 class）。我曾尝试使用谷歌搜索，但要么 Google 变得更糟，要么我没有使用正确的关键字 - 结果中充斥着 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

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

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包裹问题的方法class，但是这个解决方法有点乱。上述问题中推导的 block 函数必须 return 引用。如果它在等号的右侧，那么我们需要 return 包裹在 matrixBoost 中的子矩阵用于表达式中可能存在的其他 */+ 乘法。但是，如果在等号的左边调用呢？

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

#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;
}

实现Eigen的block或者Boost的project方法

Implementing Eigen's block or Boost's project method

c++

boost

eigen