32

Although LAPACK has some incredibly optimized code, it can still be worth it to write your own version in a few cases. The most important reason (and the reason they make you do it in your course). It's a great learning experience. You will never fully understand something like the QR algorithm (which is actually incredibly complicated to get the details ...


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


8

@ThiysSteel covers a lot, here is another perspective I find important: Even if you have available excellent implementations of any algorithm you might need, you still need to understand some of the internals in order to make good use of your library. Some examples: Any linear-algebra library can compute an inverse matrix really fast and as accuratly as ...


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


7

To clarify @ThiysSteel's good answer: The point is not to attempt to surpass the optimization of code written by very experienced people, who'd wrangled with it for decades. The point is to acquaint yourself with the issue that were addressed in that code, to appreciate the "niceties", and not be naive about related matters in the future. (It's ...


7

I used to try to optimize code via using assembly language (as opposed to C). I had some clear success, where a real-time microphone array worked with assembly language, but it would not work at all in the C language. However, since that time, more than 20 years have passed. It's very difficult to achieve this sort of improvement nowadays. Compilers have ...


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

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

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


6

This might be tangent but i feel it worth adding. A common way for user code to beat libraries is to exploit structure the libraries are not aware of. Block diagonal matrices are easier to work in structured form then dense matrix that contain the complete matrix and lots of zeros. Cyclic boundary conditions in PDEs can be solved more easily than working on ...


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

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

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

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

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


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


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