BiCGSTAB does not solve Ax=b correctly

Summary

Hi everyone, I find the BiCGSTAB function yields wrong solution to the following Ax=b equation: [{1, 0, 0, 0}, {-1/6, 2/3, -1/6, -1/3}, {0, -1/3, 1, -2/3}, {0, -1/3, -1/3, 2/3}}] * x = [0, 1, 1, 1].

It yields [0, 16.3969, 23.0687, 19] but the correct answer should be [0, 15, 18, 18]. Note that the matrix is positive definite with determinant 2/27.

Environment

• Operating System : Windows
• Architecture : x64 ...
• Eigen Version : 3.4.0
• Compiler Version : Gcc8.1.0
• Compile Flags : -std=c++11

Minimal Example

#include <iostream>
#include <vector>
#include <string>
#include <sstream>
#include <cmath>
#include <map>
#include <cassert>
#include <algorithm>
#include <functional>
#include <numeric>
#include <unordered_map>
#include <set>
#include <Eigen/Dense>
#include <Eigen/IterativeLinearSolvers>
using namespace std;

vector<double> BiCGSTAB(const vector<vector<double>>& A, Eigen::VectorXd b) {
    Eigen::MatrixXd mat(A.size(), A.size());
    for (size_t i = 0; i < A.size(); ++i)
        for (size_t j = 0; j < A.size(); ++j)
            mat(i, j) = A.at(i).at(j);
    Eigen::BiCGSTAB<Eigen::MatrixXd> solver;
    solver.compute(mat);
    Eigen::VectorXd x = solver.solve(b);
    cout << mat << endl;
    cout << b << endl;
    cout << "Converged after " << solver.iterations() << " iterations with error " << solver.error() << " " << solver.info() << endl;
    cout << x << endl;
    vector<double> ret(A.size());
    for (size_t i = 0; i < A.size(); ++i)
        ret.at(i) = x[i];
    return ret;
}
int main()
{
    vector<vector<double>> A(4, vector<double>(4, 0.0));
    Eigen::VectorXd bb(4);
    A.at(0).at(0) = 1;
    A.at(1).at(0) = -1.0 / 6;
    A.at(1).at(1) = 2.0 / 3;
    A.at(1).at(2) = -1.0 / 6;
    A.at(1).at(3) = -1.0 / 3;
    A.at(2).at(1) = -1.0 / 3;
    A.at(2).at(2) = 1;
    A.at(2).at(3) = -2.0 / 3;
    A.at(3).at(1) = -1.0 / 3;
    A.at(3).at(2) = -1.0 / 3;
    A.at(3).at(3) = 2.0 / 3;
    bb[0] = 0;
    bb[1] = 1;
    bb[2] = 1;
    bb[3] = 1;
    BiCGSTAB(A, bb);
    return 0;
}

Steps to reproduce

1. first step
2. second step
3. ...

What is the current bug behavior?

Does not solve correctly.

What is the expected correct behavior?

Yield the correct answer.

Relevant logs

Warning Messages

Benchmark scripts and results

Anything else that might help

• Have a plan to fix this issue.