# Tag Info

45

I used to implement everything myself, but lately have begun using libraries much more. I think there are several very important advantages of using a library, beyond just the issue of whether you have to write a routine yourself or not. If you use a library, you get Code that has been tested by hundreds/thousands/more users Code that will continue to be ...

37

The difference in your timings seems to be due to the manual unrolling of the unit-stride Fortran daxpy. The following timings are on a 2.67 GHz Xeon X5650, using the command ./test 1000000 10000 Intel 11.1 compilers Fortran with manual unrolling: 8.7 sec Fortran w/o manual unrolling: 5.8 sec C w/o manual unrolling: 5.8 sec GNU 4.1.2 compilers Fortran ...

34

There is substantial programmer overhead involved in linking to a library function, especially if that library is new to the programmer. It is often simpler to just rewrite simple algorithms rather than figure out the specifics of a particular library. As the algorithms become more complex this behavior switches. Python has excelled at reducing this ...

19

One of the projects I'm involved in right now is writing a flexible simulation and analysis package for a class of particle physics detectors. One of the goals of this project is to provide the code base to be used in these things for decades to come. At this point point we already have two dozen dependencies, making the build process such a nightmare that ...

19

Possibly one could start with the function $\mathtt{expm1}$ which is part of the C99 standard, and calculates $e^x-1$ accurately near $x=0$.

17

This is an instance of cancellation error. The C standard library (as of C99) includes a function called expm1 that avoids this problem. If you use expm1(x) / x instead of (exp(x) - 1.0) / x, you won't experience this issue (see graph below). The details and solution of this particular problem are discussed at length in Section 1.14.1 of Accuracy and ...

16

I'm coming late to this party, so it's hard for me to follow the back-and-forth from all above. The question is big, and I think if you are interested it could be broken up into smaller pieces. One thing I got interested in was simply the performance of your daxpy variants, and whether Fortran is slower than C on this very simple code. Running both on my ...

15

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

12

I think it is quite common, with some algorithms more likely to be re-implemented than others. There's a tricky trade-off between how annoying a library is to install, how hard it is to implement the algorithm yourself, how hard it is to optimize it, and how well the library fits your needs. Also, sometimes using a library is just overkill: I used the slow ...

12

First of all, thanks for posting this question/challenge! As a disclaimer, I'm a native C programmer with some Fortran experience, and feel most at home in C, so as such, I will focus only on improving the C version. I invite all Fortran hacks to have their go too! Just to remind newcomers about what this is about: The basic premise in this thread was that ...

11

Such an effect happens because of how the data of the int** a is stored in memory (as per C/C++). This question on StackOverflow has answers with some more details (in particular, a difference between int** and int[][] that many users noted in the comments), and how it looks like an array of arrays - it's just laid out contiguously in memory. It's worth ...

10

I think that implementing an algorithm instead of using a library can sometimes give a better understanding and control of the model. When I am coding some program for scientific computations, it's important for me to understand what I am doing. Implementing the important algorithms helps me to get a better knowledge of the problem and achieve better control ...

10

You might consider it "overkill", but PETSc's time integration package can be used with C99 complex (configure --with-scalar-type=complex). Supported methods include explicit Runge-Kutta low-memory strong stability-preserving Runge-Kutta Rosenbrock-W additive Runge-Kutta IMEX These implementations are most appropriate for high-dimensional problems such as ...

10

With the Intel compiler of any modern vintage, -O3 -vec-report3. Optimization level three guarantees that it's trying to vectorize, and the vector report will tell you what it's doing. The GNU page on vectorization says that it's on by default at optimization level 3, but I can't find the equivalent of vec-report.

9

The way I would write AXPY in Fortran is slightly different. It is the exact translation of the math. m_blas.f90 module blas interface axpy module procedure saxpy,daxpy end interface contains subroutine daxpy(x,y,a) implicit none real(8) :: x(:),y(:),a y=a*x+y end subroutine daxpy subroutine saxpy(x,y,a) ...

