Hot answers tagged

42

Newton's method may not converge for many reasons, here are some of the most common. The Jacobian is wrong (or correct in sequential but not in parallel). The linear system is not solved or is not solved accurately enough. The Jacobian system has a singularity that the linear solver is not handling. There is a bug in the function evaluation routine. The ...


25

Initial advice Always run with -ksp_converged_reason -ksp_monitor_true_residual when trying to learn why a method is not converging. Make the problem size and number of processes as small as possible to demonstrate the failure. You often gain insight by determining what small problems exhibit the behavior that is causing your method to break down and the ...


19

There are huge differences in culture, coding style, and capabilities. Probably the fundamental difference is Trilinos tries to provide an environment for solving FEM problems and PETSc provides an environment for solving sparse linear algebra problems. Why is that significant? Trilinos will provide a large number of packages concerned with separate ...


15

My advice to students is to try a direct solver in these cases. The reason is that there are two classes of reasons why a solver may not converge: (i) the matrix is wrong, or (ii) there is a problem with the solver/preconditioner. Direct solvers almost always yield something that you can compare with the solution you expect, so if the answer of the direct ...


14

Matrix Market is a terrible format for reading in parallel, therefore it is better to preprocess to a better parallel format. Your matrix size is extremely small so performance is not an issue, but the easiest and most general thing is to use Python or Matlab/Octave to write the Matrix Market file in PETSc binary format, which can be read efficiently in ...


13

It is very expensive to compute the Schur complement of a matrix and very rarely needed in practice. PETSc highly recommends avoiding algorithms that need it. The Schur complement of a matrix (dense or sparse) is essentially always dense, so begin by: forming a dense matrix $K_{ba}$, also create another dense matrix $T$ of the same size. Then either ...


13

This can arise from architectural factors: If the same memory bandwidth is available for both one or more processes, then you will see almost no speedup since SpMV and related linear algebra operations are memory bandwidth limited. It might also be the case that communication overhead overwhelms you local computation. For example, in linear iterative ...


12

Warning Solving saddle point problems involves a lot more choices than definite problems, and there are a lot more things that can go wrong. Use monitors for all levels to debug convergence, to sure that null spaces are handled correctly when auxiliary operators are singular (usually just a constant null space), and to ensure that preconditioners are ...


11

For small matrices, the condition number can be reliably computed using the singular value decomposition. Do a KSPSolve() with the matrix and run with -pc_type svd -pc_svd_monitor. For larger matrices, you can estimate the condition number using a Krylov method. For example, the Arnoldi iteration performed by GMRES incrementally computes a Hessenberg ...


11

PETSc multigrid (as a preconditioner) is quite mature and may be used with any of the KSP (iterative Krylov method) solvers in PETSc by typing: -pc_type mg However, this requires that you have some way of generating your coarse levels, such as having structured grids defined by PETSc DA objects, which will be coarsened automatically. Or, if you want to ...


11

This is a widely-held concern in the scientific programming community, and I would consider the performance uncertainty to be one of the major "myths" in computational science. As @fcruz discusses, petsc4py is a wrapper to the PETSc libraries, not a reimplementation of PETSc in Python. Therefore, you can expect any performance penalties to come from ...


10

PETSc uses BLAS for a few vector primitives, but these are generally limited by memory bandwidth and there isn't much variance in "optimization", so it tends not to make much performance difference. It also uses Lapack for some analysis such as Lanczos or Arnoldi estimates of eigenvalues and singular values, but these are generally not performance-sensitive....


10

Without taking sides the discussion about whether to use direct or iterative solvers, I just want to add two points: There exist Krylov methods for systems with multiple right-hand sides (called block Krylov methods). As an added bonus, these often have faster convergence than standard Krylov methods since the Krylov space is built from a larger collection ...


10

I've been developing a lightweight finite element library in Python 2.7 harnessing the power of NumPy arrays and SciPy sparse matrices. The general idea is that given a mesh and a finite element, you have more-or-less one-to-one correspondence between the bilinear form and a (sparse) matrix. The user can then use the resulting matrix as he or she sees fit. ...


9

