GMRES solver

[cyril]

Hello. Thanks for one more great library. I'm trying to solve a system of linear equations originating from convection equation (cell-centered finite volume method). I get the right solution for a simple test case:

Fig. 1. Test case 1. Parameters of matrix K (K C = L): dimensions 3148x3148, nonsymmetric, cond.number=152, rank=3148, det=0 (MATLAB command det(A) gives such a result)

But I get the wrong solution for a little more complex test case:

Fig. 2. Test case 2. Parameters of matrix K (K C = L): dimensions 10901x10901, nonsymmetric, cond.number=5.92E3, rank=10901, det=0

Fig. 3. Test case 2. MATLAB solution

I call ViennaCL as

#define VIENNACL_WITH_OPENMP

// ...

void mexFunction(int nlhs, mxArray *plhs[],
    int nrhs, const mxArray *prhs[])
{

    // ...

    int N_th = mxGetScalar(prhs[5]);
    omp_set_num_threads(N_th);

    // ...

    std::vector< std::map< unsigned int, ScalarType> > K(N_dof);// stiffness matrix
    std::vector<ScalarType> L(N_dof);// load vector
    std::vector<ScalarType> C(N_dof);// solution vector

    // assemble stiffness matrix and load vector

    viennacl::compressed_matrix<ScalarType> vcl_K;
    viennacl::vector<ScalarType> vcl_L(N_dof);
    viennacl::vector<ScalarType> vcl_C(N_dof);

    viennacl::copy(K, vcl_K);
    viennacl::copy(L.begin(), L.end(), vcl_L.begin());

    viennacl::linalg::ilu0_precond< viennacl::compressed_matrix<ScalarType> > vcl_ilu0(vcl_K, viennacl::linalg::ilu0_tag());

    viennacl::linalg::gmres_tag my_gmres_tag(1e-6, 100, 20);
    vcl_C = viennacl::linalg::solve(vcl_K, vcl_L, my_gmres_tag, vcl_ilu0);

    viennacl::copy(vcl_C.begin(), vcl_C.end(), C.begin());

    // ...
}

I've also tried without preconditioner and with different number of maximum iterations but without success. What should I try to understand why ViennaCL doesn't solve the linear system of test case 2?

Thanks,
Cyril

[Karl Rupp]

Hi Cyril,
in viennacl::linalg::gmres_tag my_gmres_tag(1e-6, 100, 20); you set a maximum iteration number of 100. This is likely to be insufficient for the larger problem. Try viennacl::linalg::gmres_tag my_gmres_tag(1e-8); for the default values on iterations and a tighter error bound.
Best regards,
Karli

[cyril]

Hello Karli,

I've tried your solver parameters and got the same wrong solution. I didn't mention in the first massage that I've already "played" with different parameters of the solver but without success. I've also tried to solve with transposed stiffness matrix and without preconditioner.

Why GMRES solves test case 1 though MATLAB tells me that my matrix has zero determinant?

Thanks,
Cyril

[Karl Rupp]

Hi Cyril,

ah, good catch, the determinant of zero tells us that the system matrix is singular. That is, there is a nontrivial kernel vector v for which Av = 0. For any solution x of Ax=b, the vectors x + a * v for any scalar a are also a solution of the system. Matlab might solve for the minimum norm solution, i.e. min_a ||x + av||, whereas GMRES in general does not compute the minimum-norm solution.

A possible cause of the singularity are incorrect boundary conditions (or no boundary conditions at all). For finite volume schemes I expect to see a nonzero determinant.

Best regards,
Karli

[cyril]

Thanks for the answer. I have to think over it.

[cyril]

GMRES of ViennaCL works great. The stupid error was connected with results plotting.