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

Hot answers tagged

16

As far as I know, Lapack is the only publicly available implementation of a number of algorithms (nonsymmetric dense eigensolver, pseudo-quadratic time symmetric eigensolver, fast Jacobi SVD). Most libraries that don't rely on BLAS+Lapack tend to support very primitive operations like matrix multiplication, LU factorization, and QR decomposition. Lapack ...


15

No, that is not guaranteed. If you are using a NETLIB BLAS without any optimizations, it it mostly true that the results are the same. But for any practical usage of BLAS and LAPACK one uses a highly optimized an parallel BLAS. The parallelization causes, even if it only works in parallel inside the vector registers of a CPU, that the order how the single ...


12

Consolidating the comments: No, you are very unlikely to beat a typical BLAS library such as Intel's MKL, AMD's Math Core Library, or OpenBLAS.1 These not only use vectorization, but also (at least for the major functions) use kernels that are hand-written in architecture-specific assembly language in order to optimally exploit available vector extensions (...


10

Have you looked at GNU Scientific Library's implementation? I find the source code to be sufficiently readable and the routines are well documented.


9

This can be interpreted as summing over an index of a tensor when the vector $x$ is reshaped into a box of numbers instead of a list. In particular, if $X$ is the $d\text{-by-}d$ folded version of $x$, then the operation you are doing is, \begin{align} Dx &= \mathrm{vec}\left((I \otimes \mathbf{1})\mathrm{vec}(X)\right) \\ &= \mathrm{vec}(\mathbf{1}^...


9

tl;dr Use loops My numbers indicate that ifort is smart enough to recognize the loop, forall, and do concurrent identically and achieves what I'd expect to be about 'peak' in each of those cases. gfortran, on the other hand, does a bad job (10x or more slower) with forall and do concurrent, especially as N gets large. Both ifort and gfortran seem to ...


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

The Short Answer If the two BLAS implementations are written to carry out the operations in the exact same order, and the libraries were compiled using the same compiler flags and with the same compiler, then they'll give you the same result. Floating point arithmetic is not random, so two identical implementations will give identical results. However, ...


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

I think you are overestimating the overhead of computing L. There are zero extra operations needed; the only additional cost is writing to RAM some numbers that you have already computed anyway. The algorithms commonly used (in Lapack, for instance) to compute U also compute L along the way, and you'd save 0 flops by omitting it. For instance, if you think ...


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

A notable, C language implementation of BLAS is ATLAS. Among useful features: Algebra routines implemented both as straightforward C as well as highly-optimized assembler assisted versions for multiple architectures and variants. The build system features an "auto-tuner" which compiles multiple variants of the ATLAS library to establish which one will be ...


6

BLAS routines do not typically use stable summation algorithms. In the case of gsl, you can look up its source code online - the source of gsl's sdot is contained in gsl/cblas/source_dot_r.h, and contains this loop: for (i = 0; i < N; i++) { r += X[ix] * Y[iy]; ix += incX; iy += incY; } It's just a straightforward sum. The corresponding ...


5

I can't give you a good answer for 1, but I can give you decent answers for 2 and 3. I'm not terribly familiar with the Boost interface to UMFPACK. In the C interface, normally, you call UMFPACK routines that will allocate memory. If the memory cannot be allocated because there is not enough free memory available, UMFPACK will return a null pointer. You can ...


5

There's nothing surprising about these results. Matrix multiplication is well known to be communication intensive, and you've got a relatively slow communications network between your four nodes. Using MPI between two processes on the same node is certainly faster than using MPI between processes on different nodes because you don't have the bandwidth ...


5

In comparison with things like matrix vector multiplication (in which there's no cache reuse and everything has to come out of memory), matrix-matrix multiplication allows for lots of cache reuse in a careful implementation. Performance depends on having a good implementation of BLAS and perhaps depends on how much memory bandwidth is available although is ...


5

Yes, you want to call the BLAS routine DGEMM. The place to start for how to call it from C is to look at the documentation for DGEMM, which you can find online. Then you want to understand how to call FORTRAN routines from C (DGEMM, like all the standard BLAS routines has a FORTRAN calling convention). For example, this document https://computing.llnl.gov/...


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

The various Fortran standards allow a lot of compiler dependent behaviour in terms of function binary interfaces when being called with "complicated" data types such as Fortran90 style arrays and complex numbers. This means calling code compiled with one compiler from another one is not guaranteed to do what you expect, and can lead to grabbing the wrong bit ...


5

You should definitely jump to C99, or newer(!). The C99 standard introduced the restrict keyword. Loosely speaking, with this keyword you can inform the compiler that A[i] and B[j] do not access the same memory location. In that case the compiler can generate better optimized code. For example, it makes it easier for the compiler to auto-vectorize code. ...


4

Since basically all that xTRTRS does is call xTRSM (from BLAS), why not just call that directly? xTRSM contains a SIDE parameter which lets you specify exactly what you want.


4

Folks, I found this discussion very interesting, but I was surprised to see that re-ordering the loops in the Matmul example changed the picture. I don't have an intel compiler available on my current machine, so I am using gfortran, but a rewrite of the loops in the mm_test.f90 to call cpu_time(start) do r=1,runs mat_c=0.0d0 do j=1,n ...


4

I believe that MKL has the threaded parallel and serial functions in one unified library. You can try setting OMP_NUM_THREADS or MKL_NUM_THREADS to a range of values and see how the performance varies. Setting either to 1 will give you the serial behavior.


4

Among open-source BLAS, as far as I know, OpenBLAS (http://www.openblas.net/) is the best option. The website has a DGEMM benchmark, comparing against MKL (see below) and the reference Fortran BLAS. The library is threaded and written in C and assembly. For GPUs, there's clBLAS (https://github.com/clMathLibraries/clBLAS) that implements BLAS using OpenCL. ...


4

In general, no. Leaving associativity aside, the choice of compiler flags (for example, SIMD instructions being enabled, usage of fused multiply add, etc.) or the hardware (e.g., whether extended precision is being used) may produce different results. There are some efforts to get reproducible BLAS implementations. See ReproBLAS and ExBLAS for more ...


4

Some things I can think of: use sparse matrices for Matrix1 and Matrix2 to speed up the computations of dZand dY use larger integration tolerances reltoland abstol, especially if your are searching for steady-state solutions and/or do not need a precise resolution of the transient dynamics of your system. you are solving with ode15s, which is an implicit ...


3

Not well. If serial is a common case, it is important to wrap and drop down to lapack for serial execution. I implemented this in my code. For 2013 MKL pzhegvx with $n \approx 100 (1000)$ seems to incur 30% (100%) overhead compared to zhegvx when executed in serial. This seems high to me, so I'm a little worried about my implementation. Note that I ...


3

For matrices of that size, I'm not sure if you want to use ScaLAPACK at all. If you've got the ScaLAPACK code already, it shouldn't be hard to implement your own logic to drop into LAPACK instead. At the very least, doing that will allow you to perform the experiments required to answer your own question.


3

Geoff mentioned the limitations of sparse direct solvers. To solve a 3D problem of size $N = n^3$, an optimal direct solver will need $N^{4/3}$ memory and $N^2$ time. See this question for more details. But it sounds like you have a structured grid diffusion problem. If the diffusion coefficient is smooth (or constant), you can solve the problem using ...


3

There is some good information here: https://software.intel.com/en-us/mkl-linux-developer-guide-calling-intel-mkl-functions-from-multi-threaded-applications set MKL_NUM_THREADS=n call mkl_set_num_threads(n)


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