# Tag Info

15

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

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

9

I can't answer the second half of your question as far as other implementations out there but I can provide some insight as to the challenges. For reference, I personally used ViennaCL on a nVidia GTX 560 Ti with 2GB of memory for my benchmarks. Over serial code on a mid-range i5, I saw speed-ups for dense matrix multiplications of approximately 40x. For ...

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

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

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

6

That is a classic example of cache associativity. The stride associated with that problem size is filling up certain sets causing cache eviction despite there being lots of space in other sets. Figure from Gustavo Duarte's excellent blog post on the topic See also Drepper's What every programmer should know about memory.

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

Unfortunately the desired routine is not part of BLAS but is very similar to zherk, which performs a Hermitian rank-k update. In particular, zherk supports operations which look like $$A := \alpha B B^H + \beta A,$$ where only one triangle of $A$ is updated (the traditional variable names are different in the zherk prototype, where usually $C$ is the ...

5

OpenBLAS is designed as a library to perform effectively for the most common computationally challenging problems that arise in scientific computing, linear algebra problems that take minutes or longer to solve. As @Stefano and @JedBrown hint in the comments, you are not necessarily likely to depend heavily on the performance of your BLAS library for the ...

5

Let me focus only on CUDA and BLAS. Speedup over an host BLAS implementation is not a good metric to assess throughput, since it depends on too many factors, although I agree that speedup is usually what one cares about. If you look at the benchmarks published by NVIDIA and take into account that the Tesla M2090 has 1331 Gigaflops (single precision) and ...

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

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

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

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

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

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

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