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

So there is a ton to say about this, and we will actually be putting a paper out that tries to summarize it a bit, but let me narrow it down to something that can be put into a quick StackOverflow post. I will make one statement really early and keep repeating it: you cannot untangle the efficiency of a method from the efficiency of a software. The details ...


19

See https://math.stackexchange.com/questions/861674/decompose-a-2d-arbitrary-transform-into-only-scaling-and-rotation (sorry, I would have put that in a comment but I've registered just to post this so I can't post comments yet). But since I'm writing it as an answer, I'll also write the method: $$E=\frac{m_{00}+m_{11}}{2}; F=\frac{m_{00}-m_{11}}{2}; G=\...


17

As Jed Brown mentioned, the connection between gradient descent in nonlinear optimization and time stepping of dynamical systems is rediscovered with some frequency (understandably, since it's a very satisfying connection to the mathematical mind since it links two seemingly different fields). However, it rarely turns out to be a useful connection, ...


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


16

This issue has become much more nuanced as the changes in architectures has shifted the HPC landscape. As Wolfgang Bangerth mentions one current longstanding view, I'll split my answer into basic definitiions and further details. Basic Definition A node refers to the physical box, i.e. cpu sockets with north/south switches connecting memory systems and ...


14

You can't beat an explicit formula. You can write down the formulas for the solution $x=A^{-1}b$ on a piece of paper. Let the compiler optimize things for you. Any other method will almost inevitably have if statements or for loops (e.g., for iterative methods) that will make your code slower than any straight line code.


13

Steve McConnell's book Code Complete, 2nd edition has an extensive bibliography discussing these issues from more of the standpoint of software developers than computational scientists. The book is starting to become a little dated, in that it's approaching a decade old, so it doesn't cover more recent testing methodologies like behavior-driven development. ...


13

While I haven't seen the exact formulation that you have written down here, I keep seeing talks in which people "rediscover" a connection to integrating some transient system, and proceed to write down an algorithm that is algebraically-equilavent to one form or another of an existing gradient descent or Newton-like method, and fail to cite anyone else. I ...


12

Two examples of libraries that use modern C++ constructs: Both the eigen and armadillo libraries (linear algebra) use several modern C++ constructs. For instance, they use both expression templates to simplify arithmetic expressions and can sometimes eliminate some temporaries: http://eigen.tuxfamily.org http://arma.sourceforge.net/ http://hpac.rwth-...


12

Yes, this is the standard Aubin-Nitsche (or duality) trick. The idea is to use the fact that $L^2$ is its own dual space to write the $L^2$-norm as an operator norm $$\|u\|_{L^2} = \sup_{\phi\in L^2\setminus\{0\}} \frac{(u,\phi)}{\|\phi\|_{L^2}}.$$ We thus have to estimate $(u-u_h,\phi)$ for arbitrary $\phi\in L^2$. To do that, we "lift" $u-u_h$ to $H^1_0$ ...


12

Though it is a relatively rare situation when you actually have to calculate an inverse of the matrix, not all techniques were created equally. I would use the term badly-conditioned instead of ill-conditioned. For badly conditioned matrices, you might opt in the SVD-route to calculate the inverse: $$ A=U\Sigma V^H \implies A^{-1}=V\Sigma^{-1}U^H. $$ If ...


11

I posted some benchmarks here: http://www.cecm.sfu.ca/~rpearcea/mgb.html These are for total degree orders. To solve systems you typically need to do more work. Timings are for a typical midrange desktop as of 2015 (Haswell Core i5 quad core). The fastest system on one core is Magma, which uses floating point arithmetic and SSE/AVX. Magma is the ...


10

@Pedro Gimeno "I doubt it can be any more robust than that." Challenge accepted. I noticed the usual approach is to use trig functions like atan2. Intuitively, there shouldn't be a need to use trig functions. Indeed, all the results end up as sines and cosines of arctans--which can be simplified to algebraic functions. It took quite a while, but I managed ...


10

I would suggest taking a look at Deal.II. It uses the STL, it's own iterators, shared pointers, etc. The various linear solvers can use the various matrices because of how it was designed. I haven't come across any use of move semantics, but that doesn't mean they aren't there. Here is a link.


10

For parabolic/elliptic PDE's, I highly recommend Beatrice Riviere's book: Discontinuous Galerkin methods for solving elliptic and parabolic equations: theory and implementation. For hyperbolic PDE's and general (i.e. nonlinear) conservation laws, I recommend Hesthaven & Warburton's book: Nodal discontinuous Galerkin methods: algorithms, analysis, and ...


10

I haven't worked in quantum chemistry specifically, but I've worked in other areas where high performance is a correctness requirement (along with scientific accuracy), so I think we're on the same page here. Broad but shallow knowledge of all of the above is absolutely necessary for the team as a whole. Deep knowledge can be acquired as needed, or hired as ...


9

From skimming the table of contents to the book you listed, I'd say that computational books of that type for physics (or in my case, engineering, since that is my background) tend to sacrifice depth and quality of explanation for breadth. The best book of this ilk I can think of is probably Strang's Computational Science and Engineering, because he's a ...


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


9

Since the matrix is so close to the identity, the following Neumann series will converge very rapidly: $$A^{-1} = \sum_{k=0}^\infty (I-A)^k$$ Depending on the accuracy required it might even be good enough to truncate after 2 terms: $$A^{-1} \approx I + (I - A) = 2I - A.$$ This might be slightly faster than a direct formula (as suggested in Wolfgang ...


9

The term inverse crime for a numerical test of a parameter identification method that uses data contained in the range of the discrete(!) forward operator used for the inversion (thus essentially reducing the problem to a well-posed finite-dimensional one that behaves fundamentally different from the original infinite-dimensional one -- it is important to ...


8

Kate, you asked for a reference request as a high school student interested in getting started in computational science, which I think I can tackle fairly specifically. As long as you are comfortable (or getting comfortable) with Calculus, I think there are two great self-interest textbooks for you to go, depending on your interest and access to Python and ...


8

The GSL has a 2-by-2 SVD solver underlying the QR decomposition part of the main SVD algorithm for gsl_linalg_SV_decomp. See the svdstep.c file and look for the svd2 function. The function has a few special cases, isn't exactly trivial, and looks to be doing several things to be numerically careful (e.g., using hypot to avoid overflows).


8

I have to disagree with Paul: Functional analysis is a beautiful, elegant topic with an enormous range of applications (but de gustibus... :)). But in any case, you don't need to know a lot of (pure) functional analysis to understand finite element methods: Normed vector spaces, (strong) convergence, completeness, inner products and their induced norms, ...


8

I think you start by looking at something like FEniCS. Marie Rognes has a presentation with code examples and a paper discussing the theory and implementation. libMesh is supposed to be able to do something similar for 2-manifolds in 3-space, and so is deal.II, judging from this manuscript. Developers of deal.II and FEniCS answer questions on SciComp, and ...


8

Googling benchmark polynomial systems leads to a few hits, including the University of Mannheim's Computer Algebra Benchmark Initiative. Sadly, most of these are out of date or defunct. The most active seems to be the SymbolicData Wiki, but as far as I can tell, it only collects benchmark problems, not benchmark results. Some comparisons (dating back to ...


8

One issue that you should be aware of is that NVIDIA has a market segmentation strategy in which it sells relatively inexpensive GPU's to the gaming and graphics workstation markets (GeForce and Quadro) and different higher-priced models (Tesla) to the high-performance computing market. The GPU's sold for use in gaming and graphics have limited double-...


7

When we say "numerically robust" we usually mean an algorithm in which we do things like pivoting to avoid error propagation. However, for a 2x2 matrix, you can write the result down in terms of explicit formulas -- i.e., write down formulas for the SVD elements that state the result only in terms of the inputs, rather than in terms of intermediate values ...


7

I needed an algorithm that has little branching (hopefully CMOVs) no trigonometric function calls high numerical accuracy even with 32 bit floats We want to calculate $c_1, s_1, c_2, s_2, \sigma_1$ and $\sigma_2$ as follows: $A = USV$, which can be expanded like: $ \begin{bmatrix} a & b \\ c & d \end{bmatrix} = \begin{bmatrix} c_1 & s_1 \\ -...


7

Back in the era of Intel 387 math coprocessors I had to maintain an interrupt handler for floating point exceptions. Apart from that, I agree that pretty much everyone ignores denormals (or subnormals as they are sometimes called), and this is successful in part because the IEEE 754 standard default handling of "gradual underflow" hides their (presumably ...


7

It's worth a bit of a discussion as to why denormals don't matter in typical scientific computation. Usually, we are given some input $x$ and want to compute a function $f(x)$. Due to numerical errors, we instead compute $\tilde{f}(x) \approx f(x)$. There are two kinds of accuracy typically sought: Relative: $$|f(x)-\tilde{f}(x)| < \alpha |f(x)|$$ ...


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