1

I am working on improving my C++ skills - specifically the use of templates. I am creating a filtering mathematical library that through the use of templates is compatible for different vector/matrix classes including Boost and Eigen. So far I am very happy with the library, but am encountering issues as I am expanding functionality.

Eigen and Boost are different in that for example the latter does not have multiply operators (see 1) so this presented a choice in my implementation. I decided to make sure that the matrix/vector methods and operators used within the template library use those that are defined for Eigen so that at least one library works really well with the mathematical library. In my library I would use * to multiply matrices which is unavailable for boost. For boost I created a wrapper class for both matrix and vector to allow the use of * through custom operators. Example of wrapper for matrix is below:

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

I am now trying to add Eigen's block to boost's wrapper class that allows for the following behaviour in Eigen's library:

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;

with m now equalling

2 4
3 4

My first question is can someone please link or show examples of how the block function would be coded (without considering a wrapper class even). I have tried googling but either Google has gotten worse or I am not using the right key words - the results are swarmed with operator[] results. I am having difficulty understanding. I tried searching within the Eigen library source and imagine the code must be located here but the template language is a little difficult for me to parse easily.

I know that boost also has a similar concept as shown in 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

I am not sure if adding a block function to my wrapper class is even possible. For now I have the following:

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

This will work for when the method is on the right hand side of the the equal sign. However, I am having difficulty as to how to proceed to allowing for the method being called on the left hand side. Perhaps I should make use of boost's project method. Any help would be much appreciated!

Thank you.

stantheman
  • 195
  • 8

1 Answers1

0

So I found a solution for the problem of wrapping the boost class, but the solution is a little messy. The block function as deduced in the question above must return a reference. If it is on the right hand side of the equal sign then we need to return the submatrix wrapped within matrixBoost for the other */+ multiplications that might be in the expression. However, what if this is called on the left hand side of the equal sign?

My solution was to use a boolean flag (subMatrixCopy) where I would set it to true as well as backup the full matrix value in valBack through the use of std::swap and return with the submatrix. This modification would allow for proper expression evaluation on the right hand side. For the left hand side, once the = operator is called then the boolean flag makes sure that the right hand side is properly assigned to the specified block of the backed up matrix value in 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
        }
};

and here is an example of the code in action that I used to test. It works! However, I believe the solution prevents doing certain things such as chaining 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;
}
stantheman
  • 195
  • 8