There is a lot of what I'm guessing is PETSc-specific language in your question (with which I am unfamiliar), so there might be a wrinkle here I don't quite understand, but maybe this will still be useful to get you started. Basically, you have to define the interface for the procedure, and then you can pass a pointer to a function that follows this ...


9

There are a few ways to do this: Check $PETSC_ARCH/conf/reconfigure-$PETSC_ARCH.py (where $PETSC_ARCH is expanded, i.e. for me that is 'arch-c') Depending on which information you want, you could use make getlinklibs or make getincludedirs Also, you could check $PETSC_ARCH/include/petscconf.h for all the standard #ifdef's that PETSc was configured with


9

The general advice given on the mailing list is usually about efficient assembly. The short answer is, no; this is a fundamental limitation. The longer answer is that, for one, there is no "end-all-be-all" sparse matrix implementation. In general, it will be faster to: loop over your mesh / domain, counting how many entries you'll need in your assembled ...


9

There is typically a trade-off between the amount of work you put into constructing a good preconditioner for an iterative solver and the work you save by using a good preconditioner when actually solving the linear systems. In your case, the case is pretty clear: put as much work as you can into constructing a good preconditioner because you have to solve ...


9

This is usually caused by trying to use a threaded MKL combined with MPI, resulting in over-subscription. Either explicitly configure PETSc to use non-threaded MKL or add MKL_NUM_THREADS=1 to your environment.


8

Whenever trying to parallelise a program you have to balance out a number of costs, but primarily there is The cost of running each computation The cost of any communications between those computations The cost of managing those computations If your computations are embarrassingly parallel then the communications cost will be very low (input and output ...


8

I'm a happy user of GoogleTest with a C++ MPI code in a CMake/CTest build environment: CMake automatically installs/links googletest from svn! adding tests is a one-liner! writing the tests is easy! (and google mock is very powerful!) CTest can pass command-line parameters to your tests, and exports data to CDash! This is how it works. A batch of unit-...


8

As one of the library's authors, I would of course love for deal.II to come out on top with this comparison. But I suspect it may not, and the answer lies in a factor you omit from your comparison: how long it actually takes to implement your code. Few people in academia with the skills to implement a FEM code from scratch spend more time solving PDEs than ...


7

Both PETSc and Trilinos have good algebraic multigrid methods. deal.II implements geometric multigrid methods for finite element discretizations, see for example the step-16 tutorial program.


7

If you have dense matrices with structure (e.g. fast transforms, Schur complements, etc), PETSc could be useful. In these cases, you won't be assembling the full matrix. For assembled dense systems, PETSc currently uses PLAPACK, but the matrix distribution in PETSc native format is not the best to minimize communication (for most operations). Jack Poulson, ...


7

Let me add to aterrel's good comment that Trilinos is really a big bag of (Sandia) stuff and Petsc is a more focused library. If you want to compare then you should compare PETSc's sparse solver support with Trilinos' ePetra/ML/etc sparse solver ecosystem, that do similar things. Also, PETSc supports structured grids and Sandia has historically explicitly ...


7

Petsc4py is just another way to acces PETSc but from python, or is the same to say that, petsc4py provides the bindings so that, from python, you can access PETSc data structures and routines that are meant to reduce the effort of developing parallel PDE solvers (that scale). PETSc provides several levels of abstractions to their solvers, and you can even ...


7

It is important to preallocate correctly. This is almost certainly the reason why your assembly was slow. If you are starting with a dense matrix representation, just scan through it once counting the number of nonzeros per row, then call MatSeqAIJSetPreallocation(). See this FAQ. The option MAT_IGNORE_ZERO_ENTRIES is really intended to be used when there is ...


7

Michael Pippig at the University of Chemnitz (Germany) has implemented an MPI-parallelized FFT that uses FFTW in the background. This might help you: http://www-user.tu-chemnitz.de/~mpip/software.php?lang=en It is using the algorithm proposed by Plimpton from Sandia National Labs as suggested by Eldila's comment.


7

You'll need to roll your own preconditioner. If you know the matrix, it should not be terribly difficult to implement something like an SSOR preconditioner, for example. If you know something else about the problem, for example that it comes from a PDE whose solution can be well approximated on a coarser mesh, then you can also consider constructing ...


Only top voted, non community-wiki answers of a minimum length are eligible