General matrix interative solver
Julian
julvar at tamu.edu
Thu Oct 12 13:46:32 CDT 2006
Thanks, I got it working after I switched to icc
All this time I was running test cases where I just read in an already
assembled global stiffness matrix and load vector and solved the problem
using petsc.
After I got that working, I tried to fully integrate the solver with the
FEA program.
One problem I had was with populating the matrix. The interfaces to the
other existing solvers in our code were written so that we can assign a
value directly to a memory space assigned for an element in the matrix.
Something like mat(I,J)=value; Where mat(I,J) gives the reference to the
location. But petsc does not work that way (I had already asked this on the
mailing list a couple of months back) and I understand maybe its wise to do
that, considering it handles parallel implementation as well. So, I changed
our code to work with petsc.
Now, when assembling the global matrix, sometimes I need to use ADD_VALUES
and other times INSERT_VALUES (depending on constraints)
At first, I used MatAssemblyBegin(*(Mat*)mat, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(*(Mat*)mat, MAT_FINAL_ASSEMBLY);
But then it took too long and I realized I should be using
MAT_FLUSH_ASSEMBLY. But even then it takes considerably more time in
assembling than my other solvers. So, I used a flag to call the assembly
only when I am switching between ADD_VALUES and INSERT_VALUES. That reduced
the calls to MAT_FLUSH_ASSEMBLY but it still was a slower than the other
solvers.. And this was for a small problem (1317 unknowns)
I tried it out with a larger problem ~15000 unknowns and the program just
got stuck during the assembly process ! It wouldn't crash but looks like its
doing something inside of petsc... probably memory allocations. So, I think
I need to start preallocating memory. I'm working on that now... So
hopefully, after that the assembly process won't slow down bcause of memory
allocation. But my question for now... is using the flag method to call
assembly only when needed the most efficient way to do this? Or do you have
any better ideas? One other possiblity is to have an array of flags, so as
to avoid the call to INSERT_VALUES altogether... And thereby eliminating the
need to call MAT_FLUSH_ASSEMBLY.
I'll get back to you after I get the preallocation right. Maybe that will
fix everything.
Thanks,
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 2:58 PM
> To: petsc-users at mcs.anl.gov
> Subject: RE: General matrix interative solver
>
>
> ILU is for general matrices; for symmetric you need to use
> PCSetType(pc,PCICC); or -pc_type icc
>
> Barry
>
>
> On Wed, 11 Oct 2006, Julian wrote:
>
> > Hong,
> > I made the following change to specify it to be in sbaij
> format... And
> > then proceeded to store only the upper triangle.
> >
> > Assert( mat = (double*)malloc(sizeof(Mat)) );
> > MatCreateSeqAIJ(PETSC_COMM_SELF, L, L, PETSC_DEFAULT,
> PETSC_NULL,
> > (Mat*)mat);
> > MatSetType(*(Mat*)mat,MATSBAIJ);
> >
> > But I get the following error when I try to solve:
> >
> >
> ----------------------------------------------------------------------
> > -- Petsc Release Version 2.3.1, Patch 16, Fri Jul 21
> 15:09:15 CDT 2006
> > HG
> > revision:
> > 5bc424fae3770001f0d316695a9956fde3bc58b0
> > See docs/changes/index.html for recent updates.
> > See docs/faq.html for hints about trouble shooting.
> > See docs/index.html for manual pages.
> >
> ----------------------------------------------------------------------
> > -- Unknown Name on a cygwin-cx named ALPHA2 by j0v1008 Wed Oct 11
> > 14:51:00 2006 Libraries linked from
> > /cygdrive/e/cygwin/home/petsc-2.3.1-p16/lib/cygwin-cxx-rea
> > l-debug
> > Configure run at Mon Aug 14 15:35:57 2006 Configure options
> > --with-cc="win32fe cl --nodetect" --with-cxx="win32fe cl
> --nod etect"
> > --COPTFLAGS="-MDd -Z7" --with-fc=0 --with-clanguage=cxx
> > --download-c-blas
> > -lapack=1 --with-mpi=0 --useThreads=0 --with-shared=0
> >
> ----------------------------------------------------------------------
> > -- [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!
> > [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