9

Java has been around for almost 20 years now as a major programming language, but it hasn't caught on in scientific computing so far. I think that's a good indicator for what's going to happen in the future. My take is that the issue isn't speed. Most people are probably willing to give up 20% of performance (or even a factor of 2) if they would be vastly ...

9

There is a prime missing both in your code and in the expression in your question. In the original paper, the expression is: $$\log(x/x_0) + (x-x_0) \frac{\sum^\prime_{0\leq j\leq n}p_j T_j^*(x/6)}{\sum^\prime_{0\leq j\leq n}q_j T_j^*(x/6)}.$$ This primed sum is very common and standard when dealing with Chebyshev series: it means that the first term of ...

8

One answer is that there are so many slight variations to numerical code that it is really hard to encapsulate that in a library. Take this in comparison to web software, which is often easy to install and has a clear set of inputs and outputs. I think more common is people grabbing a framework, or big library that acts like a framework (Trilinos/PETSc), ...

8

There is a C++ library which is quite mature. This is probably as close as you will get to C. I myself haven't found any usable C library yet. You could use the C++ library and still write most of your code in C using extern C { } in the C++ code.

8

Within the GNU compiler collection, you have the option -ftree-vectorizer-verbose=n where n is a number between 0 and 6 which will print information similar to icc/ifort.

8

Since you're already using C++ and your matrices are symmetric positive definite, I would perform an unpivoted $LDL^T$ factorization of $Q$ and also of $12I-Q-J$. Here I'm assuming that $12I-Q-J$ is also positive definite, otherwise the $LDL^T$ will require pivoting for numerical stability (it's also possible that even though it's not positive definite, ...

8

I think some of your issues are more important than others and some of your emphasis is misplaced. In pursuing overhead, you are in danger of making your program unmaintainable. It is easier to write a common program and direct surplus effort somewhere more interesting. I apologize for pontificating like this. If statements. From a strict programming ...

8

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

In theory, as the original authors, you're free to pick and name a standard, then expect others to follow it. In practise, if you're supporting an HPC system, then your choice is likely to be restricted among the compiler standards that the given system you're using for testing has for its tool stack (you are testing your software, aren't you?). This might ...

7

Before deciding whether or not to use libraries, I think you'd also want to figure out how much the use of a library will help your code. If you're going to be using a well-optimized library for a key computational kernel, then it's probably a lot more efficient than trying to write your own. However, if you're writing a specialized routine that's only ...

7

A simple and accurate approach in this case is to use standard differences in $\log$ space. \begin{align*} s(x) &= \log x \\ \partial_x f(x) &= \partial_s f(e^{s(x)}) \partial_x s(x) \\ &= \frac 1 x \partial_s f(e^{s(x)}) \\ \end{align*} Since $s(x)$ is equally spaced, it is easy to approximate the derivative on the bottom line using ...

7

I've had to deal with similar problems before, and my favourite solution is to use Memory-mapped I/O, albeit in C... The principle behind it is quite simple: instead of opening a file and reading from it, you load it directly to the memory and access it as if it were a huge array. The trick that makes it efficient is that the operating system doesn't ...

7

user389's answer has been deleted but let me state that I'm firmly in his camp: I fail to see what we learn by comparing micro-benchmarks in different languages. It doesn't come as much of a surprise to me that C and Fortran get pretty much the same performance on this benchmark given how short it is. But the benchmark is also boring since it can easily be ...

7

As Aron notes in his comment, rarely are matrices stored in nested pointer data structures in practice. Multiple levels of indirection (I'm told) cause considerable performance penalties, and require multiple allocations and frees (more performance issues, plus more potential to screw up and segfault). I've seen nested pointers used to store matrices only a ...

7

I would suggest to exactly duplicate the Lapack interface to the function that you need, most probably you just need dgesv. That way people that have Lapack installed can simply link to it and it will just work. For people that don't have Lapack installed, you provide your own simple implementation of this function, or possibly implement it using Eigen or ...

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