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 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. > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > >
