General matrix interative solver

[Julian]

You are right.. This is a very small problem... I eventually plan on using
it for very large problems.
But for now, I wanted to confirm that I am using the right settings and I
get an accurate enough solution by comparing it with the direct solver
solution. 

Currently, with the way I have integrated it with my fea code, I am unable
to send the command line arguments to the solver. Is there some easy way to
send the arguments to the solver during command time? 

Julian.


> -----Original Message-----
> From: owner-petsc-users at mcs.anl.gov 
> [mailto:owner-petsc-users at mcs.anl.gov] On Behalf Of Barry Smith
> Sent: Wednesday, October 11, 2006 10:22 AM
> To: petsc-users at mcs.anl.gov
> Subject: Re: General matrix interative solver
> 
> 
>   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.
> > 
> > 
> 
>