# Tag Info

17

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.

10

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

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

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

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

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

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

Good implementations of the BLAS and LAPACK routines (most importantly the BLAS routines) can be much faster than naive straight forward implementations of the same functions. However, efficient implementations typically include optimizations that are very specific to the particular computer that you're running on. Even different models of processors from ...

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

Many of these libraries have C interfaces that swap the meaning of the ordering internally without swapping the data. Also, using C-style double pointers for 2-D matrices is probably the wrong choice in the first place (due to the double lookup), so you can make your data appear column-major by indexing into a linear C array or use libraries that work around ...

3

There's nothing in BLAS or LAPACK that does this directly. You could use the dot product function in the BLAS to take the dot product of your vector and a vector of all 1's. However, it's probably more efficient to just write your own loop to compute the sum- this operation doesn't allow for any cache reuse so most compilers should be able to optimize the ...

3

Just write the loop code and forget about BLAS here. Adding two matrices is bandwidth-limited and the compiler will likely do as good a job as the BLAS implementation in this case.

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