General matrix interative solver
Hong Zhang
hzhang at mcs.anl.gov
Thu Oct 12 09:02:34 CDT 2006
> ------------------------------------------------------------------------
> [0]PETSC ERROR: MatILUFactorSymbolic() line 4573 in
> src/mat/interface/e:\cygwin\
> home\PETSC-~1.1-P\src\mat\INTERF~1\matrix.c
> [0]PETSC ERROR: No support for this operation for this object type!
> [0]PETSC ERROR: Matrix type sbaij symbolic ILU!
Using sbaij format,
you must use preconditioner icc, not ilu.
Hong
> [0]PETSC ERROR: PCSetUp_ILU() line 532 in
> src/ksp/pc/impls/factor/ilu/e:\cygwin\
> home\PETSC-~1.1-P\src\ksp\pc\impls\factor\ilu\ilu.c
> [0]PETSC ERROR: PCSetUp() line 798 in
> src/ksp/pc/interface/e:\cygwin\home\PETSC-
> ~1.1-P\src\ksp\pc\INTERF~1\precon.c
> [0]PETSC ERROR: KSPSetUp() line 234 in
> src/ksp/ksp/interface/e:\cygwin\home\PETS
> C-~1.1-P\src\ksp\ksp\INTERF~1\itfunc.c
> [0]PETSC ERROR: KSPSolve() line 334 in
> src/ksp/ksp/interface/e:\cygwin\home\PETS
> C-~1.1-P\src\ksp\ksp\INTERF~1\itfunc.c
>
> > -----Original Message-----
> > From: owner-petsc-users at mcs.anl.gov
> > [mailto:owner-petsc-users at mcs.anl.gov] On Behalf Of Hong Zhang
> > Sent: Wednesday, October 11, 2006 2:14 PM
> > To: petsc-users at mcs.anl.gov
> > Subject: RE: General matrix interative solver
> >
> >
> > 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