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 ...
- 726
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 ...
- 2,229
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
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 ...
- 11.1k
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 ...
- 186
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
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,229
6
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 ...
- 561
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 ...
- 1,773
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.
...
- 1,305
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 ...
- 3,994
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. ...
- 181
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} ...
- 1,386
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,602
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 ...
- 54.1k
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,602
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 ...
- 12.2k
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 ...
- 11.1k
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 ...
- 558
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 ...
- 18.5k
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 ...
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....
- 1,386
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}}
\...
- 594
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 ...
- 726
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
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
...
- 5,974
2
votes
Accepted
Using LAPACK to compute $B^{-1}AB^{-T}$ for thin $B$
You should be able to do this efficiently using LAPACK/BLAS, using the QR factorization [geqrf], orthogonal multiplication [ormqr], and triangular solve [trsm] routines. LAPACK and BLAS place a lot of ...
- 4,822
2
votes
BLAS operation question
Assuming we are looking at standard BLAS interface, I think you cannot really take advantage of the upper-triangular structure of matrix $U$, so you are back to a general matrix-matrix product via <...
- 8,602
2
votes
Beating typical BLAS libraries matrix multiplication performance
If we consider vectorization, are we going to beat typical libraries like OpenBLAS?
I do not fully agree with the other answer and would like to say kind of. It is true that libraries like Intel (R) ...
- 465
Only top scored, non community-wiki answers of a minimum length are eligible