General matrix interative solver

[Hong Zhang]

Julian,

When comparing the inhouse solver and petsc solver, you need
make sure

1. Timings are collected for solution process. The matrix and vector
   assembly should be excluded.
2. They should use same iterative algorithm. By default, petsc uses
   gmres with restart=30 and ilu(0) preconditioner. Petsc supports
   symmetric matrices, e.g., runtime option '-ksp_type cg -pc_type icc'
   might give better performance
3. They should start from same intial guess. By default, petsc
   initial guess is zero.

You can use '-ksp_view' to see what algorithm and options are
used in petsc.

Hong

On Wed, 11 Oct 2006, Julian wrote:

> Hello,
>
> I implemented the iterative sparse matrix solver in PetSc into my FEM code
> recently. I compared the results from a problem with 1317 unknowns. I used a
> direct solver to obtain the reference solution. I have another in-house
> sparse iterative solver that I have been using so far. It was written by
> someone else but I have access to the source for that solver.
> I find the 'error norm' in the solution by taking the square root of the sum
> of the squares of the absolute differences between the solution from the
> direct solver and the iterative solver. I am ignoring the numerical zeros in
> the solutions when doing this.
> I find that in order to get same order of the error norm (1e-13) as the
> in-house iterative solver, the petsc solver takes a much longer time and
> larger number of iterations. While the inhouse solver took less than one
> second, the petsc solver took 13 seconds. The inhouse solver took 476
> iterations whereas the petsc solver took 4738 iterations.
> I'm guessing this has to do with different setting of the solver in petsc
> such as the preconditioner etc.
> Can you tell me what the different settings are? And how to tweak them so
> that I can atleast get as good as a performance as the inhouse code ?
> Given below is how I have implemented the petsc solver:
>
> /////initialization
>     PetscInitializeNoArguments();
>         Assert( mat = (double*)malloc(sizeof(Mat)) );
>     MatCreateSeqAIJ(PETSC_COMM_SELF, L, L,
>       PETSC_DEFAULT, PETSC_NULL, (Mat*)mat);
>
> ////// this is the function I use to populate the matrix
> MatSetValue(*(Mat*)mat, ii, jj, value, ADD_VALUES);
>
> ////// this is how I actaully solve the matrix
>   MatAssemblyBegin(*(Mat*)mat, MAT_FINAL_ASSEMBLY);
>   MatAssemblyEnd(*(Mat*)mat, MAT_FINAL_ASSEMBLY);
>
>   double iter_error = 1e-10;
>   int max_iter_num = 10000;
>   int num_of_iter;
>
>         Vec rhs, x;
>     VecCreateSeqWithArray(PETSC_COMM_SELF, L, b, &rhs);
>     VecDuplicate(rhs, &x);
>
>     KSP ksp;
>     KSPCreate(PETSC_COMM_SELF, &ksp);
>     KSPSetTolerances(ksp, iter_error, PETSC_DEFAULT,
>       PETSC_DEFAULT, max_iter_num);
>     KSPSetFromOptions(ksp);
>     KSPSetOperators(ksp, *(Mat*)mat, *(Mat*)mat,
>       SAME_PRECONDITIONER);
>     KSPSolve(ksp,rhs,x);
>
>     PetscReal r_norm;
>     KSPGetResidualNorm(ksp, &r_norm);
>     KSPGetIterationNumber(ksp, &num_of_iter);
>
>         cout << "max_iter_num\t" << max_iter_num << endl;
>         cout << "iter_error\t" << iter_error << endl;
>
>         cout << "Matrix solver step " << num_of_iter << ", residual " <<
> r_norm << ".\n";
>
>         PetscScalar *p;
>     VecGetArray(x, &p);
>     for(int i=0; i<L; i++) {
>         b[i] = p[i];
>     }
>     VecRestoreArray(x, &p);
>
>     KSPDestroy(ksp);
>     VecDestroy(rhs);
>     VecDestroy(x);
>
> cout <<"Iterations for convergence="<< num_of_iter << " - Residual Norm = "
> << r_norm << endl;
>
>
>
> If this is not the typical method to be used to solve this kind of problem,
> please let me know what functions I should use.
> I should mention that the inhouse code is for symmetric matrices and from
> what I understand, the petsc solver works for general unsymmetric matrices.
> But I think for iterative solvers, it should still give around the same
> performance.
> I tested the solvers against some other problems as well, and I got the same
> performance.. In some cases, no matter how many iterations it goes through,
> the petsc solver would not go below a certain error norm whereas the inhouse
> solver would get almost exactly the same answer as the direct solver
> solution. I'm thinking the petsc solver should be able to solve this problem
> just as easily. It would be great if anyone could help me figure out the
> appropriate settings I should use in the petsc solver.
>
> Thanks,
> Julian.
>
>