General matrix interative solver

[Barry Smith]

  Julian,

   If your matrix is symmetric, positive definite then you should
ALWAYS use -ksp_type cg to use the conjugate gradient method
instead of the more expensive (and unneeded) GMRES.

  Next, in testing you can/should run with -ksp_view to have the
exact solver options printed out to the screen.

  Running with -help will have it print the current solver options.

  Run with -ksp_converged_reason to see if the solver converged or
if it failed and why.

  Next you need to determine what the inhouse solver is doing and
try to set PETSc to use the same solver.

  This problem is very, very small for an iterative solver; recommend
iterative solvers only when working with, say 20,000+ unknowns (depends
on matrix conditioning, number of nonzero etcs also). For small problems
direct solver is the way to go. You can use the PETSc direct solver
with -ksp_type preonly -pc_type lu (or for symmetric positive definite
use -pc_type cholesky).

  Let us know what progress you make, it is imperative you determine
what the in house solver is using.

   Barry

 

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.
> 
>