I'm trying to test the library and I've some files that I use to populate six square sparse matrices that I try to solve with BiCGStab and ILUT preconditioner, their sizes are 1891, 5781, 10671, 16561, 23451 and 46902.

I generate a X0 Dense Vector that contains all elements from 0 to size-1, i.e. X0_1891 = {0,1,....,1890}.

Then I find

Code: Select all

B = SparseA * X0

Finally using different kinds of solvers I should get the solutions and, in case of iterative solvers, the relative error.

I've done this program (I'm not confident with C++ so it might not be so optimized), but just the first matrix seems to converge and the time it takes to solve the problem is quite high, about 12-13 seconds, when it takes only few seconds with another Java library that I expected to be slower.

Code: Select all

#include <iostream>
#include <fstream>
#include <vector>
#include <string>
#include <sstream>
#include <chrono>
#include "Eigen/SparseCore"
#include "Eigen/Sparse"
#include "Eigen/Dense"
#include "Eigen/IterativeLinearSolvers"

using namespace Eigen;
using namespace std;

typedef Triplet<double> T;

int dim[6] = { 1891, 5781, 10671, 16561, 23451, 46902 };
string path_ns = "Matrices/square-";

string _dat = ".dat";

VectorXd Xe, B, X;

//SparseLU<SparseMatrix<double> > solver2;
void loop(SparseMatrix<double> &mat) {
   vector<T> tripletList;
   int i = 0;
   string line;
   bool first;
   bool simmetric = false;
   string sz, nz, half;
   int size, rows, cols;
   double value;
   for (i = 0; i < 6; i++) {
      cout << "Matrix " << (i + 1) << ':' << endl;
      cout << "============" << endl;
      stringstream msg;
      first = true;
      msg << path_s << dim[i] << _dat;
      ifstream iFile(msg.str());
      while (getline(iFile, line)) {
         istringstream iss(line);
         if (first) {
            first = false;
            //cout << line << endl;
            getline(iss, sz, ' ');
            getline(iss, nz, ' ');
            getline(iss, half, '\n');

            //stringstream sstream;
            //sstream<<sz<<' '<<nz<<' '<<half<<'\n';
            //cout<<sstream.str();

            size = stoi(sz);
            mat.resize(size, size);
            Xe.resize(size);
            B.resize(size);
            tripletList.reserve(stoi(nz));

            //cout<<mat.rows()<<' '<<'x'<<' '<<mat.cols()<<endl;

            if (half.compare("half"))
               simmetric = true;
         } else {
            getline(iss, sz, ' ');
            getline(iss, nz, ' ');
            getline(iss, half, '\n');
            rows = stoi(sz);
            cols = stoi(nz);
            value = stod(half);
            if (simmetric && rows != cols) {
               tripletList.push_back(T((cols - 1), (rows - 1), value));
            }
            tripletList.push_back(T((rows - 1), (cols - 1), value));
         }
      }
      iFile.close();
      auto t_start = chrono::high_resolution_clock::now();
      mat.setFromTriplets(tripletList.begin(), tripletList.end());
      auto t_end = chrono::high_resolution_clock::now();
      cout << "[" << size << "] " << "Time to populate "
            << chrono::duration_cast<chrono::milliseconds>(t_end - t_start).count()
            << " ms.\n";
      //cout<<endl;
      for (int k = 0; k < Xe.size(); k++) {
         Xe(k) = k;
      }
      t_start = chrono::high_resolution_clock::now();
      B = mat * Xe;
      t_end = chrono::high_resolution_clock::now();
      cout << "[" << size << "] " << "Time to find B by multiplication "
            << chrono::duration_cast<chrono::milliseconds>(t_end - t_start).count()
            << " ms.\n";
      t_start = chrono::high_resolution_clock::now();
      if (!simmetric) {
         BiCGSTAB<SparseMatrix<double>, IncompleteLUT<double> > solver;
         solver.setMaxIterations(100000);
         solver.setTolerance(10.0E-6);
         solver.preconditioner().setDroptol(.001);
         solver.compute(mat);
         solver.analyzePattern(mat);
         solver.factorize(mat);

         if (solver.info() != Success) {
            cout << "[" << size << "] Solver compute fail. ("
                  << solver.info() << ")\n";
            //continue;
         }
         cout<<"["<<size<<"] Start\n";

         X = solver.solve(B);
         //X = solver.solveWithGuess(B, Xe);
         if (solver.info() != Success) {
            cout << "[" << size << "] Solver solve fail. (" << solver.info()
                  << ")\n";
            //continue;
         }
         if (solver.info() == Success) {
            t_end = chrono::high_resolution_clock::now();
            cout << "[" << size << "] " << "Time to solve: "
                  << chrono::duration_cast<chrono::milliseconds>(
                        t_end - t_start).count()
                  << " ms. Estimated Error: " << solver.error() << ' '
                  << solver.iterations() << '/' << solver.maxIterations()
                  << "\n";
         }
         cout << endl;
      } else {
         ConjugateGradient<SparseMatrix<double>, Upper, IncompleteLUT<double> > solver;
         solver.setMaxIterations(100000);
         solver.setTolerance(10.0E-6);
         //solver.preconditioner().setDroptol(.001);
         solver.compute(mat);
         //solver.analyzePattern(mat);
         solver.factorize(mat);
         if (solver.info() != Success) {
            cout << "[" << size << "] Solver compute fail. ("
                  << solver.info() << ")\n";
            //continue;
         }
         cout<<"["<<size<<"] Start\n";
         X = solver.solve(B);
         //X = solver.solveWithGuess(B, Xe);
         if (solver.info() != Success) {
            cout << "[" << size << "] Solver solve fail. (" << solver.info()
                  << ")\n";
            //continue;
         }
         if (solver.info() == Success) {
            t_end = chrono::high_resolution_clock::now();
            cout << "[" << size << "] " << "Time to solve: "
                  << chrono::duration_cast<chrono::milliseconds>(
                        t_end - t_start).count()
                  << " ms. Estimated Error: " << solver.error() << ' '
                  << solver.iterations() << '/' << solver.maxIterations()
                  << "\n";
         }
         cout << endl;
      }

   }
}

