17

If you're doing celestial mechanics over long time scales, using a classical Runge-Kutta integrator will not preserve energy. In that case, using a symplectic integrator would probably be better. Boost.odeint also implements a 4th-order symplectic Runge-Kutta scheme that would work better for long time intervals. GSL does not implement any symplectic methods,...


9

FTensor is a lightweight, header only, fully templated library that includes ergonomic summation notation. It has been tested extensively in 2, 3, and 4 dimensions, but should work fine for any number of dimensions.


9

Both GNU Scientific Library (GSL) (C) and Boost Odeint (C++) feature 8th order Runge-Kutta methods. Both are opensource, and under linux and mac they should be directly available from the package manager. Under windows, it will probably be easier for you to use Boost rather than GSL. GSL is published under the GPL license, and Boost Odeint under the ...


8

Optim.jl from Julia will work with the number types that you give it, so if you make it use BigFloats then it'll do that. Local derivative based, derivative-free, global, and integrates with automatic differentiation. From Julia, it's just: using Optim rosenbrock(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2 result = optimize(rosenbrock, big.(zeros(2)), ...


7

While there are some books on the implementation of finite element methods (e.g., Mark Gockenbach's Understanding and Implementing the Finite Element Method and Yair Shapiro's Solving PDEs in C++: Numerical Methods in a Unified Object-Oriented Approach), I think it would indeed be most useful to look at (and compare!) actual large-scale libraries: There are ...


7

Prior answers to this question have covered most of the salient points, but I want to add one comment with respect to this: does MKL have the upper hand for some tasks? The MKL team is in a unique position to know about future Intel instruction sets and their implementations in specific processors. Furthermore, they have access to proprietary processor ...


7

Almost everything you can build and install in your own space. With GNU autotools, you can do something like ./configure --prefix=/path/to/your/work/space ... and then follow the usual compilation instructions. Things based on CMake and Scons have similar facilities.


7

You can calculate GFLOP rates this way, but the numbers are pretty meaningless on today's hardware: Floating point operations require a variable number of clock cycles. An addition is generally cheaper than a multiplication, but each generally takes more than one clock cycle of the 2.8 billion cycles you quite. When you have hyperthreading, you have two ...


7

I'm the primary author of many Julia libraries geared toward "architecture-specific optimizations", including LoopVectorization.jl and Octavian.jl. For BLAS-like operations, one of the most important optimizations LoopVectorization.jl does is "register tiling". While CPUs may have a huge number of actual registers (used for register ...


6

XTensor is a modern approach and is getting more and more popular. https://github.com/QuantStack/xtensor


6

I think this new taco lib is really good too. The Tensor Algebra Compiler (taco) is a C++ library that computes tensor algebra expressions on sparse and dense tensors. It uses novel compiler techniques to get performance competitive with hand-optimized kernels in widely used libraries for both sparse tensor algebra and sparse linear algebra. You can use ...


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

For what its worth, Eigen does have a Tensor class as an unsupported module. http://eigen.tuxfamily.org/dox-devel/unsupported/group_CXX11_Tensor__Module.html I haven't used it myself so can't say more about it. The Armadillo class library has a 3rd-order tensor class. http://arma.sourceforge.net/ I haven't used the tensor capabilities of Armadillo ...


6

The SINTEF Matlab Reservoir Simulation Toolbox includes a GPL-licensed AD library. The usage is mostly geared towards numerical applications in subsurface flow, but the library itself is usable for more general purposes. Here is a basic runthrough of your example as you would run it from the base directory of MRST: startup; % Load ad based module ...


6

Of course it makes sense to use the GSL (or another library for that matter) for several reasons: Don't reinvent the wheel. The work has been done, you can spend your time on more useful things. If you do decide to implement these basic things yourself, the risk that your code will probably contain some bugs and will be slower, less memory efficient etc ...


6

I would like to hear comments from users that have some practical models (e.g. black-box hyperparameter optimization) which are still needed to be solved acceptably - whether this method works or not for their models, possibly with the description of the model. Looks like you want somebody to invest what may be considerable time and energy in trying out ...


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

The deal.II library (http://www.dealii.org), while written for much larger purposes, also has a sub-library of tensor classes that likely does a lot of what you want to do. In particular, it uses templates for the dimension. (Disclaimer: I am one of the principal authors of this library.)


5

I recommend looking at different library docs: http://fenicsproject.org/book/ http://www.dealii.org/developer/


5

Ask your remote host to install what you need. We do this all the time for folks where I work. Typically they can help you out. Also, it's OK to do some of your development on your machine and then port your results to the remote machine.


5

I am a PETSc developer so take my suggestion with a grain of salt, but I would use PETSc because the problem sizes are large enough that execution overhead should be minimal, you can trivially switch between various sparse and dense solvers (1/2 sparse should be treated as dense, but it might pay off to use a sparse solver for 1/9 at your sizes), a suite of ...


5

What do people suggest for the linear system of this type? I know about trilinos, petsc, and sundials, but don't know the other alternatives or have exposure to them. Generally speaking PETSc, Trilinos, and KINSOL (from SUNDIALS) are the best-of-breed when it comes to scalability. From an ecosystem standpoint, PETSc seems most flexible, since it does not ...


5

summarizing some points: If it's a long-term integration of a non-dissapative model, a symplectic integrator is what you're looking for. Otherwise, since it's an equation of motion, Runge-Kutta Nystrom methods will be more efficient than a transformation to a first order system. There are high order RKN methods due to DP. There are some implementations, ...


5

Daniel Shapero's answer is excellent, but I felt I should add the following: Properly preconditioned iterative methods will almost certainly win here, unless your systems have a very, very special structure. There's a considerable recent literature, unfortunately mostly of the "asymptotic a priori bounds on running time" variety, on solvers for graph ...


5

You may not be able to express the weak form in deal.II as a mathematical formula, but you come pretty close. For elasticity, the bilinear form reads $$ a({\mathbf \varphi}_i, {\mathbf \varphi}_j) = \left( \lambda \nabla\cdot {\mathbf \varphi}_i, \nabla\cdot {\mathbf \varphi}_j \right)_\Omega + \left( \mu \nabla\mathbf \varphi_i, \nabla\...


5

According to [1]: "However, the IEEE-754 standard specifies nothing for elementary functions" and "Indeed, the mathematical libraries (libm) provided by operating systems do not guarantee correct rounding.". As such, according to [1], the answer is that the IEEE-754 standard does not require the exp2 function to be correctly rounded, it only recommends it. [...


5

Chris code doesn't work on my machine, so here is my solution. import julia #julia.install() #<- this is probably needed for the first time from julia.api import Julia jl = Julia(compiled_modules=False) #from julia import Base from julia import Optim pyrosen="((x,y),)->(1.0 - x)^2 + 100.0 * (y - x^2)^2; " jl.eval("setprecision(512)")...


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

I want to add something not mentioned yet here that might be useful (too long for a comment), because you say that syncing with git is a pain. There is the tool sshfs (and I believe its osx equivalent is osxfuse). If your remote machine allows ssh connections, sshfs will be able to mount the remote filesystem locally over ssh. This is useful if, as I think ...


4

This article in SIAM Review describes how to implement operator-overloading automatic differentiation in MATLAB, and gives a good introduction to the technology.


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