September 2011 - petsc-users - mailman3.cels.anl.gov

ML - zero pivot error
by Stijn A. M. Vantieghem 17 Jul '12

17 Jul '12

Dear all, I am using PETSc (Fortran interface) to solve a Poisson equation with Neumann boundary conditions. Up till now, I managed to do this with Hypre's BoomerAMG. Now, I am investigating whether I can improve the performance of my code by using ML. However, at execution I receive a zero pivot error; I tried to remove the (constant) null space with KSPSetNullSpace, but this didn't solve my problem. Do you have an idea of what I'm doing wrong? Thanks. The relevant portions of my code are as follows: !**************************************************** call KSPCreate(PETSC_COMM_WORLD,ksp,ierr) call KSPSetOperators(ksp,M,M,DIFFERENT_NONZERO_PATTERN,ierr) call MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,PETSC_NULL,sp,ierr) call KSPSetNullSpace(ksp,sp,ierr) call KSPGetPC(ksp,pc,ierr) call PCSetType(pc,PCML,ierr) call KSPSetFromOptions(ksp,ierr) ... call KSPSetInitialGuessNonzero(ksp,PETSC_TRUE,ierr) call KSPSolve(ksp,petsc_rhs,petsc_pressure,ierr) !**************************************************** and the error message is: [0]PETSC ERROR: --------------------- Error Message ------------------------------------ [0]PETSC ERROR: Detected zero pivot in LU factorization see http://www.mcs.anl.gov/petsc/petsc-as/documentation/troubleshooting.html#Ze…! [0]PETSC ERROR: Zero pivot row 0 value 3.57045e-20 tolerance 1e-12! [0]PETSC ERROR: ------------------------------------------------------------------------ ... [0]PETSC ERROR: ------------------------------------------------------------------------ [0]PETSC ERROR: MatLUFactorNumeric_SeqAIJ() line 574 in src/mat/impls/aij/seq/aijfact.c [0]PETSC ERROR: MatLUFactorNumeric() line 2587 in src/mat/interface/matrix.c [0]PETSC ERROR: PCSetUp_LU() line 158 in src/ksp/pc/impls/factor/lu/lu.c [0]PETSC ERROR: PCSetUp() line 795 in src/ksp/pc/interface/precon.c [0]PETSC ERROR: KSPSetUp() line 237 in src/ksp/ksp/interface/itfunc.c [0]PETSC ERROR: PCSetUp_MG() line 602 in src/ksp/pc/impls/mg/mg.c [0]PETSC ERROR: PCSetUp_ML() line 668 in src/ksp/pc/impls/ml/ml.c [0]PETSC ERROR: PCSetUp() line 795 in src/ksp/pc/interface/precon.c [0]PETSC ERROR: KSPSetUp() line 237 in src/ksp/ksp/interface/itfunc.c [0]PETSC ERROR: KSPSolve() line 353 in src/ksp/ksp/interface/itfunc.c Regards Stijn Stijn A.M. Vantieghem Earth and Planetary Magnetism Institute for Geophysics ETH Zürich Sonneggstrasse 5 - CH 8092 Zürich tel: +41 44 632 39 90 e-mail: stijn.vantieghem(a)erdw.ethz.ch

5 7

Re: Multiple communicators
by irfan.khan＠gatech.edu 02 Apr '12

02 Apr '12

I use seperete communicators for fluid and solid ranks in a fluid-structure interaction code. I use PETSc tools (KSP) to solve for the solid phase (FEA) which is carried out under a different communicator (FEA_Comm). Unless I am doing something wrong, I have found that the command line options: -ksp_type; -mat_partitioning_type, don't work if there are 2 communicators. But for single communicator they work. Please see the attached files containing output of -log_summary for 2 different codes (scroll down to the end). One with single communicator and the other with two communicators. Both the output have been generated with the same PETSc compilation. Futher -log_summary results are not able to detect all the events that have been registered if they are not in the WORLD communicator. Thank you Irfan ----- Original Message ----- From: "Matthew Knepley" <knepley(a)gmail.com> To: "PETSc users list" <petsc-users(a)mcs.anl.gov> Sent: Tuesday, November 24, 2009 8:03:00 PM GMT -05:00 US/Canada Eastern Subject: Re: Multiple communicators On Tue, Nov 24, 2009 at 6:59 PM, < irfan.khan(a)gatech.edu > wrote: Hello Does the procedure for providing command line options and to obtain profiling data through -log_summary change if there are multiple MPI communicators in a petsc code? There is always WORLD, so it should be fine unless I misunderstand your question. Matt If so, can somebody point me to information, references on how to provide command line options if there are multiple communicators? Thank you Irfan -- PhD Candidate G.W. Woodruff School of Mechanical Engineering Georgia Institute of Technology Atlanta, GA (30332) -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener -- PhD Candidate G.W. Woodruff School of Mechanical Engineering Georgia Institute of Technology Atlanta, GA (30332)

5 14

Re: [petsc-users] direct access of seqcusp Vec
by Shiyuan 03 Oct '11

