Hi all,

I'd like to create a dynamic matrix storage class that stores a matrix either in full or in compressed form. It will mostly be used for products of the forms

newMatrix = leftMatrix * compressedMatrix * rightMatrix
newMatrix = leftMatrix * compressedMatrix
newMatrix = compressedMatrix * rightMatrix

Here's the constructor I have so far:

Code: Select all

class CompressedMatrixStorage
{
   CompressedMatrixStorage(const MatrixXd &matrix, double epsilon)
   {
      JacobiSVD< MatrixXd > jacobiSVD(matrix, ComputeThinU | ComputeThinV);
      const VectorXd &sigma = jacobiSVD.singularValues();
      GlobalUnsignedIndex sigmaSize = sigma.size();
      GlobalUnsignedIndex remainingSingularValues = 0;
      for (GlobalUnsignedIndex i = 0; i < sigmaSize; i ++)
         if (sigma(i) > epsilon)
            remainingSingularValues ++;
         else
            break;
      
      if (remainingSingularValues * 2 < sigmaSize)
      {
         _sigma = sigma.segment(0, remainingSingularValues);
         _U = jacobiSVD.matrixU().block(0, 0, matrix.rows(), remainingSingularValues);
         _V = jacobiSVD.matrixV().block(0, 0, matrix.rows(), remainingSingularValues);
      }
      else
      {
         _fullMatrix = matrix;
      }
   }
   
protected:
   MatrixXd _fullMatrix;
   MatrixXd _U;
   MatrixXd _V;
   VectorXd _sigma;
}

but I'm not sure which would be the best way to implement the matrix products, especially when I wish to keep using the * operator. Could you point me in the right direction?

Best,

MrMage

The simplest way is to forget about expression templates and simply add operator* members for products with a template MatrixBase<Derived> object that will return a MatrixXd object by value. The implementation of operator* will have to perform a "if" do decide between the two variants.

ggael wrote:The simplest way is to forget about expression templates and simply add operator* members for products with a template MatrixBase<Derived> object that will return a MatrixXd object by value. The implementation of operator* will have to perform a "if" do decide between the two variants.

Thanks, that's the way I did it.