Tag Info

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

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

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

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

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

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

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.

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

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

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

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

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

I'm not completely clear on what you're asking here, but it looks like a solution to your problem might be the new Neighborhood Collective in MPI-3.0 Chapter 7.6 (you can find the PDF here). I'm far from an expert in these things, but the basic idea is that you attach a topology to your MPI_Communicator and get a new communicator where you can use sparse ...

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

4

PETSc is an outstanding linear algebra library. It also has modules for nonlinear solvers and, to some degree, for meshes and finite elements, but the latter aren't frequently used as there are much more comprehensive finite element packages out there. So, if you're interested in FEM calculations, my recommendation would be to go with one of the other ...

4

You can either do it from the code by doing this: KSP my_solver ; // define the KSP method PC my_prec ; // define the preconditioner /* Define the matrices to be used ... */ KSPCreate(comm,&my_solver); KSPSetType(my_solver,KSPBCGS); // Sets BiCGStab as the krylov method KSPGetPC(my_solver,&my_prec); PCSetType(my_prec,PCSOR); // Sets ...

4

The PETSc team always recommends that their users control solver options from the command line. The whole package is built with the idea of extreme flexibility in composing solvers and preconditioners, and the only way to achieve that is to use the command line scheme. You could get rid of setting the KSP and PC types in your code and use ierr = KSPCreate(...

4

No. Like many other libraries, you will spend far more time trying to isolate a piece than to simply use the whole thing. Why not just install PETSc somewhere and use it as-is?

4

We simply roll our own code in deal.II -- in essence, we tell the framework to execute tests using mpirun -np .... We had previously just used a Makefile-based testing scheme (compile, link, execute test, then compare the output with one that had previously been saved) and you can find this here: https://svn.dealii.org/branches/releases/Branch-8-0/tests/mpi/...

4

There are several MPI-enabled software packages that use the CMake set of tools for testing. The ones that I can think of off the top of my head are Trilinos, VTK and ParaView. I would think that you don't want to assume that the executable needs to be launched with mpirun and/or mpiexec. CMake has support for specifying how to properly launch the executable ...

4

Assuming your overall matrix is positive definite (it is definitely symmetric), then I would suggest looking into algebraic multigrid (AMG) methods as preconditioners. They compute hierarchies of sparsified matrices themselves. If you're already using PETSc, take a look at the hypre preconditioner. Using this may force you to actually multiply out the ...

4

If you are only interested in the smallest eigenvalue, the conjugate gradient method applied to the matrix $L$ gives you a good approximation after a reasonably small number of steps, and you won't have to solve any linear systems. The details are in Y. Saad's book on iterative methods, but here a short summary: From the coefficients that are computed in ...

4

Heh, fun question. If you read the C standard carefully, you'll find wording like (C99, 6.5.6.8) about pointer arithmetic. If both the pointer operand and the result point to elements of the same array object, or one past the last element of the array object, the evaluation shall not produce an overflow; otherwise, the behavior is undefined. As I ...

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