Hot answers tagged

15 votes
Accepted

Are BLAS implementations guaranteed to give the exact same result?

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 ...
13 votes
Accepted

Beating typical BLAS libraries matrix multiplication performance

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

Are BLAS implementations guaranteed to give the exact same result?

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 ...
  • 1,522
8 votes
Accepted

C standard for computational science

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 ...
  • 2,199
7 votes
Accepted

Getting to know about various BLAS implementations

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 ...
6 votes
Accepted

The difference between mkl_intel_lp64 vs mkl_gf_lp64 in a numerical reproducibility issue with Intel MKL

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 ...
  • 2,199
6 votes

Do BLAS routines compute their respective operations with minimum error?

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 ...
  • 11.4k
6 votes
Accepted

Functions from Scipy, Blas, or Lapack that compute only upper triangular matrix

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

BLAS, LAPACK or ATLAS for Matrix Multiplication in C

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 ...
  • 5,784
5 votes
Accepted

How does matrix-matrix product scale with multiple CPUs?

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

Distributed (MPI) matrix matrix multiplication

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

C standard for computational science

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 ...
  • 551
4 votes
Accepted

BLAS libraries for Octave or Matlab, preferrably with GPU support?

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

Are BLAS implementations guaranteed to give the exact same result?

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 ...
4 votes
Accepted

How to set up the differential equation system to speed up computation?

Some things I can think of: use sparse matrices for Matrix1 and Matrix2 to speed up the computations of ...
  • 1,450
3 votes

BLAS libraries for Octave or Matlab, preferrably with GPU support?

MATLAB already comes with Intel MKL for its BLAS implementation. There's no reason to replace that. As for using GPUs, if you make your array a gpuArray (to do ...
3 votes

Optimized parallel routine for $X' W X$ with $W$ diagonal

You haven't said anything about how much storage you have available to use for this computation. A 20 million by 500 array of double precision floats would require 80 gigabytes of RAM to store. If ...
3 votes
Accepted

Wrong result of 'ddot' from BLAS

The strange thing in your code is ddot_ being declared as extern C int, while it is actually a ...
  • 8,461
3 votes

Smart way to multiply 3 matrices

Have you considered working with the Cholesky (or low-rank) factor of $\rho(t_0)$ rather than with the matrix itself? This might reduce the number of products that you need to make, and it has the ...
3 votes

BLAS, LAPACK or ATLAS for Matrix Multiplication in C

By its anatomy the DGEMM ($C = \alpha AB + \beta C$) is one of the most optimizeable routines in computer science. For historic reason this routine is implemented in FORTRAN or the implementation ...
3 votes
Accepted

Impact of frequent row major / column major conversions using LAPACK/BLAS?

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 ...
  • 10.8k
3 votes

BLAS/LAPACK subroutine to add two matrices with different offsets and leading dimensions

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.
  • 2,041
3 votes

Danger of complex arithmetic in scientific computing

You say that the problem with complex arithmetic is that there are different ways to define the scalar product for complex vectors, compared to just one way in the real case. I think the real problem ...
3 votes
Accepted

BLAS/LAPACK Non absolute sum

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

Fast matrix multiplication with matrix elements computed on-the-fly (without forming the matrix)

If the matrix is small enough to fit into memory, then there is of course no cost associated with actually forming the elements: You will have to compute the elements at least once anyway to perform ...
3 votes

Parallel assembly of matrix

The first suggestion (without seeing the code) is to check if your memory access pattern is reasonable or, preferably, optimal for the critical parts of the code. One needs to keep in mind how the ...
  • 8,461
3 votes

Getting to know about various BLAS implementations

Abdullah, thank you for the plug for our materials. We have repackaged these as a Massive Open Online Course (MOOC) on edX titled "LAFF-On Programming for High Performance". It is free for ...
2 votes

Impact of frequent row major / column major conversions using LAPACK/BLAS?

You can use the identity C^T = (AB)^T = B^T A^T to use BLAS to compute on row-major arrays from C. In any case, the overhead of transposition is usually not a bottleneck compared to O(N^3) BLAS3 or ...
  • 2,041
2 votes
Accepted

Level 3 BLAS accelerated solver for banded linear systems.

Unfortunately, it looks like presently the only BLAS routines which take advantage of triangular/band structure are both level-2, _tbmv and ...
  • 136
2 votes

Fastest way to perform element-wise multiplication on a sparse matrix

If there is a choice in programming language, one option can be to use Julia, which has built in support for sparse matrices (via Suitsparse). The timing come out to about one and a half milliseconds ...
  • 121

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