People who code: we want your input. Take the Survey

Hot answers tagged

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

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

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


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

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


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

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

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

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


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

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

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

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

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


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

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

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

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

Continuous weak form Though I think the weak form is more fundamental, suppose we start with the strong divergence-form representation for a (first- or) second-order quasilinear PDE: find the $m$-component solution $u \in R^m(\Omega)$ $$ -\nabla \cdot f_1(u,\nabla u) + f_0(u,\nabla u) = 0 $$ on the domain $\Omega \subset R^d$ with (for simplicity) ...


4

I like this formulation for implementing Dirichlet boundary conditions in cases where elimination of boundary dofs is not convenient. If you apply effectively the same procedure to $M$, the rows and columns corresponding to Dirichlet boundary dofs will be rows and columns of the identity. This will give you new eigenvalues equal to 1 (or $\alpha$ in the ...


4

If you're using iterative methods with assembled matrices, just buy DDR channels. Don't pay attention to number of cores when loooking at the spec sheet. Within the same class of processors (e.g., a recent generation of Xeon), the achievable memory bandwidth will be proportional to the quoted peak bandwidth. Note that for very small problem sizes, you might ...


4

Unfortunately there's no tool for this. You can run each on a variety of input sizes to establish the computational complexity they appear to have, i.e. the $f$ in the $O(f(n))$ that characterizes each code. This can point you towards what underlying algorithms each is using and verify or not that something that should be $O(n)$ is actually implemented to ...


4

It is indeed possible to write an obstacle problem for an advection-diffusion equation as a variational inequality: If $a(u,v)$ is the bilinear form corresponding to your advection-diffusion equation, the corresponding obstacle problem (in your specific case) is finding $u\geq 0$ such that $$ a(u,v-u) \geq (f,v-u) \qquad \text{for all } v\geq 0.$$ In ...


4

I think you have some confusion. PETSc is not in the same league as Fenics, Libmesh, Moose etc. In fact, all of these (heavyweight) packages use PETSc for linear algebra. IMHO PETSc is as lightweight as you can get. It just requires C/Fortran compilers and Python (used only for configuration) and you can build the library in under 5 minutes on your laptop. ...


4

I can recommend nutils. nutils meets at least a few your "light-weight" requirements. it is pure python and easy to install since it only depends on standard Python libraries numpy, scipy, and matplotlib and, thus, it is well suited for interoperations. At least the developers claim that "Exposed objects are of native python type or allow for easy ...


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