int main() {
   SparseMatrix<double> mat(0, 0);
   loop(mat);
   return 0;
}


Performance wise, make sure you compiled with optimization ON. Then convergence with ILUT highly depends on the parameters you use. For symmetric problem, a simple Jacobi preconditioner (default one) might also be faster. For benchmarking purpose, you can try our bench routine: http://eigen.tuxfamily.org/dox-devel/gr ... tml#title3

ggael wrote:Performance wise, make sure you compiled with optimization ON. Then convergence with ILUT highly depends on the parameters you use. For symmetric problem, a simple Jacobi preconditioner (default one) might also be faster. For benchmarking purpose, you can try our bench routine: http://eigen.tuxfamily.org/dox-devel/gr ... tml#title3

Hello ggael and thanks for answering.
I've some preset matrices, some are simmetric, some aren't, but all of them are square and sparse.

I now have the same issue with a second program with a bit tidier code and I got the slowdowns and convergence issues again, I can't really understand where the issue is. It only solves the smaller non-symmetric matrix (square sparse 1891x1891), the others return solver.info() = 2. I don't have any issues with the symmetric ones. How should I set the solver and preconditioner to make it more likely to converge?

This is the new code:

Code: Select all

#include <iostream>
#include <fstream>
#include <vector>
#include <string>
#include <sstream>
#include <chrono>
#include <dirent.h>
#include <list>
#include "Eigen/SparseCore"
#include "Eigen/Sparse"
#include "Eigen/Dense"
#include "Eigen/IterativeLinearSolvers"

using namespace std;
using namespace Eigen;

typedef Triplet<double> T;

list<string> L;
char* path = "Matrici";
string non_ = "non-simmetrica-";
string sim_ = "simmetrica-";

static SparseMatrix<double> mat(0, 0);
static VectorXd X0, B, X;

bool is_dat(const string &str, const string &suffix) {
   return str.size() >= suffix.size()
         && str.compare(str.size() - suffix.size(), suffix.size(), suffix)
               == 0;
}
bool starts_with(const string &str, const string &prefix) {
   return str.size() >= prefix.size()
         && (str.compare(0, prefix.length(), prefix) == 0);
}
void listDir(const char* folder) {
   DIR *dir;
   struct dirent *ent;
   if ((dir = opendir(folder)) != NULL) {
      L.clear();
      while ((ent = readdir(dir)) != NULL) {
         if (ent->d_type == DT_REG) {
            L.push_back(ent->d_name);
         }
      }
      closedir(dir);
   } else {
      cout << "Error opening " << string(folder) << " dir." << endl;
      return;
   }
}

