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 ...
M.K. aka Grisu's user avatar
14 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 ...
9 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 ...
origimbo's user avatar
  • 2,249
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 ...
Tyler Olsen's user avatar
  • 1,522
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 ...
Chris Elrod's user avatar
7 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 ...
wim's user avatar
  • 571
7 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 ...
Federico Poloni's user avatar
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 ...
origimbo's user avatar
  • 2,249
5 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 ...
Laurent90's user avatar
  • 1,868
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. ...
Victor Eijkhout's user avatar
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. ...
Aditya Kashi's user avatar
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 ...
Juan M. Bello-Rivas's user avatar
4 votes
Accepted

Fast weighted vector inner product x*A*y with BLAS/LAPACK

The product $x^TAy$ can be calculated as: $\sum_{i}\sum_j x_i A_{i,j} y_j$. Easily implemented as a double for loop. You can additionally avoid a few flops by doing either: $\sum_{i}x_i \sum_j A_{i,j} ...
Thijs Steel's user avatar
  • 1,693
4 votes

Is it possible to express an arbitrary tensor contraction in terms of BLAS routines?

You can use reshape and pointwise multiply to reduce the operation to a matmul in terms of two temporary tensors $c_1,c_2$. Consider following transformations: pointwise mul to turn $a_i * b_i$ into $...
Yaroslav Bulatov's user avatar
4 votes
Accepted

What is the correct way of computing $LL^T$ in packed format with blas/lapack

The closest routine in BLAS for this would be DSYRK which performs the update C := alpha*A*A**T + beta*C, however ...
IPribec's user avatar
  • 607
4 votes

Weird runtime behavior of `scipy.linalg.solve_triangular` and `trtrs`

Thanks to Federico Poloni's comment, I have revised the script as below. It looks like scipy.linalg.solve_triangular has time complexity $O(p^2)$, because the ...
nalzok's user avatar
  • 181
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 ...
Robert van de Geijn's user avatar
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 ...
Anton Menshov's user avatar
  • 8,672
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 ...
Wolfgang Bangerth's user avatar
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 ...
Anton Menshov's user avatar
  • 8,672
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 ...
Chris Rackauckas's user avatar
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 ...
H. Rittich's user avatar
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 ...
Federico Poloni's user avatar
3 votes

Parallelize pseudo inverse of a matrix using Lapacke

LAPACK, and consequently LAPACKE are technically just interfaces. There are several implementations, some of which are already parallel. OpenBLAS and MKL should scale quite well on multicore machines....
Thijs Steel's user avatar
  • 1,693
3 votes

Automatic differentiation (AD) of a loss function which maps unitary matrix onto number

Your analytic derivative expression doesn't look right. Let's calculate the gradient (in the Wirtinger sense) $$\eqalign{ \def\l{L} \def\o{{\tt1}} \def\p{\partial} \def\grad#1#2{\frac{\p #1}{\p #2}} \...
greg's user avatar
  • 604
3 votes

How do BLAS libraries implement support for transposed matrices?

Blas kernels typically don't deal with transposed matrices. Instead, they will repack matrices (transposed or not) outside of the kernel. See https://www.cs.utexas.edu/users/flame/pubs/blis3_ipdps14....
Oscar Smith's user avatar
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 ...
aviks's user avatar
  • 121
2 votes
Accepted

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

I suggest you take a look at the Eigen C++ matrix class library, http://eigen.tuxfamily.org/index.php?title=Main_Page. I did a small test with sparse matrices of the size and sparsity you state and ...
Bill Greene's user avatar
  • 6,064
2 votes

Do BLAS routines compute their respective operations with minimum error?

Unfortunately, the "standard" BLAS implementations does not provide this feature. But due to the goals the scientific community want to achieve with the Exascale Computing serval "reproducible" BLAS ...
M.K. aka Grisu's user avatar
2 votes

Efficient change of basis real positive definite symmetric matrix

For the most part, LAPACK represents unitary matrices like $\mathbf U$ as products of householder reflectors, and provides specialized routines to work with that representation (for instance, use [...
rchilton1980's user avatar
  • 4,862

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