diff --git a/cscs-checks/apps/python/numpy_check.py b/cscs-checks/apps/python/numpy_check.py
new file mode 100644
index 0000000000..2185c465ca
--- /dev/null
+++ b/cscs-checks/apps/python/numpy_check.py
@@ -0,0 +1,87 @@
+import reframe as rfm
+import reframe.utility.sanity as sn
+
+
+class NumpyBaseTest(rfm.RunOnlyRegressionTest):
+    def __init__(self):
+        self.descr = 'Test a few typical numpy operations'
+        self.valid_prog_environs = ['PrgEnv-gnu']
+        self.modules = ['numpy/1.17.2-CrayGNU-19.10']
+        self.reference = {
+            'daint:gpu': {
+                'dot': (0.4, None, 0.05, 'seconds'),
+                'svd': (0.37, None, 0.05, 'seconds'),
+                'cholesky': (0.12, None, 0.05, 'seconds'),
+                'eigendec': (3.5, None, 0.05, 'seconds'),
+                'inv': (0.21, None, 0.05, 'seconds'),
+            },
+            'daint:mc': {
+                'dot': (0.3, None, 0.05, 'seconds'),
+                'svd': (0.35, None, 0.05, 'seconds'),
+                'cholesky': (0.1, None, 0.05, 'seconds'),
+                'eigendec': (4.14, None, 0.05, 'seconds'),
+                'inv': (0.16, None, 0.05, 'seconds'),
+            },
+            'dom:gpu': {
+                'dot': (0.4, None, 0.05, 'seconds'),
+                'svd': (0.37, None, 0.05, 'seconds'),
+                'cholesky': (0.12, None, 0.05, 'seconds'),
+                'eigendec': (3.5, None, 0.05, 'seconds'),
+                'inv': (0.21, None, 0.05, 'seconds'),
+            },
+            'dom:mc': {
+                'dot': (0.3, None, 0.05, 'seconds'),
+                'svd': (0.35, None, 0.05, 'seconds'),
+                'cholesky': (0.1, None, 0.05, 'seconds'),
+                'eigendec': (4.14, None, 0.05, 'seconds'),
+                'inv': (0.16, None, 0.05, 'seconds'),
+            },
+        }
+        self.perf_patterns = {
+            'dot': sn.extractsingle(
+                r'^Dotted two 4096x4096 matrices in\s+(?P<dot>\S+)\s+s',
+                self.stdout, 'dot', float),
+            'svd': sn.extractsingle(
+                r'^SVD of a 2048x1024 matrix in\s+(?P<svd>\S+)\s+s',
+                self.stdout, 'svd', float),
+            'cholesky': sn.extractsingle(
+                r'^Cholesky decomposition of a 2048x2048 matrix in'
+                r'\s+(?P<cholesky>\S+)\s+s',
+                self.stdout, 'cholesky', float),
+            'eigendec': sn.extractsingle(
+                r'^Eigendecomposition of a 2048x2048 matrix in'
+                r'\s+(?P<eigendec>\S+)\s+s',
+                self.stdout, 'eigendec', float),
+            'inv': sn.extractsingle(
+                r'^Inversion of a 2048x2048 matrix in\s+(?P<inv>\S+)\s+s',
+                self.stdout, 'inv', float)
+        }
+        self.sanity_patterns = sn.assert_found(r'Numpy version:\s+\S+',
+                                               self.stdout)
+        self.variables = {
+            'OMP_NUM_THREADS': '$SLURM_CPUS_PER_TASK',
+        }
+        self.executable = 'python'
+        self.executable_opts = ['np_ops.py']
+        self.num_tasks_per_node = 1
+        self.use_multithreading = False
+        self.tags = {'production'}
+        self.maintainers = ['RS', 'TR']
+
+
+@rfm.required_version('>=2.16')
+@rfm.simple_test
+class NumpyHaswellTest(NumpyBaseTest):
+    def __init__(self):
+        super().__init__()
+        self.valid_systems = ['daint:gpu', 'dom:gpu']
+        self.num_cpus_per_task = 12
+
+
+@rfm.required_version('>=2.16')
+@rfm.simple_test
+class NumpyBroadwellTest(NumpyBaseTest):
+    def __init__(self):
+        super().__init__()
+        self.valid_systems = ['daint:mc', 'dom:mc']
+        self.num_cpus_per_task = 36
diff --git a/cscs-checks/apps/python/src/np_ops.py b/cscs-checks/apps/python/src/np_ops.py
new file mode 100644
index 0000000000..dd046dcbf0
--- /dev/null
+++ b/cscs-checks/apps/python/src/np_ops.py
@@ -0,0 +1,68 @@
+import numpy as np
+from time import time
+
+
+print('Numpy version:', np.__version__)
+
+# Let's take the randomness out of random numbers (for reproducibility)
+np.random.seed(0)
+
+size = 4096
+A, B = np.random.random((size, size)), np.random.random((size, size))
+C, D = np.random.random((size * 128,)), np.random.random((size * 128,))
+E = np.random.random((int(size / 2), int(size / 4)))
+F = np.random.random((int(size / 2), int(size / 2)))
+F = np.dot(F, F.T)
+G = np.random.random((int(size / 2), int(size / 2)))
+I = (F + F.T) / 2.
+
+# Matrix multiplication
+N = 20
+t = time()
+for i in range(N):
+    np.dot(A, B)
+
+delta = time() - t
+print('Dotted two %dx%d matrices in %0.2f s.' % (size, size, delta / N))
+del A, B
+
+# Singular Value Decomposition (SVD)
+N = 3
+t = time()
+for i in range(N):
+    np.linalg.svd(E, full_matrices=False)
+
+delta = time() - t
+print("SVD of a %dx%d matrix in %0.2f s." % (size / 2, size / 4, delta / N))
+del E
+
+# Cholesky Decomposition
+N = 3
+t = time()
+for i in range(N):
+    np.linalg.cholesky(F)
+
+delta = time() - t
+print("Cholesky decomposition of a %dx%d matrix in %0.2f s." % (size / 2,
+                                                                size / 2,
+                                                                delta / N))
+
+# Eigendecomposition
+t = time()
+for i in range(N):
+    np.linalg.eig(G)
+
+delta = time() - t
+print("Eigendecomposition of a %dx%d matrix in %0.2f s." % (size / 2,
+                                                            size / 2,
+                                                            delta / N))
+
+# Matrix Inversion
+t = time()
+for i in range(N):
+    np.linalg.inv(I)
+
+delta = time() - t
+print("Inversion of a %dx%d matrix in %0.2f s." % (size / 2,
+                                                   size / 2,
+                                                   delta / N))