# Tag Info

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

20

Your matrix is of size 15,000 x 15,000, so you have 225M elements in the matrix. This makes for roughly 2GB of memory. This is much more than the cache size of your processor, so it has to be loaded completely from main memory in every matrix multiplication, making for approximately 100GB of data transfers, plus what you need for the source and destination ...

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

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

PETSc uses BLAS for a few vector primitives, but these are generally limited by memory bandwidth and there isn't much variance in "optimization", so it tends not to make much performance difference. It also uses Lapack for some analysis such as Lanczos or Arnoldi estimates of eigenvalues and singular values, but these are generally not performance-sensitive....

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

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

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

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

If you have dense matrices with structure (e.g. fast transforms, Schur complements, etc), PETSc could be useful. In these cases, you won't be assembling the full matrix. For assembled dense systems, PETSc currently uses PLAPACK, but the matrix distribution in PETSc native format is not the best to minimize communication (for most operations). Jack Poulson, ...

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

If your goal is really to squeeze as much performance out as possible, then it is important to remember: The (BLAS) library might not have been tuned for your exact system/configuration. Library developers make mistakes. A vendor-tuned BLAS library should certainly be your default approach, but if you have taken the time to time individual kernels and ...

6

I think, it is not only interesting how a compiler optimizes the code for modern hardware. Especially between GNU C and GNU Fortran the code generation can be very different. So let's consider another example to show the differences between them. Using complex numbers, the GNU C compiler produces a large overhead for nearly very basic arithmetic ...

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

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

It's important to realize that parallel dense linear algebra libraries usually focus on level 3 BLAS routines (routines that perform $O(n^3)$ work with $O(n^2)$ data) and higher-level functionality like factorizations and eigensolvers. They usually don't tune the level 1 and level 2 operations that you're referring to. Since you mentioned that you are on a ...

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

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

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

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