03 Oct '11

What about a direct access to seqaijcusp Mat? Is there a function call which return the pointed to the cusp::csr_matrix ? Also, i notice that, the flops count from --log_summary are different if we switch on/off -vec_type seqcusp and -mat_type seqaijcusp? Why is that? what is the option -cusp_synchronize for? Thanks. ------------------------------ Message: 2 Date: Wed, 28 Sep 2011 23:22:02 -0500 From: Barry Smith <bsmith(a)mcs.anl.gov> Subject: Re: [petsc-users] direct access of seqcusp Vec To: PETSc users list <petsc-users(a)mcs.anl.gov> Message-ID: <1CFA4EF7-E716-4942-BE0F-6AB2E4EF78E5(a)mcs.anl.gov> Content-Type: text/plain; charset=us-ascii VecCUSPGetArrayReadWrite() gives you access to the CUSP array that lives on the device memory. You can then do things to it using the CUSP header file definitions Barry On Sep 28, 2011, at 10:46 PM, Shiyuan wrote: > Hi, > How can we acess the data of a seqcusp Vec directly? Is it stored as consecutive 1d array in the device momery? > Thanks. > > Shiyuan ------------------------------

2 3

question about contexts
by Dominik Szczerba 01 Oct '11

01 Oct '11

I am currently declaring functions outside my C++ solver class to be used by via MatShellSetOperation like: #undef __FUNCT__ #define __FUNCT__ "StiffnessMatrixMultiplication" PetscErrorCode StiffnessMatrixMultiplication(Mat A, Vec x, Vec y) { // ... ierr = MatShellGetContext(A, (void**) &ctx);CHKERRQ(ierr); Mat foo; foo = ctx->foo; // ... } and I set/get contexts to e.g. pass 'foo' around. Is it possible in Petsc / is it a good idea / to instead declare such functions as static class member functions and then use the class members directly without any contexts shuffled around? Thanks for any hints. Dominik

3 14

MatMult a MatNest and a Vec
by Hui Zhang 01 Oct '11

01 Oct '11

Hi all, Do you know what type of the resulting vector y by MatMult of a MatNest A and a Vec x, i.e. y = Ax? Because after that I want to modify a SubVec of the result y. If it is a VectNest, I think VecNestGetSubVec should be used but not VecGetSubVector. Thank you for always helps! Hui

3 4

petsc-3.2-p3.tar.gz now available
by Satish Balay 30 Sep '11

30 Sep '11

Dear PETSc users, Patch-3 update to petsc-3.2 is now available for download. http://www.mcs.anl.gov/petsc/petsc-as/download/index.html Satish

1 0

Re: [petsc-users] direct access of seqcusp Vec
by Barry Smith 29 Sep '11

29 Sep '11

For VecNorm_MPICUSP() it called VecNorm_SeqCUSP() which logged the flops but then VecNorm_MPICUSP() logged it again. Hence they were double counted. I will push an update to petsc-dev in a few minutes that fixes this. MatMult_SeqAIJ() attempts to be more careful in the count by taking into account that each completely empty row contributes one less flop then expected which the CUSP cannot easily do that since we don't know how the CUSP actually implements the product for empty rows and are too lazy to check. I cannot explain why the VecMDot are slightly different. Since PETSc is all open source and we value improvements by our users you are free to look more closely to see if you can determine why the mdot values are slighly different. Thanks for reporting the bug, Barry On Sep 29, 2011, at 4:36 PM, Shiyuan wrote: > > > Also, i notice that, the flops count from --log_summary are different if we switch on/off -vec_type seqcusp and -mat_type seqaijcusp? > > Are you sure the exact same solver is being used and the exact same number of iterations. Check the -log_summary table for the number of times > each operation (vecdot, matmult etc etc) is done for the cusp and non-cusp case, if they are different then you will get different counts. > > I have two runs of ex19, exactly same options except -da_vec_type seqcusp and -da_mat_type are switch on/off. The count of each operation are the same but the flops count for some operations are different. > > > VecMDot 2024 1.0 4.1337e+00 1.0 2.50e+09 1.0 0.0e+00 0.0e+00 0.0e+00 17 29 0 0 0 19 29 0 0 0 605 > VecNorm 2096 1.0 9.4352e-01 1.0 1.68e+08 1.0 0.0e+00 0.0e+00 0.0e+00 4 2 0 0 0 4 2 0 0 0 178 > MatMult 2092 1.0 1.2145e+01 1.0 3.24e+09 1.0 0.0e+00 0.0e+00 0.0e+00 50 37 0 0 0 56 37 0 0 0 267 > > vs. > VecMDot 2024 1.0 1.0873e+00 1.0 2.54e+09 1.0 0.0e+00 0.0e+00 0.0e+00 17 28 0 0 0 28 28 0 0 0 2340 > VecNorm 2096 1.0 1.2240e-01 1.0 3.35e+08 1.0 0.0e+00 0.0e+00 0.0e+00 2 4 0 0 0 3 4 0 0 0 2740 > MatMult 2092 1.0 1.7572e+00 1.0 3.32e+09 1.0 0.0e+00 0.0e+00 0.0e+00 28 37 0 0 0 46 37 0 0 0 1890 > > #PETSc Option Table entries: > -cusp_synchronize > -da_grid_x 100 > -da_grid_y 100 > -da_mat_type mpiaijcusp > -da_vec_type mpicusp > -dmmg_nlevels 1 > -log_summary > -mat_no_inode > -pc_type none > -preload off > #End of PETSc Option Table entries > > The whole -log_summary is attached. > Thanks. > <Perf_ex19_cusp_0><Perf_ex19>

