General matrix interative solver
Hong Zhang
hzhang at mcs.anl.gov
Wed Oct 11 14:13:40 CDT 2006
Julian,
If you do not perform matrix reordering, you can use
petsc sbaij format, which only stores the upper triangular
entries. aij format has more algorithmic options and is
more efficient and supports matrix reordering.
Hong
On Wed, 11 Oct 2006, Julian wrote:
> I'm using '-ksp_type cg' now to solve the symmetric matric...But then do I
> need to store all the non-zero elements of a symmetric matrix ? i.e. A[i,j]
> as well as A[j,i] ?
>
> Shouldn't I be able to store just the upper or lower triangle or the matrix
> ? Is that possible with petsc ?
>
> Julian.
>
>
> > -----Original Message-----
> > From: owner-petsc-users at mcs.anl.gov
> > [mailto:owner-petsc-users at mcs.anl.gov] On Behalf Of Julian
> > Sent: Wednesday, October 11, 2006 12:37 PM
> > To: petsc-users at mcs.anl.gov
> > Subject: RE: General matrix interative solver
> >
> > Barry,
> >
> > I tried sending the commands from a file for now... And once
> > I used '-ksp_type cg', I got pretty much the same performance
> > as that of the inhouse code !
> > Now, I'll try the performance with larger cases...
> >
> > Thanks for your help!
> >
> > 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 11:25 AM
> > > To: petsc-users at mcs.anl.gov
> > > Subject: RE: General matrix interative solver
> > >
> > >
> > > PetscOptionsSetValue() but recommend putting them on the command
> > > line or in a file (pass the name of the file into
> > PetscInitialize()).
> > > Having to recompile for every option is too painful.
> > >
> > > Barry
> > >
> > >
> > > On Wed, 11 Oct 2006, Julian wrote:
> > >
> > > > Hong,
> > > >
> > > > The times are for the solution process alone and the
> > > initial guess is
> > > > the same, i.e. zero.
> > > >
> > > > As for the algorithm, they are probably different.
> > > > Can you tell me how to pass 'runtime options' like -help
> > > and -ksp_view
> > > > from within the code... I mean, at compile time.
> > > >
> > > > thanks
> > > >
> > > > >
> > > > > 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.
> > > > > >
> > > > > >
> > > > >
> > > > >
> > > >
> > > >
> > >
> >
> >
>
>
More information about the petsc-users
mailing list