long populate_matrix(const string &file) {
   vector<T> tripletList;
   int size;
   bool first = true, simmetric = false;
   string val1, val2, val3;
   int ival1, ival2;
   double dval3;
   string line;
   ifstream iFile(file);
   while (getline(iFile, line)) {
      istringstream iss(line);
      getline(iss, val1, ' ');
      getline(iss, val2, ' ');
      getline(iss, val3, '\n');
      if (first) {
         first = false;
         size = stoi(val1);
         mat.resize(size, size);
         X0.resize(size);
         X.resize(size);
         B.resize(size);
         tripletList.reserve(stoi(val2));
         if (val3.compare("half"))
            simmetric = true;
         else if (val3.compare("full")) {
            ;
         } else {
            cout << "Error with file: " << file << endl;
            return -1;
         }
      } else {
         ival1 = stoi(val1) - 1; //row
         ival2 = stoi(val2) - 1; //column
         dval3 = stod(val3); //value
         if (simmetric && ival1 != ival2) {
            tripletList.push_back(T(ival2, ival1, dval3));
         }
         tripletList.push_back(T(ival1, ival2, dval3));
      }
   }
   auto start = chrono::high_resolution_clock::now();
   mat.setFromTriplets(tripletList.begin(), tripletList.end());
   auto end = chrono::high_resolution_clock::now();
   long elapsed =
         chrono::duration_cast<chrono::milliseconds>(end - start).count();
   cout << "[" << mat.rows() << "] " << "Time to populate " << elapsed
         << " ms.\n";
   iFile.close();

   return elapsed;
}

void generateX0() {
   //cout << "Generating X0..." << endl;
   for (int k = 0; k < X0.size(); k++) {
      X0(k) = k;
   }
}
void generateB() {
   auto t_start = chrono::high_resolution_clock::now();
   B = mat * X0;
   auto t_end = chrono::high_resolution_clock::now();
   cout << "[" << mat.rows() << "] " << "Time to find B by multiplication "
         << chrono::duration_cast<chrono::milliseconds>(t_end - t_start).count()
         << " ms.\n";
}

double relError() {
   return (X0 - X).squaredNorm() / X0.squaredNorm();
}

void simmetric_test(const string &file) {
}
void non_simmetric_test(const string &file) {
   int size = mat.rows();
   BiCGSTAB<SparseMatrix<double>, IncompleteLUT<double> > solver;
   //solver.setMaxIterations(100000);
   solver.setTolerance(10.0E-6);
   solver.preconditioner().setDroptol(.001);
   solver.compute(mat);
   if (solver.info() != Success) {
      cout << "[" << size << "] Solver compute fail. (" << solver.info()
            << ")\n";
      //continue;
   }
   cout << "[" << size << "] Solver Start\n";
   auto t_start = chrono::high_resolution_clock::now();
   X = solver.solve(B);
   //X = solver.solveWithGuess(B, Xe);
   if (solver.info() != Success) {
      cout << "[" << size << "] Solver solve fail. (" << solver.info()
            << ")\n";
      //continue;
   }
   if (solver.info() == Success) {
      auto t_end = chrono::high_resolution_clock::now();
      cout << "[" << size << "] " << "Time to solve: "
            << chrono::duration_cast<chrono::milliseconds>(t_end - t_start).count()
            << " ms. Estimated Error (built-in): " << solver.error()<<" Relative Error (norm2):"<<relError() << ' '
            << solver.iterations() << '/' << solver.maxIterations() << "\n";
   }
   cout << endl;
}

void elaboraMatrice(const string &file) {
   cout << "Testing " << path << "/" << file << " ..." << endl;
   cout << "================================================" << endl;
   string fPath = string(path).append("/").append(file);
   populate_matrix(fPath);
   generateX0();
   generateB();
   if (starts_with(file, non_)) {
      non_simmetric_test(fPath);
   } else if (starts_with(file, sim_)) {
      simmetric_test(fPath);
   }
   cout << endl;
}

int main() {
   listDir(path);
   list<string>::iterator i;
   for (i = L.begin(); i != L.end(); ++i) {
      if (is_dat(*i, string(".dat")))
         elaboraMatrice(*i);

   }
   return 0;
}

Hi, All, i have the same convergence issue too. It does not converage even for small matrix such as 200 * 200. The following is my code:

Eigen::BiCGSTAB<Eigen::SparseMatrix<double> , Eigen::IncompleteLUT<double>> BCGSTCOND;

BCGSTCOND.preconditioner().setFillfactor (10000);
BCGSTCOND.preconditioner().setDroptol(.001);
BCGSTCOND.setMaxIterations(iteration_max);
BCGSTCOND.setTolerance(converageTol);
BCGSTCOND.compute(*eigenSparMatrix_);
if(BCGSTCOND.info() != Eigen::Success)
return false;
Eigen::VectorXd x = BCGSTCOND.solve(b);
if(BCGSTCOND.info() != Eigen::Success)
return false;

Regarding BiCGSTAB, you might try the current 3.2 branch (https://bitbucket.org/eigen/eigen/get/3.2.tar.gz, future 3.2.2). I fixed a regression a few hours ago. Nonetheless, note that the convergence of the BiCGSTAB algorithm is not guaranteed and it depends a lot on the kind of problem and preconditionner.