# General matrix interative solver

Julian julvar at tamu.edu
Thu Oct 12 18:25:21 CDT 2006

```I am now using MatCreateSeqSBAIJ() - I am assuming the second parameter bs=1 in my case because I am using one processor. Is that correct ?
The method of creating a matrix I used was taken from Pavel Solin's book 'Partial Differential Equations and the finite Element Method'. Is there any particular reason why is should follow this method:
mat = (double*)malloc(sizeof(Mat))
MatCreateSeqSBAIJ(PETSC_COMM_SELF, 1, L, L, PETSC_DEFAULT, nnz, (Mat*)mat);
In the examples that I have seen in the petsc website I don't see this statement where memory is allocated.

Anyway, I changed my code to the following:
PetscInitialize(0,0,"petscOptions.txt",0);

int i;
int *nnz;
nnz=new int[L];
for(i=0;i<L;i++){
nnz[i]=row[i].getNum();
}

Assert( mat = (double*)malloc(sizeof(Mat)) );
MatCreateSeqSBAIJ(PETSC_COMM_SELF, 1, L, L, PETSC_DEFAULT, nnz, (Mat*)mat);
delete [] nnz;
And I'm still getting stuck at the same place... There is a difference though.. Now the assembly process is even slower than before (when there was no preallocation)
I noticed that it is getting stuck around the 200th MAT_FLUSH_ASSEMBLY call.. I don't know if that is significant.Is there some limit on how many times MAT_FLUSH_ASSEMBLY can be called ?

This is the petsc output I get:
[0] PetscInitializePETSc successfully started: number of processors = 1
[0] PetscGetHostNameRejecting domainname, likely is NIS ALPHA2.☺0☺
[0] PetscInitializeRunning on machine: ALPHA2
[0] PetscCommDuplicateDuplicating a communicator 1 1 max tags = 100000000

I also tried this :
mat = new Mat;
MatCreateSeqSBAIJ(PETSC_COMM_SELF, 1, L, L, PETSC_DEFAULT, nnz, mat);

And got the same performance.

> -----Original Message-----
> From: owner-petsc-users at mcs.anl.gov
> [mailto:owner-petsc-users at mcs.anl.gov] On Behalf Of Satish Balay
> Sent: Thursday, October 12, 2006 5:29 PM
> To: petsc-users at mcs.anl.gov
> Subject: RE: General matrix interative solver
>
>         MatCreateSeqAIJ(PETSC_COMM_SELF, L, L, PETSC_DEFAULT,
> nnz,(Mat*)mat);
>         MatSetType(*(Mat*)mat,MATSBAIJ);
>
> The above is a fuzzy way of creating a SBAIJ matrix. Use: either:
>
>     MatCreateSeqSBAIJ()
>
> or
>    MatCreate()
>    MatSetType(SBAIJ)
>    MatSeqSBAIJSetPreallocation()
>
> What you've done is:
>      Create SeqAIJ
>      set Preallocation for SeqAIJ
>      change MatType [SeqAIJ preallocation is of no use for SeqSBAIJ]
>
> Satish
>
> On Thu, 12 Oct 2006, Julian wrote:
>
> > Ok, so it looks like it is not preallocating.... This is
> what is get
> > when I do a '-info':
> >
> > [0] PetscInitializePETSc successfully started: number of
> processors =
> > 1 [0] PetscGetHostNameRejecting domainname, likely is NIS
> > ALPHA2.?????????????????
> >
> ??????????????????????????????????????????????????????????????
> ??????????????
> > ????
> > ????????????????????????????????????????????????M
> > [0] PetscInitializeRunning on machine: ALPHA2 [0]
> > PetscCommDuplicateDuplicating a communicator 1 1 max tags =
> 100000000
> > [0] MatSetUpPreallocationWarning not preallocating matrix storage
> >
> >
> > This is what I have been doing:
> >
> > 	PetscInitialize(0,0,"petscOptions.txt",0);
> >
> > 	int i;
> > 	int *nnz;
> > 	nnz=new int[L];
> > 	for(i=0;i<L;i++){
> > 		nnz[i]=row[i].getNum();
> > 	}
> >
> > 	Assert( mat = (double*)malloc(sizeof(Mat)) );
> > 	MatCreateSeqAIJ(PETSC_COMM_SELF, L, L, PETSC_DEFAULT, nnz,
> > (Mat*)mat);
> >
> > 	MatSetType(*(Mat*)mat,MATSBAIJ);
> > 	delete [] nnz;
> >
> > Is this not the right method to do preallocate?
> >
> > Julian.
> >
> > > -----Original Message-----
> > > From: owner-petsc-users at mcs.anl.gov
> > > [mailto:owner-petsc-users at mcs.anl.gov] On Behalf Of Satish Balay
> > > Sent: Thursday, October 12, 2006 2:33 PM
> > > To: petsc-users at mcs.anl.gov
> > > Subject: RE: General matrix interative solver
> > >
> > > On Thu, 12 Oct 2006, Julian wrote:
> > >
> > > > I think I got the preallocation right.. But it still
> gets stuck at
> > > > pretty much the same spot.
> > >
> > > Can you veify with '-info' - and look for number of
> mallocs in the
> > > verbose output.
> > >
> > > Satish
> > > > When it slowed to almost pretty much a full stop, I did a
> > > break all in
> > > > debug mode and it looks like its stuck in a memset
> function which
> > > > according to the call stack is called down the line from
> > > the following petsc function:
> > > > PetscErrorCode PETSC_DLLEXPORT PetscMallocAlign(size_t mem,int
> > > > line,const char func[],const char file[],const char dir[],void**
> > > > result)
> > > >
> > > > I am using the symmetric matrix (MATSBAIJ) ..and I'm
> not sure how
> > > > I should specify the nnz when I create the matrix Is it
> the number
> > > > of non-zeroes in each row in the upper triangle (cos
> that is all
> > > > I'm storing, right?) or is it the non zeroes in each
> row for the
> > > > whole matrix?
> > > >
> > > > Julian.
> > > >
> > > > > -----Original Message-----
> > > > > From: owner-petsc-users at mcs.anl.gov
> > > > > [mailto:owner-petsc-users at mcs.anl.gov] On Behalf Of Julian
> > > > > Sent: Thursday, October 12, 2006 1:47 PM
> > > > > To: petsc-users at mcs.anl.gov
> > > > > Subject: RE: General matrix interative solver
> > > > >
> > > > > Thanks, I got it working after I switched to icc
> > > > >
> > > > > All this time I was running test cases where I just
> > > > > 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
> > > > > 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
> > > > > > > -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
> > > > > > > > > > > > > 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.
> > > > > > > > > > > > > >
> > > > > > > > > > > > > >
> > > > > > > > > > > > >
> > > > > > > > > > > > >
> > > > > > > > > > > >
> > > > > > > > > > > >
> > > > > > > > > > >
> > > > > > > > > >
> > > > > > > > > >
> > > > > > > > >
> > > > > > > > >
> > > > > > > >
> > > > > > > >
> > > > > > >
> > > > > > >
> > > > > >
> > > > > >
> > > > >
> > > >
> > > >
> > >
> > >
> >
> >
>

```