Hi all,

I 'm a new user on this forum, so let me start by saying thanks for the developing of the library and also for documentation and support. It's really very helpful.

I'm using the the Eigen 3.2.0 for my thesis and I 've stuck on a problem called "reverse communication". In Krylov methods we only need to multiply a matrix by a vector. In case that we don't have the matrix implicitly (that's my case) we can provide a code that returns the result of the multiplication and that's all.

So I created a wrapper class that provides an operator*() for the A matrix and I would use it for the GMRES algorithm (currently unsupported) . In a similar thread https://forum.kde.org/viewtopic.php?f=74&t=110556&p=261852&hilit=bicgstab#p261852
I found that I have to create rows() and cols() methods but it seems that there more much to do.

After the compiler complained about scalar and index types, I 've added templates :

Code: Select all

#ifndef CUSTOMMAT_HPP_
#define CUSTOMMAT_HPP_

#include <iostream>
#include <fstream>
#include <eigen3/Eigen/Sparse>
#include <string>
#include <stdexcept>

template<typename _Scalar, int _Options, typename _Index>
class CustomMat{

private:
      Eigen::SparseMatrix<double> A;

public:
      typedef _Scalar Scalar;
      typedef _Index Index;
      CustomMat(std::string path);
      int rows();
      int cols();
      Eigen::VectorXd operator*(Eigen::VectorXd &v);
      void print();
};


template<typename _Scalar, int _Options, typename _Index>
CustomMat<_Scalar, _Options, _Index>::CustomMat(std::string path)
   :A(){

      unsigned n,nnz;
      std::ifstream fs(path.c_str());
      if (!fs)
         throw std::runtime_error(std::string("Could not open file"));
      fs >> n;
      fs >> nnz;
      A.resize(n,n);

      typedef Eigen::Triplet<double> T;
      std::vector<T> tripletList;
      tripletList.reserve(nnz);
      for (unsigned z = 0; z < nnz; z++) {
         int i,j;
         double data;
         fs >> i;
         i--;
         fs >> j;
         j--;
         fs >> data;
         tripletList.push_back(T(i,j,data));
      }
      A.setFromTriplets(tripletList.begin(), tripletList.end());
      fs.close();
}

template<typename _Scalar, int _Options, typename _Index>
int CustomMat<_Scalar, _Options, _Index>::rows(){
   return A.rows();
}

template<typename _Scalar, int _Options, typename _Index>
int CustomMat<_Scalar, _Options, _Index>::cols(){
   return A.cols();
}

template<typename _Scalar, int _Options, typename _Index>
void CustomMat<_Scalar, _Options, _Index>::print(){
   std::cout<< A << std::endl;
}

template<typename _Scalar, int _Options, typename _Index>
Eigen::VectorXd CustomMat<_Scalar, _Options, _Index>::operator*(Eigen::VectorXd &v){
   return A*v;
}

int main() {
   CustomMat<int,0,int> A("Desktop/data.txt");
   Eigen::MatrixXd *bPtr;
   bPtr = initiator("/Desktop/vector.txt");

   std::cout<<"Matrices read"<<std::endl;
   A.print();
   Eigen::VectorXd x(bPtr->rows());

   Eigen::GMRES < CustomMat<int,0,int> > solver(A);
   solver.setTolerance(1e-6);
   x = solver.solve(*bPtr);

   std::cout << "#iterations:     " << solver.iterations() << std::endl;
   std::cout << "estimated error: " << solver.error()      << std::endl;

   printToFile("Desktop/sol-c.txt",x);

   return 0;
}


This is a test class and the matrix is explicitly declared, but in my thesis, I don't have that matrix. So this class is just for testing reasons. Now the compiler complains:

Code: Select all

/usr/include/eigen3/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h:28:43: error: no type named ‘RealScalar’ in ‘Eigen::IterativeSolverBase<Eigen::GMRES<CustomMat<int, 0, int> > >::MatrixType {aka class CustomMat<int, 0, int>}’
   typedef typename MatrixType::RealScalar RealScalar;

/usr/include/eigen3/unsupported/Eigen/src/IterativeSolvers/GMRES.h:267:15: error: no members matching ‘Eigen::GMRES<CustomMat<int, 0, int> >::Base {aka Eigen::IterativeSolverBase<Eigen::GMRES<CustomMat<int, 0, int> > >}::m_error’ in ‘Eigen::GMRES<CustomMat<int, 0, int> >::Base {aka class Eigen::IterativeSolverBase<Eigen::GMRES<CustomMat<int, 0, int> > >}’
   using Base::m_error;

and some more


Is there any easier way to implement this?