1 0

Question about ksp ex3.c
by Matija Kecman 29 Sep '11

29 Sep '11

Hi all, I have been playing around with ksp/examples/tutorials/ex3.c, a brief summary of this example follows: The Laplace equation in two-dimensions on a unit square [0,1]x[0,1] is solved using a structured mesh of linear quadrilateral finite elements, m elements in each direction. Dirichlet boundary conditions are specified: u_hat = y on the boundary, which means the exact solution is u_star = y. The global system matrix dimension is specified and PETSc decides how to split the matrix amongst processors. The DOFs are numbered 'naturally' from the origin, left to right, bottom to top. Boundary conditions are applied by zeroing the rows of the matrix corresponding to the boundary DOFs, inserting 1's on the diagonal entries of these rows and inserting the boundary conditions into the global RHS vector entries corresponding to these DOFs. I ran the following uniprocess tests (to run the direct solve tests I had to comment the line 'ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);CHKERRQ(ierr);' but otherwise I made no changes to the example code): ---Tests with 5x5 element mesh (36 DOFs)--- $ mpirun -np 1 ./ex3 -ksp_type gmres Norm of error 1.73622e-06 Iterations 3 $ mpirun -np 1 ./ex3 -ksp_type bcgs Norm of error 1.5157e-07 Iterations 2 $ mpirun -np 1 ./ex3 -ksp_type cg Norm of error 1.88557e-06 Iterations 3 $ mpirun -np 1 ./ex3 -ksp_type gmres -pc_type none Norm of error < 1.e-12 Iterations 7 $ mpirun -np 1 ./ex3 -ksp_type bcgs -pc_type none Norm of error 0.00109538 Iterations 6 $ mpirun -np 1 ./ex3 -ksp_type cg -pc_type none Norm of error < 1.e-12 Iterations 7 $ mpirun -np 1 ./ex3 -ksp_type preonly -pc_type lu Norm of error < 1.e-12 Iterations 1 ---Tests with 100x100 element mesh (10201 DOFs)--- $ mpirun -np 1 ./ex3 -ksp_type gmres -m 100 Norm of error 2.90652e-05 Iterations 32 $ mpirun -np 1 ./ex3 -ksp_type bcgs -m 100 Norm of error 3.92733e-05 Iterations 24 $ mpirun -np 1 ./ex3 -ksp_type cg -m 100 Norm of error 2.4332e-05 Iterations 32 $ mpirun -np 1 ./ex3 -ksp_type gmres -pc_type none -m 100 Norm of error 0.570146 Iterations 0 $ mpirun -np 1 ./ex3 -ksp_type bcgs -pc_type none -m 100 Norm of error 0.570146 Iterations 0 $ mpirun -np 1 ./ex3 -ksp_type cg -pc_type none -m 100 Norm of error 0.570146 Iterations 0 $ mpirun -np 1 ./ex3 -ksp_type preonly -pc_type lu -m 100 Norm of error < 1.e-12 Iterations 1 I had two questions: 1. I was wondering if someone could explain why the system cannot be solved using for a mesh of 100x100 elements without a using preconditioner? The way the boundary conditions are applied causes the system matrix to lose symmetry but GMRes and BiCGStab should work for non-symmetric systems. 2. I'm also unsure about why the CG method seems to be able to solve this non-symmetric system in some cases, are there some cases where this is possible? Thanks, Matija

4 8

Matrix Reordering
by behzad baghapour 29 Sep '11

29 Sep '11

Dear all, How can I use the result ( rperm and cperm ) of MatGetOrdering( Mat mat, const MatOrderingType type, IS *rperm, IS *cperm ) as an array of integers in order to permute my own Matrix? Thanks a lot, Behzad -- ================================== Behzad Baghapour Ph.D. Candidate, Mechecanical Engineering University of Tehran, Tehran, Iran https://sites.google.com/site/behzadbaghapour Fax: 0098-21-88020741 ==================================

2 1

Memory usage of KSP (or any PETSc object)
by Mohamad M. Nasr-Azadani 28 Sep '11

28 Sep '11

Hi, I was trying to find the exact (or more or less exact) memory consumption of any PETSc object. I was wondering how accurate is this method: *PetscMemoryGetMaximumUsage(&mem1);* *.* *(creation of any PETSc object).* *.* * * *PetscMemoryGetMaximumUsage(&mem2);* *Memory = mem2 - mem1; * * * *Thanks in advance, * *Mohamad* * *

4 7