Skip to content

Instantly share code, notes, and snippets.

@N0rbert
Last active February 17, 2024 14:55
Show Gist options
  • Star 2 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save N0rbert/cfda101b8f0aa326df1edb6beee0076d to your computer and use it in GitHub Desktop.
Save N0rbert/cfda101b8f0aa326df1edb6beee0076d to your computer and use it in GitHub Desktop.
This test script is a result of a discussion on AskUbuntu (https://askubuntu.com/q/1265756/66509 ) about using Intel MKL in Ubuntu
#!/bin/sh
# This test script is a result of a discussion on AskUbuntu (https://askubuntu.com/q/1265756/66509 ) about using Intel MKL in Ubuntu
#
# Modern Ubuntu versions include Intel MKL libraries since Ubuntu 19.10 (such as https://packages.ubuntu.com/focal/libmkl-full-dev )
# and it is expected that this library may be used by some scientific applicaitons like Octave, Scilab and others.
#
# Test method:
# 1. Execute this script with default mathematical libraries, save the results
# 2. Install the Intel MKL library with `sudo apt-get install libmkl-full-dev` confirming its usage as default math libraries alternative
# 3. Execute this script again after Intel MKL installation to compare the results
# Matrix size for all tests
size=500
# Known maximal eigenvalue for size 500
eigvalue_500=16.914886
# Octave test program is written by Archisman Panigrahi, adapted by Norbert
cat << EOF > /tmp/test_octave.m
% Octave test program is written by Archisman Panigrahi, adapted by Norbert
M = $size;
i = 1:M;
c = sin(i' + i.^2);
tic;
g = eig(c);
toc
m = max(real(g))
if (M == 500)
assert (m, $eigvalue_500, 1e-6)
%Correct result is ans = 16.915
%With MKL in Ubuntu 20.04, I get random numbers of order 10^5 - 10^6, which changes every time
endif
EOF
# Scilab
cat << EOF > /tmp/test_scilab.sce
mode(0);
M = $size;
i = 1:M;
c = sin( repmat(i', 1, length(i)) + repmat(i, length(i), 1).^2 );
tic();
g = spec(c);
t = toc()
m = max(real(g))
if M==500 then
assert_checkalmostequal(m, $eigvalue_500, 1e-6);
end
EOF
# Julia
cat << EOF > /tmp/test_julia.jl
using LinearAlgebra
M = $size;
c = [ sin(a + b.^2) for a in 1:M, b in 1:M]
t0 = time();
g = LinearAlgebra.eigvals(c);
t = time() - t0
println("t=", t)
m = maximum(real(g))
println("m=", m)
if M == 500
@assert isapprox(m, $eigvalue_500, rtol=1e-6)
end
EOF
# Python 3 with Numpy
cat << EOF > /tmp/test_numpy.py
import numpy
from time import time
M = $size;
c = numpy.empty((M, M));
for i in range(1,M+1):
for j in range(1,M+1):
c[i-1][j-1] = numpy.sin(i + numpy.power(j, 2.))
t0 = time();
g = numpy.linalg.eigvals(c);
t = time() - t0;
print("t = ", t)
m = numpy.max(numpy.real(g));
print("m = ", m)
if M == 500:
numpy.testing.assert_almost_equal(m, $eigvalue_500, 6)
EOF
# Python 3 with Torch and CUDA
cat << EOF > /tmp/test_torch.py
import torch
import numpy
from time import time
M = 500;
c = numpy.empty((M, M));
for i in range(1,M+1):
for j in range(1,M+1):
c[i-1][j-1] = numpy.sin(i + numpy.power(j, 2.))
t0 = time();
c_gpu = torch.tensor(c)
g = torch.linalg.eigvals(c_gpu)
t = time() - t0;
print("t = ", t)
m = torch.max(torch.real(g));
print("m = ", m)
if M == 500:
numpy.testing.assert_almost_equal(m, $eigvalue_500, 6)
EOF
# R
cat << EOF > /tmp/test_r.R
# needs r-cran-assertthat deb-package
M = $size;
c = matrix(nrow=M, ncol=M)
for (i in 1:M) {
for (j in 1:M) {
c[i, j] = sin(i + j**2);
}
}
t0 = Sys.time();
g = eigen(c);
t = Sys.time() - t0;
cat(sprintf("t = %f\n", t));
m = max(Re(g[["values"]]));
cat(sprintf("m = %f\n", m));
if (M == 500) {
assertthat::validate_that(abs(m - $eigvalue_500) < 1e-6);
}
EOF
###### Call test programs
echo
echo "> Scilab"
scilab-cli -nb -noatomsautoload -f /tmp/test_scilab.sce -quit
echo
echo "> Julia"
julia /tmp/test_julia.jl
echo
echo "> Python 3 with NumPy"
python3 /tmp/test_numpy.py
echo
echo "> Python 3 with Torch and CUDA"
python3 /tmp/test_torch.py
echo
echo "> R (will produce error for MKL)"
Rscript /tmp/test_r.R
echo
echo "> Octave (will produce error for MKL)"
octave-cli /tmp/test_octave.m
echo
echo "> Octave (with MKL error fixed)"
# need to define MKL_THREADING_LAYER environment variable
dpkg -l | grep --silent mkl && export MKL_THREADING_LAYER=gnu
octave-cli /tmp/test_octave.m
echo
echo "> R (with MKL error fixed)"
Rscript /tmp/test_r.R
@archisman-panigrahi
Copy link

archisman-panigrahi commented Aug 11, 2023

Turns out that Julia has MKL. One just need to do using MKL after installing the MKL package. I will run the script in my workstation once again

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment