Hot answers tagged

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

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


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


9

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

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


6

Assuming that your structures are actually 3D (rather than only thin features, perhaps discretized with shell elements) and that the model is larger than a few hundred thousand dofs, direct solvers become impractical, especially if you only need to solve each problem once. Additionally, unless the structure is always "close" to a Dirichlet boundary, you will ...


6

Ordering is only significant in load balance/communication and suitability for preconditioning. The Krylov method does not care about the order or even whether the matrix entries are stored. In practice, a bad ordering may require much more communication than otherwise necessary when multiplying by a matrix. See the section of the PETSc User's Manual on "...


6

Does your dense Hamiltonian have structure (e.g., sparse plus low-rank) or is it unstructured dense (no compressed representation available)? Is it well-conditioned? You can use Elemental, either on its own or through PETSc's interface, for the parallel dense linear algebra. I'm not aware of a suitable library implementation of the exponential, but you ...


6

You can configure with separate libraries (--with-single-library=0) and link only the ones you need (e.g., -lpetscmat -lpetscvec -lpetscsys), but this is generally a waste of effort. If you use static libraries, then only the parts you reference go into the binary (if you're trying to squeeze the last megabyte out of a memory-constrained environment). PETSc ...


6

I would recommend looking at a time-dependent example because PETSc can provide a lot more diagnostics if you formulate at that level. For example, you could use a Rosenbrock method, an additive Runge-Kutta IMEX method, or others. Some involve Newton iteration, but that is not required and is not always the best approach. To use those methods, you can ...


5

First, sparse direct is completely different from sparse iterative. You cannot reliably predict performance if you don't have a good understanding of what your code is doing. For sparse MatMult, MatSolve, MatSOR, and similar kernels, you have an arithmetic intensity of no more than 1 flop/4 bytes of memory bandwidth. Meanwhile, most recent multicore chips ...


5

I have used Thrust in my linked cluster expansion project. Depending on the situation, Thrust can perform as well as or better than a low level implementation that you roll yourself (in particular, the reduce kernel has been working quite well for me). However Thrust's generic nature and flexibility means it sometimes has to do a lot of extra copying, array ...


5

There's no need to use transfer to emulate void * in a modern Fortran code. Instead, just use the ISO_C_BINDING intrinsic module, which is supported by all mainstream Fortran compilers. This module makes it very simple to interface between Fortran and C, with some very minor caveats. One can use the C_LOC and C_FUNLOC functions to get C pointers to ...


5

When there is a PetscCopyMode parameter, the behavior is explicit. We would like to be explicit in every instance, but that would become very cumbersome in the interface. Create() gives back a new reference, whereas Get() returns a borrowed reference and should be given back using Restore(). The exceptions come up in Set() methods. They should all ...


5

Your idea of keeping track of which i,j interactions you have found can work, I think that's the "CS trick" that you and Stefano M are referring to. This amounts to constructing your sparse matrix in list of lists format. Not sure how much CS you have so I apologize if this is already known to you: in a linked list data structure, every entry stores a ...


5

If you specify your mesh as a DMPlex and your data layout as a PetscSection, then DMCreateMatrix() will give you the correctly preallocated matrix automatically. Here are PETSc examples for the Poisson Problem and Stokes Problem.


5

The development version of scipy has recently added expm_multiply which computes the action of the matrix exponential without explicitly computing the matrix exponential itself. This uses the algorithm in "Computing the Action of the Matrix Exponential..." http://eprints.ma.man.ac.uk/1591/ https://github.com/scipy/scipy/pull/2456 The implementation has ...


5

Some thoughts from someone who has worked a fair amount in compiled languages, and has done a tiny bit of FVM: Typically, if you have experience programming in C, you sketch out a high-level description (pseudocode) of what you would like to do. Then you look for libraries that might implement the data structures and capabilities you need for your high-...


5

You say that you want an MPI version. Then you need to study the literature, as the distributed memory variant of matrix-matrix product are not a simple parallellization of the sequential version. The Cannon algorithm is pretty cute if you're on a square processor grid. In each step you rotate the input matrix rows and columns, so that in the end each ...


5

Your problem is too small. You have to consider that to get good efficiency, each processor has to have enough work to offset the cost of communication. In other words, there is a threshold how many degrees of freedom there have to be below which efficiency deteriorates. I don't know where that threshold is for MUMPS. But, to give just one example, for ...


4

In general, the Set methods (and PETSc) pass-by-reference, so it is not safe to destroy the object that you have passed to the method. Methods such as PCASMSetLocalSubdomains are the exceptions, not the rules.


4

I would say the canonical choice for this problem would be the Conjugate Gradient solver plus an algebraic multigrid preconditioner. For PETSc, hypre/boomeramg or ML would be the obvious preconditioner choices. All of these components when used through PETSc scale very well to thousands or tens of thousands of processors if the problem is large enough (at ...


4

If "make test" reports something like Running test examples to verify correct installation Using PETSC_DIR=/home/nate/src/PETSc-3.3 and PETSC_ARCH=arch-linux2-c-opt C/C++ example src/snes/examples/tutorials/ex19 run successfully with 1 MPI process C/C++ example src/snes/examples/tutorials/ex19 run successfully with 2 MPI processes Fortran example src/snes/...


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