Repository: alan-turing-institute/mogp-emulator (language: python)

Issue:
Feature/gpu meanfunc
`feature/gpu-meanfunc` branch was created to add minimal mean-function implementation on GaussianProcessGPU (and associated C++/Cuda), while leaving the `feature/gpu` branch with minimal set of features that would be ready to merge with main.
However, this branch currently has several fixes that should really be on the `feature/gpu` branch to be merged:
* Fix to sumLogDiag for N>32
* Allow fit_GP_MAP to accept GaussianProcessGPU objects
* Allow MultiOutputGPU to accept GaussianProcessGPU objects (though currently doesn't work, as we can't yet pickle the GaussianProcessGPU).
Fix threadblock size calculation in some of the CUDA kernel calls

Rename variables in mogp_gpu to be consistent with the python implementation



Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index 8746f4ac..c06e4496 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,6 +1,9 @@
 *.pyc
 *.DS_Store
+*~
 docs/_build/*
+mogp_gpu/obj
+mogp_gpu/bin
 mogp_emulator/benchmarks/*.png
 mogp_emulator/benchmarks/*.pdf
 mogp_emulator/benchmarks/*.npy
diff --git a/docs/demos/gp_demoGPU.rst b/docs/demos/gp_demoGPU.rst
new file mode 100644
index 00000000..a7111dd7
--- /dev/null
+++ b/docs/demos/gp_demoGPU.rst
@@ -0,0 +1,16 @@
+.. _gpdemoGPU:
+
+Gaussian Process Demo (GPU)
+=========================
+This demo illustrates a simple example of fitting a GP emulator to results of the projectile
+problem discussed in the :ref:`tutorial`, using the GPU implementation of the emulator.
+
+Note that in order for this to work, it must be run on a machine with an Nvidia GPU, and with
+CUDA libraries available.  It also depends on Eigen and pybind.
+
+The example uses Maximum Likelihood Estimation with a Squared Exponential kernel, which is currently
+the only kernel supported by the GPU implementation.
+
+
+.. literalinclude::
+   ../../mogp_emulator/demos/gp_demo_gpu.py
diff --git a/docs/implementation/GaussianProcessGPU.rst b/docs/implementation/GaussianProcessGPU.rst
new file mode 100644
index 00000000..09304361
--- /dev/null
+++ b/docs/implementation/GaussianProcessGPU.rst
@@ -0,0 +1,15 @@
+.. _GaussianProcessGPU:
+
+**********************************
+The ``GaussianProcessGPU`` Class
+**********************************
+
+.. automodule:: mogp_emulator.GaussianProcessGPU
+    :noindex:
+
+.. autoclass:: mogp_emulator.GaussianProcessGPU.GaussianProcessGPU
+    :members:
+    
+    .. automethod:: __init__
+
+.. automodule:: libgpgpu
diff --git a/docs/implementation/implementation.rst b/docs/implementation/implementation.rst
index 2ac42695..df4186fe 100644
--- a/docs/implementation/implementation.rst
+++ b/docs/implementation/implementation.rst
@@ -8,6 +8,7 @@ mogp_emulator Implementation Details
    :caption: Class and Function Implementation Details:
 
    GaussianProcess
+   GaussianProcessGPU
    MultiOutputGP
    fitting
    MeanFunction
diff --git a/docs/index.rst b/docs/index.rst
index 54b6d03f..31330c02 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -35,6 +35,7 @@ details, and some included benchmarks.
    demos/mice_demos
    demos/historymatch_demos
    demos/kdr_demos
+   demos/gp_demoGPU
    demos/gp_demoR
    demos/excalibur_workshop_demo
 
diff --git a/mogp_emulator/GaussianProcess.py b/mogp_emulator/GaussianProcess.py
index f05cc0cd..bcca8ea8 100644
--- a/mogp_emulator/GaussianProcess.py
+++ b/mogp_emulator/GaussianProcess.py
@@ -6,7 +6,12 @@
 from scipy.optimize import OptimizeResult
 from mogp_emulator.linalg import jit_cholesky, pivot_cholesky, pivot_cho_solve
 
-class GaussianProcess(object):
+
+class GaussianProcessBase(object):
+    pass
+
+
+class GaussianProcess(GaussianProcessBase):
     """Implementation of a Gaussian Process Emulator.
 
     This class provides a representation of a Gaussian Process
diff --git a/mogp_emulator/GaussianProcessGPU.py b/mogp_emulator/GaussianProcessGPU.py
new file mode 100644
index 00000000..91a941f3
--- /dev/null
+++ b/mogp_emulator/GaussianProcessGPU.py
@@ -0,0 +1,416 @@
+"""
+extends GaussianProcess with an (optional) GPU implementation
+"""
+
+import os
+import numpy as np
+from mogp_emulator.Kernel import SquaredExponential
+
+import mogp_emulator.LibGPGPU as LibGPGPU
+
+from mogp_emulator.GaussianProcess import GaussianProcessBase, PredictResult
+
+
+class GPUUnavailableError(RuntimeError):
+    """Exception type to use when a GPU, or the GPU library, is unavailable"""
+    pass
+
+
+def ndarray_coerce_type_and_flags(arr):
+    """
+    Helper function for the GaussianProcessGPU methods that call
+    CUDA/C++ functions (those wrapped by _dense_gpgpu) and that take
+    numpy arrays as arguments.  Ensures that an array is of the
+    correct type for this purpose.
+
+    Takes an array or array-like.
+
+    Returns an ndarray with the same data, that has:
+    - dtype of float64
+    - flags.writable
+    - flags.c_contiguous
+
+    The array returned may reference the original array, be a copy of
+    it, or be newly constructed from the argument.
+    """
+
+    arr = np.array(arr, copy=False)
+
+    # cast into float64, just in case we were given integers and ensure contiguous (C type)
+    arr_contiguous_float64 = np.ascontiguousarray(arr.astype(np.float64, copy=False))
+    if not arr_contiguous_float64.flags['WRITEABLE']:
+        return np.copy(arr_contiguous_float64)
+    else:
+        return arr_contiguous_float64
+
+
+class GaussianProcessGPU(GaussianProcessBase):
+
+    """
+    This class implements the same interface as
+    :class:`mogp_emulator.GaussianProcess.GaussianProcess`, but using a GPU if available.
+    Will raise a RuntimeError if a CUDA-compatible GPU, GPU-interface library libgpgpu
+    could not be found.
+    """
+
+    def __init__(self, inputs, targets, mean=None, kernel=SquaredExponential(), priors=None,
+                 nugget="adaptive", inputdict = {}, use_patsy=True, max_batch_size=2000):
+
+        if not LibGPGPU.HAVE_LIBGPGPU:
+            raise RuntimeError("Cannot construct GaussianProcessGPU: "
+                               "The GPU library (libgpgpu) could not be loaded")
+
+        elif not LibGPGPU.gpu_usable():
+            raise RuntimeError("Cannot construct GaussianProcessGPU: "
+                               "A compatible GPU could not be found")
+
+        inputs = ndarray_coerce_type_and_flags(inputs)
+        if inputs.ndim == 1:
+            inputs = np.reshape(inputs, (-1, 1))
+        assert inputs.ndim == 2
+
+        targets = ndarray_coerce_type_and_flags(targets)
+        assert targets.ndim == 1
+        assert targets.shape[0] == inputs.shape[0]
+
+        self._inputs = inputs
+        self._targets = targets
+        self._max_batch_size = max_batch_size
+
+        if mean:
+            raise ValueError("GPU implementation requires mean to be None")
+
+        if isinstance(kernel, str):
+            if kernel == "SquaredExponential":
+                kernel = SquaredExponential()
+            else:
+                raise ValueError("GPU implementation requires kernel to be SquaredExponential")
+        elif kernel and not isinstance(kernel, SquaredExponential):
+                raise ValueError("GPU implementation requires kernel to be SquaredExponential()")
+        self.kernel = kernel
+        self.nugget = nugget
+
+        self._densegp_gpu = None
+        self._init_gpu()
+
+
+    def _init_gpu(self):
+        """
+        Instantiate the DenseGP_GPU C++/CUDA class, if it doesn't already exist.
+        """
+        if not self._densegp_gpu:
+            self._densegp_gpu = LibGPGPU.DenseGP_GPU(self._inputs, self._targets, self._max_batch_size)
+
+
+    @property
+    def inputs(self):
+        """
+        Returns inputs for the emulator as a numpy array
+
+        :returns: Emulator inputs, 2D array with shape ``(n, D)``
+        :rtype: ndarray
+        """
+        return self._densegp_gpu.inputs()
+
+    @property
+    def targets(self):
+        """
+        Returns targets for the emulator as a numpy array
+
+        :returns: Emulator targets, 1D array with shape ``(n,)``
+        :rtype: ndarray
+        """
+        return self._densegp_gpu.targets()
+
+    @property
+    def n(self):
+        """
+        Returns number of training examples for the emulator
+
+        :returns: Number of training examples for the emulator object
+        :rtype: int
+        """
+        return self._densegp_gpu.n()
+
+    @property
+    def D(self):
+        """
+        Returns number of inputs (dimensions) for the emulator
+
+        :returns: Number of inputs for the emulator object
+        :rtype: int
+        """
+        return self._densegp_gpu.D()
+
+    @property
+    def n_params(self):
+        """
+        Returns number of hyperparameters
+
+        Returns the number of hyperparameters for the emulator. The number depends on the
+        choice of mean function, covariance function, and nugget strategy, and possibly the
+        number of inputs for certain choices of the mean function.
+
+        :returns: Number of hyperparameters
+        :rtype: int
+        """
+        return self._densegp_gpu.n_params()
+
+    @property
+    def nugget_type(self):
+        """
+        Returns method used to select nugget parameter
+
+        Returns a string indicating how the nugget parameter is treated, either ``"adaptive"``,
+        ``"fit"``, or ``"fixed"``. This is automatically set when changing the ``nugget``
+        property.
+
+        :returns: Current nugget fitting method
+        :rtype: str
+        """
+        return self._nugget_type.__str__().split(".")[1]
+
+    @property
+    def nugget(self):
+        """
+        See :func:`mogp_emulator.GaussianProcess.GaussianProcess.nugget`
+        """
+        return self._nugget
+
+    @nugget.setter
+    def nugget(self, nugget):
+        if not isinstance(nugget, (str, float)):
+            try:
+                nugget = float(nugget)
+            except TypeError:
+                raise TypeError("nugget parameter must be a string or a non-negative float")
+
+        if isinstance(nugget, str):
+            if nugget == "adaptive":
+                self._nugget_type = LibGPGPU.nugget_type(0)
+            elif nugget == "fit":
+                self._nugget_type = LibGPGPU.nugget_type(1)
+            else:
+                raise ValueError("nugget must be a float set to 'adaptive', 'fit', or 'fixed'")
+            self._nugget = None
+        else:
+            if nugget < 0.:
+                raise ValueError("nugget parameter must be non-negative")
+            self._nugget_type = LibGPGPU.nugget_type(2) #fixed
+            self._nugget = float(nugget)
+
+    @property
+    def theta(self):
+        """
+        Returns emulator hyperparameters
+        see
+        :func:`mogp_emulator.GaussianProcess.GaussianProcess.theta`
+
+        :type theta: ndarray
+        """
+        if not self._densegp_gpu.theta_fit_status():
+            return None
+        theta = np.zeros(self.n_params)
+        self._densegp_gpu.get_theta(theta)
+        return theta
+
+    @theta.setter
+    def theta(self, theta):
+        """
+        Fits the emulator and sets the parameters (property-based setter
+        alias for ``fit``)
+
+        See :func:`mogp_emulator.GaussianProcess.GaussianProcess.theta`
+
+        :type theta: ndarray
+        :returns: None
+        """
+        if theta is None:
+            self._densegp_gpu.reset_theta_fit_status()
+            self.current_logpost = None
+        else:
+            self.fit(theta)
+
+    @property
+    def L(self):
+        """
+        Return the lower triangular Cholesky factor.
+
+        :returns: np.array
+        """
+        result = np.zeros((self.n, self.n))
+        self._densegp_gpu.get_cholesky_lower(result)
+        return np.tril(result.transpose())
+
+    def get_K_matrix(self):
+        """
+        Returns current value of the inverse covariance matrix as a numpy array.
+
+        Does not include the nugget parameter, as this is dependent on how the
+        nugget is fit.
+        """
+        result = np.zeros((self.n, self.n))
+        self._densegp_gpu.get_K(result)
+        return result
+
+    def fit(self, theta):
+        """
+        Fits the emulator and sets the parameters.
+
+        Implements the same interface as
+        :func:`mogp_emulator.GaussianProcess.GaussianProcess.fit`
+        """
+        theta = ndarray_coerce_type_and_flags(theta)
+
+        nugget_size = self.nugget if self.nugget_type == "fixed" else 0.
+        self._densegp_gpu.fit(theta, self._nugget_type, nugget_size)
+        if self.nugget_type == "adaptive" or self.nugget_type == "fit":
+            self._nugget = self._densegp_gpu.get_jitter()
+        invQt_result = np.zeros(self.n)
+        self._densegp_gpu.get_invQt(invQt_result)
+        self.invQt = invQt_result
+        self.current_logpost = self.logposterior(theta)
+
+    def logposterior(self, theta):
+        """
+        Calculate the negative log-posterior at a particular value of the hyperparameters
+
+        See :func:`mogp_emulator.GaussianProcess.GaussianProcess.logposterior`
+
+        :param theta: Value of the hyperparameters. Must be array-like with shape ``(n_params,)``
+        :type theta: ndarray
+        :returns: negative log-posterior
+        :rtype: float
+        """
+        if self.theta is None or not np.allclose(theta, self.theta, rtol=1.e-10, atol=1.e-15):
+            self.fit(theta)
+
+        return self._densegp_gpu.get_logpost()
+
+    def logpost_deriv(self, theta):
+        """
+        Calculate the partial derivatives of the negative log-posterior
+
+        See :func:`mogp_emulator.GaussianProcess.GaussianProcess.logpost_deriv`
+
+        :param theta: Value of the hyperparameters. Must be array-like with shape ``(n_params,)``
+        :type theta: ndarray
+        :returns: partial derivatives of the negative log-posterior with respect to the
+                  hyperparameters (array with shape ``(n_params,)``)
+        :rtype: ndarray
+        """
+        theta = np.array(theta, copy=False)
+
+        assert theta.shape == (self.n_params,), "bad shape for new parameters"
+
+        if self.theta is None or not np.allclose(theta, self.theta, rtol=1.e-10, atol=1.e-15):
+            self.fit(theta)
+
+        result = np.zeros(self.n_params)
+        self._densegp_gpu.logpost_deriv(result)
+        return result
+
+    def logpost_hessian(self, theta):
+        """
+        Calculate the Hessian of the negative log-posterior
+
+        See :func:`mogp_emulator.GaussianProcess.GaussianProcess.logpost_hessian`
+
+        :param theta: Value of the hyperparameters. Must be array-like with shape
+                      ``(n_params,)``
+        :type theta: ndarray
+        :returns: Hessian of the negative log-posterior (array with shape
+                  ``(n_params, n_params)``)
+        :rtype: ndarray
+        """
+        raise GPUUnavailableError(
+            "The Hessian calculation is not currently implemented in the GPU version of MOGP."
+        )
+
+
+    def predict(self, testing, unc=True, deriv=True, include_nugget=True):
+        """
+        Make a prediction for a set of input vectors for a single set of hyperparameters.
+        This method implements the same interface as
+        :func:`mogp_emulator.GaussianProcess.GaussianProcess.predict`
+        """
+
+        if self.theta is None:
+            raise ValueError("hyperparameters have not been fit for this Gaussian Process")
+
+        testing = ndarray_coerce_type_and_flags(testing)
+
+        if self.D == 1 and testing.ndim == 1:
+            testing = np.reshape(testing, (-1, 1))
+        elif testing.ndim == 1:
+            testing = np.reshape(testing, (1, len(testing)))
+        assert testing.ndim == 2
+
+        n_testing, D = np.shape(testing)
+
+        assert D == self.D
+
+        means = np.zeros(n_testing)
+
+        if unc:
+            variances = np.zeros(n_testing)
+            for i in range(0, n_testing, self._max_batch_size):
+                self._densegp_gpu.predict_variance_batch(
+                    testing[i:i+self._max_batch_size],
+                    means[i:i+self._max_batch_size],
+                    variances[i:i+self._max_batch_size])
+            if include_nugget:
+                variances += self._nugget
+        else:
+            for i in range(0, n_testing, self._max_batch_size):
+                self._densegp_gpu.predict_batch(
+                    testing[i:i+self._max_batch_size],
+                    means[i:i+self._max_batch_size])
+            variances = None
+        if deriv:
+            deriv_result = np.zeros((n_testing,self.D))
+            for i in range(0, n_testing, self._max_batch_size):
+                self._densegp_gpu.predict_deriv(
+                    testing[i:i+self._max_batch_size],
+                    deriv_result[i:i+self._max_batch_size])
+        else:
+            deriv_result = None
+        return PredictResult(mean=means, unc=variances, deriv=deriv_result)
+
+
+    def __call__(self, testing):
+        """A Gaussian process object is callable: calling it is the same as
+        calling `predict` without uncertainty and derivative
+        predictions, and extracting the zeroth component for the
+        'mean' prediction.
+        """
+        return (self.predict(testing, unc=False, deriv=False)[0])
+
+
+    def __str__(self):
+        """
+        Returns a string representation of the model
+
+        :returns: A string representation of the model
+        (indicates number of training examples and inputs)
+        :rtype: str
+        """
+        return ("Gaussian Process with " + str(self.n) + " training examples and " +
+                str(self.D) + " input variables")
+
+
+    ## __setstate__ and __getstate__ for pickling: don't pickle "_dense_gp_gpu",
+    ## and instead reinitialize this when unpickling
+    ## (Pickling is required to use multiprocessing.)
+    def __setstate__(self, state):
+        self.__dict__ = state
+        self.init_gpu()
+        if self._theta is not None:
+            self.fit(self._theta)
+
+
+    def __getstate__(self):
+        copy_dict = self.__dict__.copy()
+        del copy_dict["_densegp_gpu"]
+        copy_dict["_densegp_gpu"] = None
+        return copy_dict
diff --git a/mogp_emulator/HistoryMatching.py b/mogp_emulator/HistoryMatching.py
index e128f132..a3149948 100644
--- a/mogp_emulator/HistoryMatching.py
+++ b/mogp_emulator/HistoryMatching.py
@@ -1,5 +1,5 @@
 import numpy as np
-from mogp_emulator.GaussianProcess import GaussianProcess, PredictResult
+from mogp_emulator.GaussianProcess import GaussianProcessBase, PredictResult
 
 class HistoryMatching(object):
     r"""
@@ -80,10 +80,10 @@ def __init__(self, gp=None, obs=None, coords=None, expectations=None, threshold=
         numbers of these to allow the number of dimensions and/or number of coordinates to
         be computed also causes these to be set at this stage.
 
-        :param gp: (Optional) ``GaussianProcess`` object used to make predictions in
+        :param gp: (Optional) ``GaussianProcessBase`` object used to make predictions in
                    history matching. Optional, can instead provide predictions directly
                    via the ``expectations`` argument. Default is ``None``
-        :type gp: GaussianProcess or None
+        :type gp: GaussianProcessBase or None
         :param obs: (Optional) Observations against which the predictions will be
                     compared. If provided, must either be a float (assumes no
                     uncertainty in the observations) or a list of two floats
@@ -438,17 +438,15 @@ def check_gp(self, gp):
         Checks if the provided argument is consistent with expectations for a GP.
 
         Returns a boolean that is True if the provided quantity is consistent with
-        the requirements for a gaussian process, i.e. is of type ``GaussianProcess``.
+        the requirements for a gaussian process, i.e. is of type ``GaussianProcessBase``.
 
         :param gp: Input GP object to be checked.
-        :type gp: GaussianProcess
+        :type gp: GaussianProcessBase
         :returns: Boolean indicating if provided object is a ``GaussianProcess``
         :rtype: bool
         """
-        if gp is None: return False
-        if isinstance(gp, GaussianProcess):
-            return True
-        return False
+        return isinstance(gp, GaussianProcessBase)
+
 
 
     def check_obs(self, obs):
diff --git a/mogp_emulator/LibGPGPU.py b/mogp_emulator/LibGPGPU.py
new file mode 100644
index 00000000..939629a3
--- /dev/null
+++ b/mogp_emulator/LibGPGPU.py
@@ -0,0 +1,14 @@
+"""
+Wrapper for the optional GPU interface library, libgpgpu.so
+"""
+
+HAVE_LIBGPGPU = False
+try:
+    from libgpgpu import *
+    HAVE_LIBGPGPU = True
+except ModuleNotFoundError:
+    pass
+
+
+def gpu_usable():
+    return HAVE_LIBGPGPU and have_compatible_device()
diff --git a/mogp_emulator/MultiOutputGP.py b/mogp_emulator/MultiOutputGP.py
index 8fea2927..23458f9d 100644
--- a/mogp_emulator/MultiOutputGP.py
+++ b/mogp_emulator/MultiOutputGP.py
@@ -1,11 +1,18 @@
 from multiprocessing import Pool
 import platform
 import numpy as np
-from mogp_emulator.GaussianProcess import GaussianProcess, PredictResult
+from mogp_emulator.GaussianProcess import (
+    GaussianProcessBase,
+    GaussianProcess,
+    PredictResult
+)
+from mogp_emulator.GaussianProcessGPU import GaussianProcessGPU
 from mogp_emulator.MeanFunction import MeanBase
 from mogp_emulator.Kernel import Kernel, SquaredExponential, Matern52
 from mogp_emulator.Priors import Prior
 
+
+
 class MultiOutputGP(object):
     """Implementation of a multiple-output Gaussian Process Emulator.
 
@@ -32,11 +39,18 @@ class MultiOutputGP(object):
     """
 
     def __init__(self, inputs, targets, mean=None, kernel=SquaredExponential(), priors=None,
-                 nugget="adaptive", inputdict={}, use_patsy=True):
+                 nugget="adaptive", inputdict={}, use_patsy=True, use_gpu=False):
         """
         Create a new multi-output GP Emulator
         """
 
+        # if use_gpu is selected, check whether we found the GPU .so file, and raise error if not
+        self.use_gpu = use_gpu
+        if self.use_gpu:
+            self.GPClass = GaussianProcessGPU
+        else:
+            self.GPClass = GaussianProcess
+
         # check input types and shapes, reshape as appropriate for the case of a single emulator
         inputs = np.array(inputs)
         targets = np.array(targets)
@@ -92,7 +106,7 @@ def __init__(self, inputs, targets, mean=None, kernel=SquaredExponential(), prio
         assert isinstance(nugget, list), "nugget must be a string, float, or a list of strings and floats"
         assert len(nugget) == self.n_emulators
 
-        self.emulators = [ GaussianProcess(inputs, single_target, m, k, p, n, inputdict, use_patsy)
+        self.emulators = [ self.GPClass(inputs, single_target, m, k, p, n, inputdict, use_patsy)
                            for (single_target, m, k, p, n) in zip(targets, mean, kernel, priors, nugget)]
 
 
@@ -196,9 +210,9 @@ def predict(self, testing, unc=True, deriv=True, include_nugget=True,
         if allow_not_fit:
             predict_method = _gp_predict_default_NaN
         else:
-            predict_method = GaussianProcess.predict
-            
-        if platform.system() == "Windows":
+            predict_method = self.GPClass.predict
+
+        if platform.system() == "Windows" or self.use_gpu:
             predict_vals = [predict_method(gp, testing, unc, deriv, include_nugget)
                             for gp in self.emulators]
         else:
@@ -228,7 +242,7 @@ def __call__(self, testing, process=None):
 
         return self.predict(testing, unc=False, deriv=False, processes=processes)[0]
 
-    
+
     def get_indices_fit(self):
         """Returns the indices of the emulators that have been fit
 
@@ -253,7 +267,7 @@ def get_indices_fit(self):
         return [idx for (failed_fit, idx)
                 in zip([em.theta is None for em in self.emulators],
                        list(range(len(self.emulators)))) if not failed_fit]
-    
+
     def get_indices_not_fit(self):
         """Returns the indices of the emulators that have not been fit
 
@@ -328,7 +342,7 @@ def get_emulators_not_fit(self):
                 in zip([em.theta is None for em in self.emulators],
                        self.emulators) if failed_fit]
 
-        
+
     def __str__(self):
         """Returns a string representation of the model
 
@@ -387,15 +401,15 @@ def _gp_predict_default_NaN(gp, testing, unc, deriv, include_nugget):
               arrays are replaced by ``None``.
     :rtype: tuple
     """
-    
-    assert isinstance(gp, GaussianProcess)
+
+    assert isinstance(gp, GaussianProcessBase)
 
     try:
         return gp.predict(testing, unc, deriv, include_nugget)
     except ValueError:
-        
+
         n_predict = testing.shape[0]
-        
+
         if unc:
             unc_array = np.array([np.nan]*n_predict)
         else:
diff --git a/mogp_emulator/__init__.py b/mogp_emulator/__init__.py
index 370cffc5..b46c6790 100644
--- a/mogp_emulator/__init__.py
+++ b/mogp_emulator/__init__.py
@@ -1,6 +1,11 @@
 from .version import version as __version__
 
 from .GaussianProcess import GaussianProcess
+try:
+    from .GaussianProcessGPU import GaussianProcessGPU
+    from .GaussianProcessGPU import GPUUnavailableError
+except:
+    pass
 from .MultiOutputGP import MultiOutputGP
 from .fitting import fit_GP_MAP
 from .MeanFunction import MeanFunction
diff --git a/mogp_emulator/demos/gp_demo_gpu.py b/mogp_emulator/demos/gp_demo_gpu.py
new file mode 100644
index 00000000..edc6c773
--- /dev/null
+++ b/mogp_emulator/demos/gp_demo_gpu.py
@@ -0,0 +1,65 @@
+import numpy as np
+import mogp_emulator
+from projectile import simulator, print_results
+
+# GP example using the projectile demo on a GPU
+
+# To run this demo you must be on a machine with an Nvidia GPU, and with
+# CUDA libraries available.  There are also dependencies on eigen and pybind11
+# If you are working on a managed cluster, these may be available via commands
+#
+# module load cuda/11.2
+# module load py-pybind11-2.2.4-gcc-5.4.0-tdtz6iq
+# module load gcc/7
+# module load eigen
+#
+# You should then be able to compile the cuda code at the same time as installing the mogp_emulator package, by doing (from the main mogp_emulator/ directory:
+# pip install .
+# (note that if you don't have write access to the global directory
+# (e.g. if you are on a cluster such as CSD3), you should add the
+# `--user` flag to this command)
+
+
+# define some common variables
+
+n_samples = 20
+n_preds = 10
+
+# Experimental design -- requires a list of parameter bounds if you would like to use
+# uniform distributions. If you want to use different distributions, you
+# can use any of the standard distributions available in scipy to create
+# the appropriate ppf function (the inverse of the cumulative distribution).
+# Internally, the code creates the design on the unit hypercube and then uses
+# the distribution to map from [0,1] to the real parameter space.
+
+ed = mogp_emulator.LatinHypercubeDesign([(-5., 1.), (0., 1000.)])
+
+# sample space
+
+inputs = ed.sample(n_samples)
+
+# run simulation
+
+targets = np.array([simulator(p) for p in inputs])
+
+###################################################################################
+
+# Basic example -- fit GP using MLE and Squared Exponential Kernel and predict
+
+print("Example: Basic GP")
+
+# create GP and then fit using MLE
+# the only difference between this and the standard CPU implementation
+# is to use the GaussianProcessGPU class rather than GaussianProcess.
+
+gp = mogp_emulator.GaussianProcessGPU(inputs, targets)
+
+gp = mogp_emulator.fit_GP_MAP(gp)
+
+# create 20 target points to predict
+
+predict_points = ed.sample(n_preds)
+
+means, variances, derivs = gp.predict(predict_points)
+
+print_results(predict_points, means)
diff --git a/mogp_emulator/fitting.py b/mogp_emulator/fitting.py
index f6675bc6..800ed180 100644
--- a/mogp_emulator/fitting.py
+++ b/mogp_emulator/fitting.py
@@ -2,12 +2,14 @@
 import scipy.stats
 from scipy.linalg import LinAlgError
 from scipy.optimize import minimize
-from mogp_emulator.GaussianProcess import GaussianProcess
-from mogp_emulator.MultiOutputGP import MultiOutputGP
 from multiprocessing import Pool
 from functools import partial
 import platform
 
+from mogp_emulator.GaussianProcess import GaussianProcessBase, GaussianProcess
+from mogp_emulator.MultiOutputGP import MultiOutputGP
+
+
 def fit_GP_MAP(*args, n_tries=15, theta0=None, method="L-BFGS-B",
                skip_failures=True, refit=False, **kwargs):
     """Fit one or more Gaussian Processes by attempting to minimize the
@@ -140,9 +142,10 @@ def fit_GP_MAP(*args, n_tries=15, theta0=None, method="L-BFGS-B",
     if len(args) == 1:
         gp = args[0]
         if isinstance(gp, MultiOutputGP):
-            gp =  _fit_MOGP_MAP(gp, n_tries, theta0, method, refit, **kwargs)
-        elif isinstance(gp, GaussianProcess):
+            gp = _fit_MOGP_MAP(gp, n_tries, theta0, method, refit, **kwargs)
+        elif isinstance(gp, GaussianProcessBase):
             gp = _fit_single_GP_MAP(gp, n_tries, theta0, method, **kwargs)
+
         else:
             raise TypeError("single arg to fit_GP_MAP must be a GaussianProcess or MultiOutputGP instance")
     elif len(args) < 2:
@@ -163,7 +166,7 @@ def fit_GP_MAP(*args, n_tries=15, theta0=None, method="L-BFGS-B",
             except AssertionError:
                 raise ValueError("Bad values for *args in fit_GP_MAP")
 
-    if isinstance(gp, GaussianProcess):
+    if isinstance(gp, GaussianProcessBase):
         if gp.theta is None:
             raise RuntimeError("GP fitting failed")
     else:
@@ -177,7 +180,7 @@ def fit_GP_MAP(*args, n_tries=15, theta0=None, method="L-BFGS-B",
 
     return gp
 
-        
+
 def _fit_single_GP_MAP(gp, n_tries=15, theta0=None, method='L-BFGS-B', **kwargs):
     """Fit hyperparameters using MAP for a single GP
 
@@ -187,7 +190,7 @@ def _fit_single_GP_MAP(gp, n_tries=15, theta0=None, method='L-BFGS-B', **kwargs)
 
     """
 
-    assert isinstance(gp, GaussianProcess)
+    assert isinstance(gp, GaussianProcessBase)
 
     n_tries = int(n_tries)
     assert n_tries > 0, "number of attempts must be positive"
@@ -281,6 +284,7 @@ def _fit_MOGP_MAP(gp, n_tries=15, theta0=None, method='L-BFGS-B',
 
     n_tries = int(n_tries)
 
+
     if refit:
         emulators_to_fit = gp.emulators
         indices_to_fit = list(range(len(gp.emulators)))
@@ -289,8 +293,9 @@ def _fit_MOGP_MAP(gp, n_tries=15, theta0=None, method='L-BFGS-B',
         indices_to_fit = gp.get_indices_not_fit()
         emulators_to_fit = gp.get_emulators_not_fit()
         thetavals = [ theta0[idx] for idx in indices_to_fit]
-    
-    if platform.system() == "Windows":
+
+    if platform.system() == "Windows" or gp.use_gpu:
+
         fit_MOGP = [fit_GP_MAP(emulator, n_tries=n_tries, theta0=t0, method=method, **kwargs)
                     for (emulator, t0) in zip(emulators_to_fit, thetavals)]
     else:
diff --git a/mogp_gpu/Makefile b/mogp_gpu/Makefile
new file mode 100644
index 00000000..bde99270
--- /dev/null
+++ b/mogp_gpu/Makefile
@@ -0,0 +1,48 @@
+# Makefile for standalone C++ tests of CUDA functionality.
+# Note that this Makefile is NOT used when building mogp_emulator.
+# (The logic in setup.py that builds the Extension takes care of that).
+# The targets in this Makefile are mainly standalone tests of low-level functionality.
+#
+# If running on CSD3, run the following commands to ensure dependencies are in place:
+#
+# module load anaconda  #(not needed here, but needed for mogp_emulator)
+# module load py-pybind11-2.2.4-gcc-5.4.0-tdtz6iq #(only needed for mogp_emulator)
+# module load cuda/11.2
+# module load gcc/7
+# module load eigen
+#
+
+CXX=g++
+CXXFLAGS=-std=c++14 -O3
+NVCCFLAGS=--compiler-options -Wall,-Wextra -arch=sm_60 --generate-code arch=compute_60,code=sm_60 --generate-code arch=compute_37,code=sm_37
+EXTRAFLAGS=-Xcompiler -fPIC
+
+CUDA_INC=/usr/local/cuda/include
+PYBIND_INC=$(shell python3 -m pybind11 --includes)
+
+all: lib/libgpgpu.so
+
+clean:
+	rm -rf obj lib bin
+
+lib/libgpgpu.so:  src/gp_gpu.cu obj/cov_gpu.o obj/util.o
+	mkdir -p lib
+	nvcc $(CXXFLAGS) $(NVCCFLAGS) $(EXTRAFLAGS) -shared $(PYBIND_INC) $< obj/cov_gpu.o obj/util.o -o $@ -lcusolver -lcublas -lcudart
+
+obj/%.o: src/%.cu src/%.hpp
+	mkdir -p obj
+	nvcc $(CXXFLAGS) $(NVCCFLAGS) $(EXTRAFLAGS) -c $< -o $@
+
+bin/test/gpu_basic: test/gpu_basic.cu obj/cov_gpu.o obj/util.o
+	mkdir -p bin/test
+	nvcc $(CXXFLAGS) $(NVCCFLAGS) $< obj/cov_gpu.o obj/util.o -o $@
+
+bin/test/gpu_test_utils: test/gpu_test_utils.cu obj/util.o
+	mkdir -p bin/test
+	nvcc $(CXXFLAGS) $(NVCCFLAGS) $< obj/util.o -o $@
+
+bin/test/gpu_cholesky: test/gpu_cholesky.cu
+	mkdir -p bin/test
+	nvcc $(CXXFLAGS) $(NVCCFLAGS) -lcusolver -lcublas $< -o $@
+
+.PHONY: all clean
diff --git a/mogp_gpu/src/cov_gpu.cu b/mogp_gpu/src/cov_gpu.cu
new file mode 100644
index 00000000..dbe2b20b
--- /dev/null
+++ b/mogp_gpu/src/cov_gpu.cu
@@ -0,0 +1,182 @@
+#include "cov_gpu.hpp"
+#include <stdio.h>
+
+// Covariance device function
+__device__ REAL cov_val_d(int Ninput, REAL *x_d, REAL *y_d, REAL *theta_d)
+{
+    REAL s = 0.0;
+    for (unsigned int i=0; i < Ninput; i++)
+    {
+        REAL d_i = x_d[i] - y_d[i];
+        s += d_i * d_i * exp(theta_d[i]);
+    }
+    return exp(-0.5 * s + theta_d[Ninput]);
+}
+
+////////////////////
+__global__ void cov_val_kernel(REAL *result_d, int Ninput, REAL *x_d,
+                               REAL *y_d, REAL *theta_d)
+{
+    *result_d = cov_val_d(Ninput, x_d, y_d, theta_d);
+}
+
+void cov_val_gpu(REAL *result_d, int Ninput, REAL *x_d, REAL *y_d,
+                 REAL *theta_d)
+{
+    cov_val_kernel<<<1,1>>>(result_d, Ninput, x_d, y_d, theta_d);
+}
+
+////////////////////
+__global__ void cov_diag_kernel(REAL *result_d, int N, int Ninput, REAL *xnew_d,
+                                REAL *xs, REAL *theta_d)
+{
+    int i = blockIdx.x * blockDim.x + threadIdx.x;
+    if (i < N) {
+        result_d[i] = cov_val_d(Ninput, xnew_d + Ninput * i, xs + Ninput * i,
+                                theta_d);
+    }
+}
+
+void cov_diag_gpu(REAL *result_d, int N, int Ninput, REAL *xnew_d, REAL *xs_d,
+                  REAL *theta_d)
+{
+    const int threads_per_block = 256;
+    cov_diag_kernel<<<10, threads_per_block>>>(
+        result_d, N, Ninput, xnew_d, xs_d, theta_d);
+}
+
+////////////////////
+__global__ void cov_all_kernel(REAL *result_d, int N, int Ninput, REAL *xnew_d,
+                               REAL *xs_d, REAL *theta_d)
+{
+    int i = blockIdx.x * blockDim.x + threadIdx.x;
+    if (i < N) {
+        result_d[i] = cov_val_d(Ninput, xnew_d, xs_d + Ninput * i, theta_d);
+    }
+}
+
+void cov_all_gpu(REAL *result_d, int N, int Ninput, REAL *xnew_d, REAL *xs_d,
+                 REAL *theta_d)
+{
+    const int threads_per_block = 256;
+    cov_all_kernel<<<10, threads_per_block>>>(
+        result_d, N, Ninput, xnew_d, xs_d, theta_d);
+}
+
+////////////////////
+__global__ void cov_batch_kernel(REAL *result_d, int Nnew, int N, int Ninput,
+                                 REAL *xsnew_d, REAL *xs_d, REAL *theta_d)
+{
+    int i = blockIdx.x * blockDim.x + threadIdx.x;
+    int j = blockIdx.y * blockDim.y + threadIdx.y;
+    if (i < N && j < Nnew)
+    {
+        result_d[j + Nnew * i] =
+            cov_val_d(Ninput, xsnew_d + Ninput * j, xs_d + Ninput * i, theta_d);
+    }
+}
+
+void cov_batch_gpu(REAL *result_d, int Nnew, int N, int Ninput, REAL *xsnew_d,
+                   REAL *xs_d, REAL *theta_d)
+{
+    dim3 threads_per_block(8, 32);
+    dim3 blocks(250, 625);
+    cov_batch_kernel<<<blocks, threads_per_block>>>(
+        result_d, Nnew, N, Ninput, xsnew_d, xs_d, theta_d);
+}
+
+
+////////////////////
+__device__ void cov_deriv_x(REAL *result_d, int Ninput,
+			    const REAL *x_d, const REAL *y_d,
+			    const REAL *theta_d)
+{
+    REAL s = 0.0;
+    for (unsigned int i=0; i < Ninput; i++)
+    {
+        REAL d_i = x_d[i] - y_d[i];
+        REAL a = d_i * exp(theta_d[i]);
+	result_d[i] = a;
+        s += d_i * a;
+    }
+    REAL c = -exp(-0.5*s + theta_d[Ninput]);
+    for (unsigned int i=0; i < Ninput; i++)
+    {
+        result_d[i] *= c;
+    }
+}
+
+__global__ void cov_deriv_x_batch_kernel(
+    REAL *result_d, int Ninput, int Nx, int Ny, const REAL *xs_d,
+    const REAL *ys_d, const REAL *theta_d)
+{
+    int i = blockIdx.x * blockDim.x + threadIdx.x;
+    int j = blockIdx.y * blockDim.y + threadIdx.y;
+    if (i < Nx && j < Ny)
+    {
+        cov_deriv_x(result_d + Ninput * (Nx * j + i),
+		    Ninput, xs_d + Ninput * i, ys_d + Ninput * j,
+		    theta_d);
+    }
+}
+
+void cov_deriv_x_batch_gpu(
+    REAL *result_d, int Ninput, int Nx, int Ny, const REAL *xs_d,
+    const REAL *ys_d, const REAL *theta_d)
+{
+    const int Bx = 16, By = 16;
+    dim3 threads_per_block(Bx, By);
+    dim3 blocks((Nx + Bx - 1)/Bx, (Ny + By - 1)/By);
+    cov_deriv_x_batch_kernel<<<blocks, threads_per_block>>>(
+        result_d, Ninput, Nx, Ny, xs_d, ys_d, theta_d);
+}
+
+
+////////////////////
+__device__ void cov_deriv_theta(REAL *result_d, int result_stride, int Ninput,
+                          const REAL *x_d, const REAL *y_d,
+                          const REAL *theta_d)
+{
+    REAL s = 0.0;
+    REAL exp_thetaN = exp(theta_d[Ninput]);
+    for (unsigned int i=0; i < Ninput; i++)
+    {
+        REAL d_i = x_d[i] - y_d[i];
+        REAL a = d_i * d_i * exp(theta_d[i]);
+        s += a;
+        result_d[i * result_stride] = -0.5 * a;
+    }
+    s = exp_thetaN * exp(-0.5 * s);
+    for (unsigned int i=0; i < Ninput; i++)
+    {
+        result_d[i * result_stride] *= s;
+    }
+    result_d[Ninput * result_stride] = s;
+}
+
+__global__ void cov_deriv_theta_batch_kernel(
+    REAL *result_d, int Ninput, int Nx, int Ny, const REAL *xs_d,
+    const REAL *ys_d, const REAL *theta_d)
+{
+    int i = blockIdx.x * blockDim.x + threadIdx.x;
+    int j = blockIdx.y * blockDim.y + threadIdx.y;
+    const int result_stride = Nx * Ny;
+    if (i < Nx && j < Ny)
+    {
+        cov_deriv_theta(result_d + Ny * i + j,
+			result_stride,
+			Ninput, xs_d + Ninput * j, ys_d + Ninput * i,
+			theta_d);
+    }
+}
+
+void cov_deriv_theta_batch_gpu(
+    REAL *result_d, int Ninput, int Nx, int Ny, const REAL *xs_d,
+    const REAL *ys_d, const REAL *theta_d)
+{
+    const int Bx = 16, By = 16;
+    dim3 threads_per_block(Bx, By);
+    dim3 blocks((Nx + Bx - 1)/Bx, (Ny + By - 1)/By);
+    cov_deriv_theta_batch_kernel<<<blocks, threads_per_block>>>(
+        result_d, Ninput, Nx, Ny, xs_d, ys_d, theta_d);
+}
diff --git a/mogp_gpu/src/cov_gpu.hpp b/mogp_gpu/src/cov_gpu.hpp
new file mode 100644
index 00000000..d1460b0a
--- /dev/null
+++ b/mogp_gpu/src/cov_gpu.hpp
@@ -0,0 +1,37 @@
+#ifndef COV_GPU_HPP
+#define COV_GPU_HPP
+
+typedef double REAL;
+#define CUBLASDOT cublasDdot
+#define CUBLASGEMV cublasDgemv
+
+// Squared exponential covariance function
+
+// Performs a single evaluation of the covariance function, for inputs
+// x and y; that is, Cov(x(:), y(:); theta).
+void cov_val_gpu(REAL *result_d, int Ninput, REAL *x_d, REAL *y_d,
+                 REAL *theta_d);
+
+// Computes the vector k_i = Cov(xnew(:), xs(i,:); theta).
+void cov_all_gpu(REAL *result_d, int N, int Ninput, REAL *xnew_d, REAL *xs_d,
+                 REAL *theta_d);
+
+// Computes the vector k_i = Cov(xnew(i,:), xs(i,:); theta).
+void cov_diag_gpu(REAL *result_d, int N, int Ninput, REAL *xnew_d, REAL *xs_d,
+                  REAL *theta_d);
+
+// Computes the submatrix K_ij = Cov(xsnew(i,:), xs(j,:); theta)
+void cov_batch_gpu(REAL *result_d, int Nnew, int N, int Ninput, REAL *xsnew_d,
+                   REAL *xs_d, REAL *theta_d);
+
+// Computes the matrix (dC/dtheta)(xs(i,:), ys(j,:); theta)
+void cov_deriv_theta_batch_gpu(REAL *result_d, int Ninput, int Nx, int Ny,
+			       const REAL *xs_d, const REAL *ys_d,
+			       const REAL *theta_d);
+
+// Computes the matrix (dC/dx)(xs(i,:), ys(j,:); theta)
+void cov_deriv_x_batch_gpu(REAL *result_d, int Ninput, int Nx, int Ny,
+			   const REAL *xs_d, const REAL *ys_d,
+			   const REAL *theta_d);
+
+#endif
diff --git a/mogp_gpu/src/gp_gpu.cu b/mogp_gpu/src/gp_gpu.cu
new file mode 100644
index 00000000..b39cc5bd
--- /dev/null
+++ b/mogp_gpu/src/gp_gpu.cu
@@ -0,0 +1,238 @@
+#include "gp_gpu.hpp"
+#include "pybind11/pybind11.h"
+#include "pybind11/stl.h"
+#include "pybind11/eigen.h"
+
+namespace py = pybind11;
+
+bool have_compatible_device(void)
+{
+    int device_count = 0;
+    cudaError_t err = cudaGetDeviceCount(&device_count);
+
+    return (err == cudaSuccess && device_count > 0);
+
+}
+
+
+PYBIND11_MODULE(libgpgpu, m) {
+    py::class_<DenseGP_GPU>(m, "DenseGP_GPU")
+        .def(py::init<mat_ref, vec_ref, unsigned int>())
+
+////////////////////////////////////////
+        .def("n", &DenseGP_GPU::get_n,
+             "Returns the number of training examples used for the emulator")
+
+////////////////////////////////////////
+        .def("D", &DenseGP_GPU::get_D,
+             "Returns the dimension of the emulator inputs")
+
+////////////////////////////////////////
+        .def("inputs", &DenseGP_GPU::get_inputs,
+             "Returns the array of inputs points used by the emulator, with "
+             "shape ``(n, D)``")
+
+////////////////////////////////////////
+        .def("targets", &DenseGP_GPU::get_targets,
+             "Returns the array of targets used by the emulator, with shape "
+             "``(n,)``")
+
+////////////////////////////////////////
+        .def("n_params", &DenseGP_GPU::get_n_params,
+             "Returns the number of hyperparameters required for the emulator.")
+
+////////////////////////////////////////
+        .def("theta_fit_status", &DenseGP_GPU::get_theta_fit_status,
+             "Has the emulator been fit?  If `false`, "
+             ":func:``GaussianProcessGPU.theta`` will be ``None``")
+
+////////////////////////////////////////
+        .def("reset_theta_fit_status", &DenseGP_GPU::reset_theta_fit_status,
+             "Marks this emulator as having not been fit: "
+             ":func:`theta_fit_status` will be `false`")
+
+////////////////////////////////////////
+        .def("get_theta", &DenseGP_GPU::get_theta,
+             "Returns the emulator hyperparameters")
+
+////////////////////////////////////////
+        .def("predict", &DenseGP_GPU::predict,
+             R"(Single point predictive mean.
+
+:param testing: The input point
+
+:returns: The predicted mean
+
+Currently unused by :class`GaussianProcessGPU`, but useful for testing)",
+             py::arg("testing"))
+
+////////////////////////////////////////
+        .def("predict_variance", &DenseGP_GPU::predict_variance,
+             R"(Single point combined predictive mean and variance.
+
+:param testing: The input point to predict, with shape `(Nbatch, D)`
+                (`Nbatch` is determined from the shape of the input).
+                `Nbatch` must be less than `testing_size`.
+:param var: (Output) The variance prediction is written to `var[0]`
+
+:returns: The predicted mean
+
+Currently unused by :class`GaussianProcessGPU`, but useful for testing.
+)",
+             py::arg("testing"),
+             py::arg("var"))
+
+////////////////////////////////////////
+        .def("predict_batch", &DenseGP_GPU::predict_batch,
+             R"(Batched predictive means.
+
+:param testing: The input point to predict, with shape `(Nbatch, D)`
+                (`Nbatch` is determined from the shape of the input).
+                `Nbatch` must be less than `testing_size`.
+:param result: (Output) The predicted mean at each `testing` point.
+               Shape `(Nbatch, D)`.
+
+:returns: None
+)",
+             py::arg("testing"),
+             py::arg("result"))
+
+////////////////////////////////////////
+        .def("predict_variance_batch", &DenseGP_GPU::predict_variance_batch,
+             R"(Batched predictive means and variances.
+
+:param testing: The input point to predict, with shape `(Nbatch, D)`
+                (`Nbatch` is determined from the shape of the input).
+                `Nbatch` must be less than `testing_size`.
+:param mean: (Output) The predicted mean at each `testing` point.
+             Shape `(Nbatch, D)`.
+:param var: (Output) The predicted variance at each `testing` point.
+            Shape `(Nbatch, D)`.
+
+:returns: None
+)",
+             py::arg("testing"),
+             py::arg("mean"),
+             py::arg("var"))
+
+////////////////////////////////////////
+        .def("predict_deriv", &DenseGP_GPU::predict_deriv,
+             R"(Batched prediction of derivatives.
+
+:param testing: The input point to predict, with shape `(Nbatch, D)`
+                (`Nbatch` is determined from the shape of the input).
+                `Nbatch` must be less than `testing_size`.
+:param result: (Output) The predicted gradient of the mean with
+               respect to the inputs.  Shape `(Nbatch, D)`.
+
+:returns: None
+)",
+             py::arg("testing"),
+             py::arg("result"))
+
+////////////////////////////////////////
+
+        .def("fit", &DenseGP_GPU::fit,
+             R"(Sets theta and the nugget parameters, and updates the internal
+state of the emulator accordingly.
+
+:param theta: Values for the hyperparameters (length ``n_params``)
+:param nugget: The interpretation of the nugget: adaptive, fixed or fit
+:param nugget_size: The (initial) nugget value
+
+:returns: None
+)",
+             py::arg("theta"),
+             py::arg("nugget"),
+             py::arg("nugget_size")
+            )
+
+////////////////////////////////////////
+        .def("get_K", &DenseGP_GPU::get_K,
+             R"(The current covariance matrix
+
+:param K_h: (Output) The matrix, of shape ``(n, n)``.
+
+:returns: None
+)",
+             py::arg("K_h"))
+
+////////////////////////////////////////
+        .def("get_invQ", &DenseGP_GPU::get_invQ,
+             R"(The inverse of the current covariance matrix
+
+:param invQ_h: (Output) The matrix, of shape ``(n, n)``.
+
+:returns: None
+)",
+             py::arg("invQ_h"))
+
+////////////////////////////////////////
+        .def("get_invQt", &DenseGP_GPU::get_invQt,
+             R"(The inverse covariance matrix times the target vector.
+
+:param invQt_h: (Output) The resulting vector, of shape ``(n,)``
+
+:returns: None
+)",
+             py::arg("invQt_h"))
+
+////////////////////////////////////////
+        .def("get_logpost", &DenseGP_GPU::get_logpost,
+	     "Return the log posterior"
+            )
+
+////////////////////////////////////////
+        .def("get_jitter", &DenseGP_GPU::get_jitter,
+             "The value of jitter used for an adaptive nugget.")
+
+////////////////////////////////////////
+        .def("get_cholesky_lower", &DenseGP_GPU::get_cholesky_lower,
+             R"(Return the (lower) Cholesky factor of the covariance matrix.
+
+This is the output of (cu)BLAS dpotrf.  The lower triangular entries correspond
+to the factor; the upper triangular entries contain junk values.
+
+:param result: (Output) the resulting Cholesky factor.  Shape ``(n,n)``.
+               Lower triangular entries set.
+
+:returns: None
+ )",
+             py::arg("result"))
+
+
+////////////////////////////////////////
+        .def("logpost_deriv", &DenseGP_GPU::logpost_deriv,
+             R"(Calculate the partial derivatives of the negative log-posterior
+likelihood of the hyperparameters, from the current state of the emulator.)
+
+:param result: (Output) the resulting derivatives, shape ``(D+1,)``
+
+:returns: `None`
+)",
+             py::arg("result"));
+
+
+
+    py::enum_<nugget_type>(m, "nugget_type")
+        .value("adaptive", NUG_ADAPTIVE)
+        .value("fixed", NUG_FIXED)
+        .value("fit", NUG_FIT);
+
+    m.def("have_compatible_device", &have_compatible_device);
+
+    m.doc() = R"(
+The libgpgpu library
+--------------------
+
+.. currentmodule:: libgpgpu
+
+The members of this class are
+
+.. autoclass:: libgpgpu.DenseGP_GPU
+    :members:
+
+)";
+
+    // Include have_compatible_device in the docstring
+}
diff --git a/mogp_gpu/src/gp_gpu.hpp b/mogp_gpu/src/gp_gpu.hpp
new file mode 100644
index 00000000..07d87331
--- /dev/null
+++ b/mogp_gpu/src/gp_gpu.hpp
@@ -0,0 +1,815 @@
+#ifndef GP_GPU_HPP
+#define GP_GPU_HPP
+
+/*
+This file contains the C++ implementation of the Gaussian Process class.
+(The corresponding .cu file contains the pybind11 bindings allowing the functions
+to be called from Python).
+
+The key methods of the DenseGP_GPU class are:
+  fit:  update the theta hyperparameters
+  predict_batch:  make a prediction on a set of testing points
+  predict_variance_batch: make a prediction on a set of testing points, including variance
+  predict_deriv: get the derivative of prediction on a set of testing points.
+These in turn use CUDA kernels defined in the file cov_gpu.cu
+*/
+
+#include <iostream>
+
+#include <algorithm>
+#include <string>
+#include <sstream>
+#include <assert.h>
+#include <stdexcept>
+#include <cuda_runtime.h>
+#include <cublas_v2.h>
+#include <cusolverDn.h>
+#include <thrust/device_ptr.h>
+#include <thrust/device_malloc.h>
+#include <thrust/device_free.h>
+#include <thrust/device_vector.h>
+#include <thrust/functional.h>
+#include <thrust/transform_reduce.h>
+#include <thrust/copy.h>
+#include <cub/cub.cuh>
+
+#include <Eigen/Dense>
+
+#include "strided_range.hpp"
+#include "cov_gpu.hpp"
+#include "util.hpp"
+
+#define WARP_SIZE 32
+#define FULL_MASK 0xffffffff
+
+#define USE_SHUFFLE_SUM_IMPL 0
+
+// ----------------------------------------------
+// ----   Utility functions ---------------------
+
+/// Extract the raw pointer from a device vector
+template <typename T>
+T *dev_ptr(thrust::device_vector<T>& dv)
+{
+    return dv.data().get();
+}
+
+
+/// Fail if a recent cusolver call did not succeed
+inline void check_cusolver_status(cusolverStatus_t status, int info_h)
+{
+    if (status || info_h) {
+	std::string msg;
+	std::stringstream smsg(msg);
+	smsg << "Error in potrf: return code " << status << ", info " << info_h;
+	throw std::runtime_error(smsg.str());
+    }
+}
+
+/// Can a usable CUDA capable device be found?
+bool have_compatible_device(void);
+
+// ----------------------------------------------
+
+
+// --------------------------------------------------------
+// CUDA kernel for summing log diagonal elements of a matrix.
+// Two methods implemented - warp reduction or cub::DeviceReduce
+
+#if USE_SHUFFLE_SUM_IMPL
+// Implementation of sum_log_diag using warp reduction,
+// based on https://devblogs.nvidia.com/faster-parallel-reductions-kepler/
+//
+// Will fail for N > 1024
+
+
+__inline__ __device__
+double warp_reduce_sum(double val)
+{
+    for (int offset = WARP_SIZE/2; offset > 0; offset /= 2)
+        val += __shfl_down_sync(FULL_MASK, val, offset);
+    return val;
+}
+
+__global__
+void sum_log_diag_kernel(int N, double *A, double *result)
+{
+    int i = threadIdx.x;
+
+    static __shared__ double work[WARP_SIZE];
+    double log_diag = 0.0;
+
+    if (i < WARP_SIZE) work[i] = 0.0;
+    if (i < N) {
+        log_diag = log(A[i * (N+1)]);
+    }
+
+    int laneIdx = i % WARP_SIZE;
+    int warpIdx = i / WARP_SIZE;
+
+    __syncthreads();
+    log_diag = warp_reduce_sum(log_diag);
+
+    if (laneIdx == 0) work[warpIdx] = log_diag;
+
+    __syncthreads();
+
+    log_diag = work[laneIdx];
+
+    if (warpIdx == 0) log_diag = warp_reduce_sum(log_diag);
+
+    if (i == 0) *result = 2.0 * log_diag;
+}
+
+//
+// The unused arguments are needed for the CUB implementation
+void sum_log_diag(int N, double *A, double *result, double *, size_t)
+{
+    // The number of thread blocks *must* be 1
+    // N can be at most 1024 (WARP_SIZE**2)
+    int thread_blocks = 1;
+    int threads_per_block = N;
+    sum_log_diag_kernel<<<thread_blocks, threads_per_block>>>(N, A, result);
+}
+
+#else
+// ----------------------------------------
+
+// Implementation of sum_log_diag using cub::DeviceReduce
+
+struct LogSq : public thrust::unary_function<double, double>
+{
+    __host__ __device__ double operator()(double x) const { return 2.0 * log(x); }
+};
+
+void sum_log_diag(int N, double *A, double *result, double *work, size_t work_size)
+{
+    auto transform_it = thrust::make_transform_iterator(A, LogSq());
+    auto sr = make_strided_range(transform_it, transform_it + N*N, N+1);
+
+    cub::DeviceReduce::Sum(work, work_size, sr.begin(), result, N);
+}
+
+// Implementation of trace using cub::DeviceReduce
+
+void trace(int N, double *A, double *result, double *work, size_t work_size)
+{
+    auto sr = make_strided_range(A, A + N*N, N+1);
+    cub::DeviceReduce::Sum(work, work_size, sr.begin(), result, N);
+}
+
+// ----------------------------------------
+#endif // USE_SHUFFLE_SUM_IMPL
+
+
+// ----------------------------------------
+// ------- Some useful typedefs for vectors and matrices
+
+typedef int obs_kind;
+typedef typename Eigen::Matrix<REAL, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> mat;
+typedef typename Eigen::Ref<Eigen::Matrix<REAL, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> > mat_ref;
+typedef typename Eigen::Matrix<REAL, Eigen::Dynamic, 1> vec;
+typedef typename Eigen::Ref<Eigen::Matrix<REAL, Eigen::Dynamic, 1> > vec_ref;
+// ----------------------------------------
+// enum to allow the python code to select the type of "nugget"
+enum nugget_type {NUG_ADAPTIVE, NUG_FIT, NUG_FIXED};
+
+
+
+//////////////////////////////////////////////////////////////////////////////////////
+////////////////////////// The Gaussian Process  class ///////////////////////////////
+//////////////////////////////////////////////////////////////////////////////////////
+
+// By convention, members ending "_d" are allocated on the device
+class DenseGP_GPU {
+    // batch size (for allocation)
+    unsigned int testing_size;
+
+    // Number of training points, dimension of input
+    unsigned int n, D;
+
+    // inputs
+    mat inputs;
+
+    // targets
+    vec targets;
+
+    // is the nugget adapted, fixed, or fitted?
+    nugget_type nug_type;
+
+    // flag for whether we have fit hyperparameters
+    bool theta_fitted;
+
+    // handle for CUBLAS calls
+    cublasHandle_t cublasHandle;
+
+    // handle for cuSOLVER calls
+    cusolverDnHandle_t cusolverHandle;
+
+    // covariance matrix on device
+    thrust::device_vector<REAL> K_d;
+
+    // inverse covariance matrix on device
+    thrust::device_vector<REAL> invQ_d;
+
+    // lower triangular Cholesky factor on device
+    thrust::device_vector<REAL> chol_lower_d;
+
+    // log determinant of covariance matrix (on host)
+    double logdetC;
+
+    // adaptive nugget jitter
+    double jitter;
+
+    // hyperparameters on the device
+    thrust::device_vector<REAL> theta_d;
+
+    // inputs, on the device, row major order
+    thrust::device_vector<REAL> inputs_d;
+
+    // targets, on the device
+    thrust::device_vector<REAL> targets_d;
+
+    // precomputed product, used in the prediction
+    thrust::device_vector<REAL> invQt_d;
+
+    // preallocated work array (length n) on the CUDA device
+    thrust::device_vector<REAL> work_d;
+
+    // buffer for Cholesky factorization
+    thrust::device_vector<REAL> potrf_buffer_d;
+
+    // buffer for cub::DeviceReduce::Sum
+    thrust::device_vector<REAL> sum_buffer_d;
+
+    // size of sum_buffer_d
+    size_t sum_buffer_size_bytes;
+
+    // more work arrays
+    thrust::device_vector<REAL> testing_d, work_mat_d, result_d, kappa_d, invQk_d;
+
+public:
+    int get_n(void) const
+    {
+        return n;
+    }
+
+    int get_D(void) const
+    {
+        return D;
+    }
+
+    mat get_inputs(void) const
+    {
+        return inputs;
+    }
+
+    vec get_targets(void) const
+    {
+        return targets;
+    }
+
+    int get_n_params(void) const
+    {
+        return D + 2;
+    }
+
+    void get_theta(vec_ref theta)
+    {
+        thrust::copy(theta_d.begin(), theta_d.end(), theta.data());
+    }
+
+    double get_jitter(void) const
+    {
+        return jitter;
+    }
+
+    // make a single prediction (mainly for testing - most use-cases will use predict_batch or predict_deriv_batch)
+    double predict(mat_ref testing)
+    {
+        // On entry: the number of points to predict (number of rows
+        // of testing) is assumed to be D
+
+        thrust::device_vector<REAL> testing_d(testing.data(), testing.data() + D);
+        cov_all_gpu(dev_ptr(work_d), n, D, dev_ptr(testing_d), dev_ptr(inputs_d),
+                    dev_ptr(theta_d));
+
+        REAL result = std::numeric_limits<double>::quiet_NaN();
+        CUBLASDOT(cublasHandle, n, dev_ptr(work_d), 1, dev_ptr(invQt_d), 1,
+                  &result);
+
+        return double(result);
+    }
+
+    // variance of a single prediction (mainly for testing - most use-cases will use predict_variance_batch)
+    double predict_variance(mat_ref testing, vec_ref var)
+    {
+        if (var.size() < testing.rows()) {
+            throw std::runtime_error("predict_variance: the result buffer passed was "
+                                     "too small to hold the variance");
+        }
+
+        // value prediction
+        thrust::device_vector<REAL> testing_d(testing.data(), testing.data() + D);
+        cov_all_gpu(dev_ptr(work_d), n, D, dev_ptr(testing_d), dev_ptr(inputs_d),
+                    dev_ptr(theta_d));
+
+        REAL result = std::numeric_limits<double>::quiet_NaN();
+        CUBLASDOT(cublasHandle, n, dev_ptr(work_d), 1, dev_ptr(invQt_d), 1,
+                  &result);
+
+        // variance prediction
+        thrust::device_vector<REAL> kappa_d(1);
+        cov_val_gpu(dev_ptr(kappa_d), D, dev_ptr(testing_d), dev_ptr(testing_d),
+                    dev_ptr(theta_d));
+
+        double zero(0.0);
+        double one(1.0);
+        thrust::device_vector<REAL> invQk_d(n);
+
+        cublasDgemv(cublasHandle, CUBLAS_OP_N, n, n, &one, dev_ptr(invQ_d), n,
+                    dev_ptr(work_d), 1, &zero, dev_ptr(invQk_d), 1);
+
+        CUBLASDOT(cublasHandle, n, dev_ptr(work_d), 1, dev_ptr(invQt_d),
+                  1, &result);
+        CUBLASDOT(cublasHandle, n, dev_ptr(work_d), 1, dev_ptr(invQk_d),
+                  1, var.data());
+
+        double kappa;
+        thrust::copy(kappa_d.begin(), kappa_d.end(), &kappa);
+
+        cudaDeviceSynchronize();
+
+        var = kappa - var.array();
+
+        return REAL(result);
+    }
+
+    // Use the GP emulator to calculate a prediction on testing points, without calculating variance or derivative
+    void predict_batch(mat_ref testing, vec_ref result)
+    {
+        // Assumes on entry that testing has shape (Nbatch, D), and that result
+        // contains space for Nbatch result values
+
+        int Nbatch = testing.rows();
+
+        if (result.size() < testing.rows()) {
+            throw std::runtime_error("predict_batch: the result buffer passed was "
+                                     "too small to hold the result");
+        }
+
+        if (Nbatch > testing_size) {
+            throw std::runtime_error("predict_variance_batch: More test points were passed "
+                                     "than the maximum batch size");
+        }
+
+        REAL zero(0.0);
+        REAL one(1.0);
+        thrust::device_vector<REAL> testing_d(testing.data(), testing.data() + Nbatch * D);
+        thrust::device_vector<REAL> result_d(Nbatch);
+
+        cov_batch_gpu(dev_ptr(work_mat_d), Nbatch, n, D,
+                      dev_ptr(testing_d), dev_ptr(inputs_d), dev_ptr(theta_d));
+
+        cublasStatus_t status =
+            cublasDgemv(cublasHandle, CUBLAS_OP_N, Nbatch, n, &one,
+                        dev_ptr(work_mat_d), Nbatch, dev_ptr(invQt_d), 1, &zero,
+                        dev_ptr(result_d), 1);
+
+        cudaDeviceSynchronize();
+
+        thrust::copy(result_d.begin(), result_d.end(), result.data());
+    }
+
+    // Use the GP emulator to calculate a prediction on testing points, also calculating variance
+    void predict_variance_batch(mat_ref testing, vec_ref mean, vec_ref var)
+    {
+        REAL zero(0.0);
+        REAL one(1.0);
+        REAL minus_one(-1.0);
+        int Nbatch = testing.rows();
+
+        if (var.size() < testing.rows() || mean.size() < testing.rows()) {
+            throw std::runtime_error("predict_variance_batch: The result buffer passed was "
+                                     "too small to hold the variance");
+        }
+
+        if (Nbatch > testing_size) {
+            throw std::runtime_error("predict_variance_batch: More test points were passed "
+                                     "than the maximum batch size");
+        }
+
+
+        thrust::device_vector<REAL> testing_d(testing.data(), testing.data() + Nbatch * D);
+        thrust::device_vector<REAL> mean_d(Nbatch);
+
+        // compute predictive means for the batch
+        cov_batch_gpu(dev_ptr(work_mat_d), Nbatch, n, D, dev_ptr(testing_d),
+                      dev_ptr(inputs_d), dev_ptr(theta_d));
+
+        cublasStatus_t status =
+            cublasDgemv(cublasHandle, CUBLAS_OP_N, Nbatch, n, &one,
+                        dev_ptr(work_mat_d), Nbatch, dev_ptr(invQt_d), 1, &zero,
+                        dev_ptr(mean_d), 1);
+
+        // compute predictive variances for the batch
+        cov_diag_gpu(dev_ptr(kappa_d), Nbatch, D, dev_ptr(testing_d),
+                     dev_ptr(testing_d), dev_ptr(theta_d));
+
+        cublasDgemm(cublasHandle, CUBLAS_OP_N, CUBLAS_OP_T,
+                    n,
+                    Nbatch,
+                    n,
+                    &one,
+                    dev_ptr(invQ_d), n,
+                    dev_ptr(work_mat_d), Nbatch,
+                    &zero,
+                    dev_ptr(invQk_d), n);
+
+        // result accumulated into 'kappa'
+        status = cublasDgemmStridedBatched(
+            cublasHandle, CUBLAS_OP_N, CUBLAS_OP_N,
+            1, // m
+            1, // n
+            n, // k
+            // A (m x k), B (k x n), C (m x n)
+            &minus_one, // alpha
+            dev_ptr(work_mat_d), Nbatch, 1, // A, lda, strideA
+            dev_ptr(invQk_d), n, n, // B, ldb, strideB (= covariances "k")
+            &one,
+            dev_ptr(kappa_d), 1, 1, // C, ldc, strideC
+            Nbatch);
+
+        cudaDeviceSynchronize();
+
+        // copy back means
+        thrust::copy(mean_d.begin(), mean_d.end(), mean.data());
+
+        // copy back variances
+        thrust::copy(kappa_d.begin(), kappa_d.begin() + Nbatch, var.data());
+    }
+
+    // Use the GP emulator to calculate derivative of  prediction on testing points
+    void predict_deriv(mat_ref testing, mat_ref result)
+    {
+        int Nbatch = testing.rows();
+
+        if (result.rows() < testing.rows() || result.cols() != D) {
+            throw std::runtime_error("predict_deriv: the result buffer passed was "
+                                     "the wrong shape to hold the result");
+        }
+
+        if (Nbatch > testing_size) {
+            throw std::runtime_error("predict_variance_batch: More test points were passed "
+                                     "than the maximum batch size");
+        }
+
+        REAL zero(0.0);
+        REAL one(1.0);
+        thrust::device_vector<REAL> testing_d(testing.data(), testing.data() + Nbatch * D);
+        thrust::device_vector<REAL> result_d(Nbatch*D);
+
+        cov_deriv_x_batch_gpu(dev_ptr(work_mat_d), D, Nbatch, n,
+			      dev_ptr(testing_d), dev_ptr(inputs_d), dev_ptr(theta_d));
+
+        cublasStatus_t status =
+            cublasDgemv(cublasHandle, CUBLAS_OP_N,
+                        D * Nbatch, n, // nrows, ncols
+                        &one, dev_ptr(work_mat_d), D * Nbatch, // alpha, A, lda
+                        dev_ptr(invQt_d), 1, // x, incx
+                        &zero, dev_ptr(result_d), 1); // beta, y, incy
+
+        cudaDeviceSynchronize();
+
+        thrust::copy(result_d.begin(), result_d.end(), result.data());
+    }
+
+    // perform Cholesky factorization of the matrix currently stored in work_mat_d
+    int calc_cholesky_factors(void)
+    {
+        // On entry: assumes that work_mat_d contains the inverse covariance matrix invQ_d
+
+        thrust::device_vector<int> info_d(1);
+        int info_h;
+        cusolverStatus_t status;
+
+	// compute Cholesky factors
+	status = cusolverDnDpotrf(cusolverHandle, CUBLAS_FILL_MODE_LOWER, n,
+				  dev_ptr(work_mat_d), n, dev_ptr(potrf_buffer_d),
+				  potrf_buffer_d.size(), dev_ptr(info_d));
+
+	thrust::copy(info_d.begin(), info_d.end(), &info_h);
+
+	if (status != CUSOLVER_STATUS_SUCCESS || info_h < 0) {
+            std::string msg;
+            std::stringstream smsg(msg);
+            smsg << "Error in potrf: return code " << status << ", info " << info_h;
+            throw std::runtime_error(smsg.str());
+
+	}
+	return info_h;
+    }
+
+    // return the lower triangular matrix (actually it will be the whole matrix
+    // but only the lower triangle will be correct!)
+    void get_cholesky_lower(mat_ref result)
+    {
+        thrust::copy(chol_lower_d.begin(), chol_lower_d.end(), result.data());
+    }
+
+    // Update the hyperparameters, and K, invQ and invQt which depend on them
+    void fit(vec_ref theta, nugget_type nugget, double nugget_size=0.0)
+    {
+
+        nug_type = nugget;
+        thrust::copy(theta.data(), theta.data() + D + 2, theta_d.begin());
+
+	// calculate covariance matrix, put result into K_d
+        cov_batch_gpu(dev_ptr(K_d), n, n, D, dev_ptr(inputs_d),
+                      dev_ptr(inputs_d), dev_ptr(theta_d));
+
+	/// copy the covariance matrix K_d into work_mat_d
+	thrust::copy(K_d.begin(), K_d.end(), work_mat_d.begin());
+
+        thrust::device_vector<int> info_d(1);
+        int info_h;
+        cusolverStatus_t status;
+	int factorisation_status;
+
+	// for adaptive nugget start with a nugget of zero and increase by small amount
+	// until we find a value where factorization succeeds.
+	if (nug_type == NUG_ADAPTIVE) {
+	    double mean_diag;
+	    const int max_tries = 5;
+	    int itry = 0;
+	    while (itry < max_tries) {
+		// K_d holds the covariance matrix - work with a copy
+		// in work_mat_d, in case the factorization fails
+		thrust::copy(K_d.begin(), K_d.end(), work_mat_d.begin());
+
+		// if the first attempt at factorization failed, add a
+		// small term to the diagonal, increasing each iteration
+		// until the factorization succeeds
+		if (itry >= 1) {
+		    if (itry == 1) {
+			// find mean of (absolute) diagonal elements (diagonal
+			// elements should all be positive)
+			cublasDasum(cublasHandle, n, dev_ptr(K_d), n+1, &mean_diag);
+			mean_diag /= n;
+			jitter = 1e-6 * mean_diag;
+		    }
+		    add_diagonal(n, jitter, dev_ptr(work_mat_d));
+		}
+
+		factorisation_status = calc_cholesky_factors();
+		if (factorisation_status == 0) {
+                    break;
+		}
+		jitter *= 10;
+		itry++;
+	    }
+
+	    // if none of the factorization attempts succeeded:
+	    if (itry == max_tries) {
+		std::string msg;
+		std::stringstream smsg(msg);
+		smsg << "All attempts at factorization failed. Last return code " << factorisation_status;
+		throw std::runtime_error(smsg.str());
+	    }
+	// for fixed nugget, add "nugget_size" to the diagonal of the matrix.
+	} else if (nug_type == NUG_FIXED) {
+	    add_diagonal(n, nugget_size, dev_ptr(work_mat_d));
+
+	    factorisation_status = calc_cholesky_factors();
+	    if (factorisation_status != 0) {
+                throw std::runtime_error("Unable to factorize matrix using fixed nugget");
+	    }
+
+	} else if (nug_type == NUG_FIT) {
+	    // set to exp(last-element-of-theta)
+	    jitter = exp( theta(theta.size()-1) );
+
+	    add_diagonal(n, jitter, dev_ptr(work_mat_d));
+	    factorisation_status = calc_cholesky_factors();
+	    if (factorisation_status != 0) {
+                throw std::runtime_error("Unable to factorize matrix using fitted nugget");
+	    }
+	} else throw std::runtime_error("Unrecognized nugget_type");
+
+
+        // get the inverse covariance matrix invQ by solving the system of linear eqns
+	//    work_mat_d . invQ_d = I
+	// where work_mat_d is holding the current covariance matrix.
+
+        identity_device(n, dev_ptr(invQ_d));
+
+        status = cusolverDnDpotrs(cusolverHandle, CUBLAS_FILL_MODE_LOWER, n, n,
+                                  dev_ptr(work_mat_d), n, dev_ptr(invQ_d), n,
+                                  dev_ptr(info_d));
+
+        thrust::copy(info_d.begin(), info_d.end(), &info_h);
+        check_cusolver_status(status, info_h);
+
+
+        // invQt - product of inverse covariance matrix with the target values
+        thrust::copy(targets_d.begin(), targets_d.end(), invQt_d.begin());
+        status = cusolverDnDpotrs(cusolverHandle, CUBLAS_FILL_MODE_LOWER, n, 1,
+                                  dev_ptr(work_mat_d), n, dev_ptr(invQt_d), n,
+                                  dev_ptr(info_d));
+
+        thrust::copy(info_d.begin(), info_d.end(), &info_h);
+        check_cusolver_status(status, info_h);
+
+        // logdetC - sum the log of the diagonal elements of the Cholesky factor of covariance matrix (in work_mat_d)
+        thrust::device_vector<double> logdetC_d(1);
+
+        sum_log_diag(n, dev_ptr(work_mat_d), dev_ptr(logdetC_d), dev_ptr(sum_buffer_d), sum_buffer_size_bytes);
+
+        thrust::copy(logdetC_d.begin(), logdetC_d.end(), &logdetC);
+
+	//copy work_mat_d into the lower triangular Cholesky factor
+	thrust::copy(work_mat_d.begin(), work_mat_d.begin()+n*n, chol_lower_d.begin());
+
+	//set the flag to say we have fitted theta
+	theta_fitted = true;
+    }
+
+    bool get_theta_fit_status(void)
+    {
+        return theta_fitted;
+    }
+
+    void reset_theta_fit_status(void)
+    {
+        theta_fitted = false;
+    }
+
+    void get_K(mat_ref K_h)
+    {
+        thrust::copy(K_d.begin(), K_d.end(), K_h.data());
+    }
+
+    void get_invQ(mat_ref invQ_h)
+    {
+        thrust::copy(invQ_d.begin(), invQ_d.end(), invQ_h.data());
+    }
+
+    void get_invQt(mat_ref invQt_h)
+    {
+        thrust::copy(invQt_d.begin(), invQt_d.end(), invQt_h.data());
+    }
+
+    double get_logpost(void)
+    {
+        double result;
+        CUBLASDOT(cublasHandle, n, dev_ptr(targets_d), 1, dev_ptr(invQt_d), 1,
+                  &result);
+
+        result += logdetC + n * log(2.0 * M_PI);
+
+        return 0.5 * result;
+    }
+
+    void logpost_deriv(vec_ref result)
+    {
+        double zero(0.0);
+        double one(1.0);
+        double half(0.5);
+        double m_half(-0.5);
+
+        const int Ntheta = D + 1;
+
+        // Compute
+        //   \pderiv{C_{jk}}{theta_i}
+        //
+        // The length of work_mat_d is n * testing_size
+        // The derivative above has  n * n * Ntheta components
+        // The following assumes that testing_size > Ntheta * n
+        cov_deriv_theta_batch_gpu(dev_ptr(work_mat_d),
+				  D, n, n,
+				  dev_ptr(inputs_d), dev_ptr(inputs_d),
+				  dev_ptr(theta_d));
+
+        // Compute
+        //   \deriv{logpost}{theta_i}
+        //     =   0.5 Tr(inv(C)*\pderiv{C}{\theta_i}           (Term 1)
+        //       - 0.5 (invQ_ts)^T \pderiv{C}{\theta_i} invQt  (Term 2)
+
+        // Working memory on device:
+        //   result_d:
+        //       accumulated result
+        //   work_mat_d (until Term 2 step 1):
+        //       covariance derivatives
+        //   work_mat_d (after Term 2 step 1):
+        //       temporary product in term 2 (okay to overwrite derivatives at this point)
+
+        // ***** Term 1 *****
+        // Treat 'jk' as a new 1-d index in C_{jk} and \pderiv{C_{jk}}{\theta}, called 'l' below.
+        // That is:
+        //    R_i = \pderiv{vec(C)_l}{\theta_i} vec(inv(C))_l
+        //
+        // Compute with gemv (y = \alpha A x + \beta y)
+        cublasDgemv(cublasHandle, CUBLAS_OP_T, // handle, op (transpose)
+                    n * n, Ntheta, // nrows, ncols (N.B. column-major order!)
+                    &half, dev_ptr(work_mat_d), n * n, // alpha, A, lda
+                    dev_ptr(invQ_d), 1, // x, incx
+                    &zero, // beta
+                    dev_ptr(result_d), 1); // y, incy
+
+        // ***** Term 2 *****
+        // First step:
+        //   S_ij = \pderiv{C_{jk}}{\theta_i} invQt_k
+        //
+        // Treat \pderiv{C_{jk}}{\theta_i} as (n * Ntheta) * n, treating 'ij' as a single index
+        // and take mat-vec product with invQt:
+        cublasDgemv(cublasHandle, CUBLAS_OP_T, // handle, op (transpose)
+                    n, n * Ntheta, // nrows, ncols (N.B. column-major order!)
+                    &one, dev_ptr(work_mat_d), n, // alpha, A, lda
+                    dev_ptr(invQt_d), 1, // x, incx
+                    &zero, // beta
+                    dev_ptr(work_mat_d), 1); // y, incy
+
+        // Second step:
+        //   R_i += 0.5 S_ij invQt_j
+        cublasDgemv(cublasHandle, CUBLAS_OP_T, // handle, op (transpose)
+                    n, Ntheta, // nrows, ncols (N.B. column-major order!)
+                    &m_half, dev_ptr(work_mat_d), n, // alpha, A, lda
+                    dev_ptr(invQt_d), 1, // x, incx
+                    &one, // beta
+                    dev_ptr(result_d), 1); // y, incy
+
+        thrust::copy(result_d.begin(), result_d.begin() + Ntheta, result.data());
+
+	// fitted nugget - last element of theta is log(nugget)
+	// partial deriv is 0.5*nugget*(trace(invQ) - invQt.invQt)
+	if (nug_type == NUG_FIT) {
+
+	  // trace of invQ
+	  REAL tr_invQ = 0.;
+	  thrust::device_vector<double> tr_invQ_d(1);
+	  trace(n, dev_ptr(invQ_d), dev_ptr(tr_invQ_d), dev_ptr(sum_buffer_d), sum_buffer_size_bytes);
+	  thrust::copy(tr_invQ_d.begin(), tr_invQ_d.end(), &tr_invQ);
+	  // invQt dot invQt
+	  REAL invQtSq = std::numeric_limits<double>::quiet_NaN();
+	  CUBLASDOT(cublasHandle, n, dev_ptr(invQt_d), 1, dev_ptr(invQt_d), 1,
+		    &invQtSq);
+
+	  // set the last element of result, putting it all together
+	  result(result.size()-1) = 0.5 * jitter * (tr_invQ - invQtSq);
+	}
+    }
+
+    // constructor
+    DenseGP_GPU(mat_ref inputs_, vec_ref targets_, unsigned int testing_size_)
+        : testing_size(testing_size_)
+        , n(inputs_.rows())
+        , D(inputs_.cols())
+        , K_d(n * n, 0.0)
+        , invQ_d(n * n, 0.0)
+        , chol_lower_d(n * n, 0.0)
+        , theta_d(D + 2)
+        , inputs(inputs_)
+        , targets(targets_)
+	, nug_type(NUG_ADAPTIVE)
+        , theta_fitted(false)
+        , inputs_d(inputs_.data(), inputs_.data() + D * n)
+        , targets_d(targets_.data(), targets_.data() + n)
+        , logdetC(0.0)
+        , jitter(0.0)
+        , invQt_d(n, 0.0)
+        , testing_d(D * testing_size, 0.0)
+        , work_d(n, 0.0)
+        , work_mat_d(n * testing_size, 0.0)
+        , sum_buffer_size_bytes(0)
+        , result_d(testing_size, 0.0)
+        , kappa_d(testing_size, 0.0)
+        , invQk_d(testing_size * n, 0.0)
+    {
+        cublasStatus_t cublas_status = cublasCreate(&cublasHandle);
+        if (cublas_status != CUBLAS_STATUS_SUCCESS)
+        {
+            throw std::runtime_error("CUBLAS initialization error\n");
+        }
+        cusolverStatus_t cusolver_status = cusolverDnCreate(&cusolverHandle);
+        if (cusolver_status != CUSOLVER_STATUS_SUCCESS)
+        {
+            throw std::runtime_error("cuSolver initialization error\n");
+        }
+
+        // CUBLAS potrf workspace
+        int potrfBufferSize;
+        cusolverDnDpotrf_bufferSize(cusolverHandle, CUBLAS_FILL_MODE_LOWER,
+                                    n, dev_ptr(potrf_buffer_d), n,
+                                    &potrfBufferSize);
+        potrf_buffer_d.resize(potrfBufferSize);
+
+        // The following call determines the size of the cub::DeviceReduce::Sum workspace:
+        // sum_buffer_size_bytes is 0 before this call, and the size of sum_buffer_d afterwards.
+        // The end iterators are supplied but are not used.
+        cub::DeviceReduce::Sum(dev_ptr(sum_buffer_d), sum_buffer_size_bytes,
+                               sum_buffer_d.end(), sum_buffer_d.end(), n);
+        sum_buffer_d.resize(sum_buffer_size_bytes);
+    }
+};
+
+#endif
diff --git a/mogp_gpu/src/strided_range.hpp b/mogp_gpu/src/strided_range.hpp
new file mode 100644
index 00000000..626ce4fb
--- /dev/null
+++ b/mogp_gpu/src/strided_range.hpp
@@ -0,0 +1,88 @@
+#ifndef STRIDED_RANGE_HPP
+#define STRIDED_RANGE_HPP
+
+#include <thrust/iterator/counting_iterator.h>
+#include <thrust/iterator/transform_iterator.h>
+#include <thrust/iterator/permutation_iterator.h>
+#include <thrust/functional.h>
+#include <thrust/device_vector.h>
+#include <iostream>
+
+// ----------------------------------------
+// The `strided_range` class, from a portion of the file found at:
+//
+//      https://github.com/thrust/thrust/blob/master/examples/strided_range.cu
+//
+//  Copyright NVIDIA Corporation
+//
+//  Licensed under the Apache License, Version 2.0 (the "License")
+//
+//  You may obtain a copy of the License at
+//
+//      http://www.apache.org/licenses/LICENSE-2.0
+//
+//  Unless required by applicable law or agreed to in writing, software
+//  distributed under the License is distributed on an "AS IS" BASIS,
+//  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+//  See the License for the specific language governing permissions and
+//  limitations under the License.
+// 
+template <typename Iterator>
+class strided_range
+{
+    public:
+
+    typedef typename thrust::iterator_difference<Iterator>::type difference_type;
+
+    struct stride_functor : public thrust::unary_function<difference_type,difference_type>
+    {
+        difference_type stride;
+
+        stride_functor(difference_type stride)
+            : stride(stride) {}
+
+        __host__ __device__
+        difference_type operator()(const difference_type& i) const
+        { 
+            return stride * i;
+        }
+    };
+
+    typedef typename thrust::counting_iterator<difference_type>                   CountingIterator;
+    typedef typename thrust::transform_iterator<stride_functor, CountingIterator> TransformIterator;
+    typedef typename thrust::permutation_iterator<Iterator,TransformIterator>     PermutationIterator;
+
+    // type of the strided_range iterator
+    typedef PermutationIterator iterator;
+
+    // construct strided_range for the range [first,last)
+    strided_range(Iterator first, Iterator last, difference_type stride)
+        : first(first), last(last), stride(stride) {}
+   
+    iterator begin(void) const
+    {
+        return PermutationIterator(first, TransformIterator(CountingIterator(0), stride_functor(stride)));
+    }
+
+    iterator end(void) const
+    {
+        return begin() + ((last - first) + (stride - 1)) / stride;
+    }
+    
+    protected:
+    Iterator first;
+    Iterator last;
+    difference_type stride;
+};
+
+// ----------------------------------------
+
+template <typename Iterator>
+strided_range<Iterator> make_strided_range(
+    Iterator first, Iterator last,
+    typename thrust::iterator_difference<Iterator>::type stride)
+{
+    return strided_range<Iterator>(first, last, stride);
+}
+
+#endif // STRIDED_RANGE_HPP
diff --git a/mogp_gpu/src/util.cu b/mogp_gpu/src/util.cu
new file mode 100644
index 00000000..8ed513d1
--- /dev/null
+++ b/mogp_gpu/src/util.cu
@@ -0,0 +1,29 @@
+__global__ void identity_kernel(int N, double *A)
+{
+    int i = blockDim.x*blockIdx.x + threadIdx.x;
+    int j = blockDim.y*blockIdx.y + threadIdx.y;
+    if (i < N && j < N) A[i + N * j] = (i == j)?1.0:0.0;
+}
+
+void identity_device(int N, double *A)
+{
+    const int Nx = 8, Ny = 8;
+    dim3 threads_per_block(Nx, Ny);
+    dim3 blocks((N + Nx - 1)/Nx, (N + Ny - 1)/Ny);
+    identity_kernel<<<blocks, threads_per_block>>>(N, A);
+}
+
+__global__ void add_diagonal_kernel(int N, double b, double *A)
+{
+    int i = blockDim.x*blockIdx.x + threadIdx.x;
+    int j = blockDim.y*blockIdx.y + threadIdx.y;
+    if (i < N && j < N && i == j) A[i + N * j] += b;
+}
+
+void add_diagonal(int N, double b, double *A)
+{
+    const int Nx = 8, Ny = 8;
+    dim3 threads_per_block(Nx, Ny);
+    dim3 blocks((N + Nx - 1)/Nx, (N + Ny - 1)/Ny);
+    add_diagonal_kernel<<<blocks, threads_per_block>>>(N, b, A);
+}
diff --git a/mogp_gpu/src/util.hpp b/mogp_gpu/src/util.hpp
new file mode 100644
index 00000000..545f308e
--- /dev/null
+++ b/mogp_gpu/src/util.hpp
@@ -0,0 +1,7 @@
+#ifndef UTIL_HPP
+#define UTIL_HPP
+
+void identity_device(int N, double *A);
+void add_diagonal(int N, double b, double *A);
+
+#endif
diff --git a/setup.py b/setup.py
index 9b77e53a..55ceb83e 100644
--- a/setup.py
+++ b/setup.py
@@ -1,5 +1,11 @@
+import os
 import setuptools
 
+from distutils.extension import Extension
+from distutils.command.build_ext import build_ext
+
+import numpy as np
+
 # version information
 MAJOR = 0
 MINOR = 5
@@ -26,9 +32,99 @@ def write_version_file(version):
 
 write_version_file(version)
 
+
 with open("README.md", "r") as fh:
     long_description = fh.read()
 
+
+def get_cuda_config():
+    """
+    Find the nvcc executable in the path, and
+    return a dictionary with keys 'home', 'nvcc', 'include', 'lib64'.
+    """
+    path = os.environ["PATH"]
+    exe = "nvcc"
+    binpath = None
+    for location in path.split(os.pathsep):
+        if os.path.exists(os.path.join(location, exe)):
+            binpath = os.path.abspath(os.path.join(location, exe))
+            break
+    if not binpath:
+        print("Unable to find 'nvcc' in path")
+        return {}
+    home = os.path.dirname(os.path.dirname(binpath))
+    include = os.path.join(home, "include")
+    lib64 = os.path.join(home, "lib64")
+    return {
+        "home": home,
+        "nvcc": binpath,
+        "include": include,
+        "lib64": lib64
+    }
+
+
+def customize_compiler_for_nvcc(self):
+    """
+    Taken from SO
+    https://stackoverflow.com/questions/10034325/can-python-distutils-compile-cuda-code
+    by Robert T. McGibbon.
+    Modify the _compile of distutils extension builder.
+    """
+    # add .cu to list of extensions compile knows about
+    self.src_extensions += ['.cu',".hpp"]
+    # save references to original methods
+    default_compiler_so = self.compiler_so
+    super = self._compile
+    # redefine _compile method.
+    def _compile(obj, src, ext, cc_args, extra_postargs, pp_opts):
+        if os.path.splitext(src)[1] == '.cu':
+            # use nvcc for .cu files
+            cuda_config = get_cuda_config()
+            if not "nvcc" in cuda_config.keys():
+                return
+            self.set_executable('compiler_so', cuda_config['nvcc'])
+            # use subset of extra_postargs
+            postargs = extra_postargs['nvcc']
+        else:
+            postargs = extra_postargs['gcc']
+
+        super(obj, src, ext, cc_args, postargs, pp_opts)
+        # reset default compiler_so
+        self.compiler_so = default_compiler_so
+    # replace the class default _compile method with this one
+    self._compile = _compile
+
+# Run the custom compiler
+class custom_build_ext(build_ext):
+    def build_extensions(self):
+        customize_compiler_for_nvcc(self.compiler)
+        build_ext.build_extensions(self)
+
+cuda_config = get_cuda_config()
+
+# we only want to add the mogp_gpu extension if we have cuda compiler
+ext_modules = []
+if len(cuda_config) > 0:
+    import pybind11
+    pybind_include = pybind11.get_include()
+    numpy_include = np.get_include()
+    ext = Extension("libgpgpu",
+                    sources=["mogp_gpu/src/gp_gpu.cu",
+                             "mogp_gpu/src/cov_gpu.cu",
+                             "mogp_gpu/src/util.cu"],
+                    library_dirs=[cuda_config["lib64"]],
+                    libraries=["cudart","cublas","cusolver"],
+                    runtime_library_dirs=[cuda_config["lib64"]],
+                    extra_compile_args={"gcc":["-std=c++14"],
+                                    "nvcc": ["--compiler-options",
+                                             "-O3,-Wall,-shared,-std=c++14,-fPIC",
+                                             "-arch=sm_60",
+                                             "--generate-code","arch=compute_37,code=sm_37",
+                                             "--generate-code","arch=compute_60,code=sm_60"
+                                    ]},
+                    include_dirs=[numpy_include, pybind_include, cuda_config["include"],"mogp_gpu/src"])
+    ext_modules.append(ext)
+
 setuptools.setup(name='mogp_emulator',
       version=version,
       description='Tool for Multi-Output Gaussian Process Emulators',
@@ -44,4 +140,7 @@ def write_version_file(version):
       author_email='edaub@turing.ac.uk',
       packages=setuptools.find_packages(),
       license=['MIT'],
-      install_requires=['numpy', 'scipy'])
+      ext_modules=ext_modules,
+      cmdclass={"build_ext": custom_build_ext},
+      install_requires=['numpy', 'scipy'],
+      zip_safe=False)


Test patch:
diff --git a/mogp_emulator/tests/test_GaussianProcess.py b/mogp_emulator/tests/test_GaussianProcess.py
index 09870c64..5ac0a47d 100644
--- a/mogp_emulator/tests/test_GaussianProcess.py
+++ b/mogp_emulator/tests/test_GaussianProcess.py
@@ -2,12 +2,16 @@
 import numpy as np
 import pytest
 from numpy.testing import assert_allclose
+from ..LibGPGPU import gpu_usable
 from ..GaussianProcess import GaussianProcess, PredictResult
+from ..GaussianProcessGPU import GaussianProcessGPU
 from ..MeanFunction import ConstantMean, LinearMean, MeanFunction
 from ..Kernel import SquaredExponential, Matern52
 from ..Priors import NormalPrior, GammaPrior, InvGammaPrior
 from scipy import linalg
 
+GPU_NOT_FOUND_MSG = "A compatible GPU could not be found or the GPU library (libgpgpu) could not be loaded"
+
 @pytest.fixture
 def x():
     return np.array([[1., 2., 3.], [4., 5., 6.]])
@@ -48,6 +52,28 @@ def test_GaussianProcess_init(x, y):
 
     assert str(gp.mean) == "c + c*x[0]"
 
+
+@pytest.mark.skipif(not gpu_usable(), reason=GPU_NOT_FOUND_MSG)
+def test_GaussianProcessGPU_init(x, y):
+    "Test function for correct functioning of the init method of GaussianProcess"
+
+    gp = GaussianProcessGPU(x, y)
+    assert_allclose(x, gp.inputs)
+    assert_allclose(y, gp.targets)
+    assert gp.D == 3
+    assert gp.n == 2
+    assert gp.nugget == None
+    assert gp.nugget_type == "adaptive"
+    gp = GaussianProcessGPU(y, y)
+    assert gp.inputs.shape == (2, 1)
+
+    gp = GaussianProcessGPU(x, y, nugget=1.e-12)
+    assert_allclose(gp.nugget, 1.e-12)
+
+    gp = GaussianProcessGPU(x, y, kernel="SquaredExponential")
+    assert isinstance(gp.kernel, SquaredExponential)
+
+
 def test_GP_init_failures(x, y):
     "Tests that GaussianProcess fails correctly with bad inputs"
 
@@ -72,6 +98,33 @@ def test_GP_init_failures(x, y):
     with pytest.raises(ValueError):
         gp = GaussianProcess(x, y, nugget="a")
 
+
+@pytest.mark.skipif(not gpu_usable(), reason=GPU_NOT_FOUND_MSG)
+def test_GPGPU_init_failures(x, y):
+    "Tests that GaussianProcessGPU fails correctly with bad inputs"
+
+    with pytest.raises(AssertionError):
+        gp = GaussianProcessGPU(np.ones((2, 2, 2)), y)
+
+    with pytest.raises(AssertionError):
+        gp = GaussianProcessGPU(x, x)
+
+    with pytest.raises(AssertionError):
+        gp = GaussianProcessGPU(np.ones((2, 3)), np.ones(3))
+
+    with pytest.raises(ValueError):
+        gp = GaussianProcessGPU(x, y, mean=1)
+
+    with pytest.raises(ValueError):
+        gp = GaussianProcessGPU(x, y, kernel="blah")
+
+    with pytest.raises(ValueError):
+        gp = GaussianProcessGPU(x, y, kernel=1)
+
+    with pytest.raises(ValueError):
+        gp = GaussianProcessGPU(x, y, nugget="a")
+
+
 def test_GaussianProcess_n_params(x, y):
     "test the get_n_params method of GaussianProcess"
 
@@ -81,6 +134,15 @@ def test_GaussianProcess_n_params(x, y):
     gp = GaussianProcess(x, y, mean="x[0]")
     assert gp.n_params == 2 + x.shape[1] + 2
 
+
+@pytest.mark.skipif(not gpu_usable(), reason=GPU_NOT_FOUND_MSG)
+def test_GaussianProcessGPU_n_params(x, y):
+    "test the get_n_params method of GaussianProcessGPU"
+
+    gp = GaussianProcessGPU(x, y)
+    assert gp.n_params == x.shape[1] + 2
+
+
 def test_GaussianProcess_nugget(x, y):
     "Tests the get_nugget method of GaussianProcess"
 
@@ -103,7 +165,37 @@ def test_GaussianProcess_nugget(x, y):
     gp.nugget = "pivot"
     assert gp.nugget is None
     assert gp.nugget_type == "pivot"
-    
+
+    with pytest.raises(TypeError):
+        gp.nugget = [1]
+
+    with pytest.raises(ValueError):
+        gp.nugget = "blah"
+
+    with pytest.raises(ValueError):
+        gp.nugget = -1.
+
+
+@pytest.mark.skipif(not gpu_usable(), reason=GPU_NOT_FOUND_MSG)
+def test_GaussianProcessGPU_nugget(x, y):
+    "Tests the get_nugget method of GaussianProcessGPU"
+
+    gp = GaussianProcessGPU(x, y)
+    assert gp.nugget is None
+    assert gp.nugget_type == "adaptive"
+
+    gp.nugget = "fit"
+    assert gp.nugget is None
+    assert gp.nugget_type == "fit"
+
+    gp.nugget = 1.
+    assert_allclose(gp.nugget, 1.)
+    assert gp.nugget_type == "fixed"
+
+    gp.nugget = 0
+    assert_allclose(gp.nugget, 0.)
+    assert gp.nugget_type == "fixed"
+
     with pytest.raises(TypeError):
         gp.nugget = [1]
 
@@ -150,6 +242,39 @@ def test_GaussianProcess_theta(x, y, mean, nugget, sn):
     assert_allclose(invQt_expect, gp.invQt)
     assert_allclose(logpost_expect, gp.current_logpost)
 
+@pytest.mark.skipif(not gpu_usable(), reason=GPU_NOT_FOUND_MSG)
+@pytest.mark.parametrize("nugget,sn", [(0., 1.), ("adaptive", 0.)]) # ("fit", np.log(1.e-6))])
+def test_GaussianProcessGPU_theta(x, y, nugget, sn):
+    "test the theta property of GaussianProcess (effectively the same as fit)"
+
+    # zero mean, zero nugget
+
+    gp = GaussianProcessGPU(x, y, nugget=nugget)
+
+    theta = np.ones(gp.n_params)
+    theta[-1] = sn
+
+    gp.theta = theta
+
+    if nugget == "adaptive" or nugget == 0.:
+        assert gp.nugget == 0.
+        noise = 0.
+    else:
+        assert_allclose(gp.nugget, np.exp(sn))
+        noise = np.exp(sn)*np.eye(x.shape[0])
+    Q = gp.kernel.kernel_f(x, x, theta[:-1]) + noise
+
+    L_expect = np.linalg.cholesky(Q)
+    invQt_expect = np.linalg.solve(Q, y)
+    logpost_expect = 0.5*(np.log(np.linalg.det(Q)) +
+                          np.dot(y, invQt_expect) +
+                          gp.n*np.log(2.*np.pi))
+
+    assert_allclose(L_expect, gp.L)
+    assert_allclose(invQt_expect, gp.invQt)
+    assert_allclose(logpost_expect, gp.current_logpost)
+
+
 def test_GaussianProcess_theta_pivot():
     "test that pivoting works as expected"
 
@@ -169,7 +294,7 @@ def test_GaussianProcess_theta_pivot():
     assert_allclose(gp1.L, gp2.L)
     assert_allclose(gp1.invQt, gp2.invQt[gp2.P])
     assert np.array_equal(gp2.P, [0, 2, 1])
-    
+
 def test_GaussianProcess_priors(x, y):
     "test that priors are set properly"
 
@@ -253,6 +378,7 @@ def test_GaussianProcess_fit_logposterior(x, y, mean, nugget, sn):
     assert_allclose(logpost_expect, gp.current_logpost)
     assert_allclose(logpost_expect, gp.logposterior(theta))
 
+
 def test_GaussianProcess_logposterior(x, y):
     "test logposterior method of GaussianProcess"
 
@@ -278,6 +404,32 @@ def test_GaussianProcess_logposterior(x, y):
     assert_allclose(invQt_expect, gp.invQt)
     assert_allclose(logpost_expect, gp.current_logpost)
 
+@pytest.mark.skipif(not gpu_usable(), reason=GPU_NOT_FOUND_MSG)
+def test_GaussianProcessGPU_logposterior(x, y):
+    "test logposterior method of GaussianProcessGPU"
+
+    # logposterior already tested, but check that parameters are re-fit if changed
+
+    gp = GaussianProcessGPU(x, y, nugget = 0.)
+
+    theta = np.ones(gp.n_params)
+    gp.fit(theta)
+
+    theta = np.zeros(gp.n_params)
+
+    Q = gp.kernel.kernel_f(x, x, theta[:-1])
+
+    L_expect = np.linalg.cholesky(Q)
+    invQt_expect = np.linalg.solve(Q, y)
+    logpost_expect = 0.5*(np.log(np.linalg.det(Q)) +
+                          np.dot(y, invQt_expect) +
+                          gp.n*np.log(2.*np.pi))
+
+    assert_allclose(logpost_expect, gp.logposterior(theta))
+    assert_allclose(gp.L, L_expect)
+    assert_allclose(invQt_expect, gp.invQt)
+    assert_allclose(logpost_expect, gp.current_logpost)
+
 @pytest.fixture
 def dx():
     return 1.e-6
@@ -303,8 +455,31 @@ def test_GaussianProcess_logpost_deriv(x, y, dx, mean, nugget, sn):
 
     assert_allclose(deriv, gp.logpost_deriv(theta), atol=1.e-7, rtol=1.e-5)
 
+@pytest.mark.skipif(not gpu_usable(), reason=GPU_NOT_FOUND_MSG)
+@pytest.mark.parametrize("nugget,sn", [(0., 1.), ("adaptive", 1.),
+                                            ("fit", np.log(1.e-6))])
+def test_GaussianProcessGPU_logpost_deriv(x, y, dx, nugget, sn):
+    "test logposterior derivatives for GaussianProcessGPU via finite differences"
+
+    gp = GaussianProcessGPU(x, y, nugget=nugget)
+
+    n = gp.n_params
+    theta = np.ones(n)
+    theta[-1] = sn
+
+    deriv = np.zeros(n)
+
+    for i in range(n):
+        dx_array = np.zeros(n)
+        dx_array[i] = dx
+        deriv[i] = (gp.logposterior(theta) - gp.logposterior(theta - dx_array))/dx
+
+    assert_allclose(deriv, gp.logpost_deriv(theta), atol=1.e-7, rtol=1.e-5)
+
 @pytest.mark.parametrize("mean,nugget,sn", [(None, 0., 1.),
                                             (None, "pivot", 1.), ("x[0]", "fit", np.log(1.e-6))])
+
+
 def test_GaussianProcess_logpost_hessian(x, y, dx, mean, nugget, sn):
     "test the hessian method of GaussianProcess with finite differences"
 
@@ -498,6 +673,95 @@ def test_GaussianProcess_predict(x, y, dx):
 
     assert_allclose(gp(x_test), mu_expect)
 
+@pytest.mark.skipif(not gpu_usable(), reason=GPU_NOT_FOUND_MSG)
+def test_GaussianProcessGPU_predict(x, y, dx):
+    "test the predict method of GaussianProcessGPU"
+
+    # zero mean
+
+    gp = GaussianProcessGPU(x, y, nugget=0.)
+    theta = np.ones(gp.n_params)
+
+    gp.fit(theta)
+
+    x_test = np.array([[2., 3., 4.]])
+
+    mu, var, deriv = gp.predict(x_test)
+
+    K = gp.kernel.kernel_f(x, x, theta[:-1])
+    Ktest = gp.kernel.kernel_f(x_test, x, theta[:-1])
+
+    mu_expect = np.dot(Ktest, gp.invQt)
+    var_expect = np.exp(theta[-2]) - np.diag(np.dot(Ktest, np.linalg.solve(K, Ktest.T)))
+
+    D = gp.D
+
+    deriv_expect = np.zeros((1, D))
+
+    for i in range(D):
+        dx_array = np.zeros(D)
+        dx_array[i] = dx
+        deriv_expect[0, i] = (gp.predict(x_test)[0] - gp.predict(x_test - dx_array)[0])/dx
+
+    assert_allclose(mu, mu_expect)
+    assert_allclose(var, var_expect)
+    assert_allclose(deriv, deriv_expect, atol=1.e-7, rtol=1.e-7)
+
+    # check that reshaping works as expected
+
+    x_test = np.array([2., 3., 4.])
+
+    mu, var, deriv = gp.predict(x_test)
+
+    assert_allclose(mu, mu_expect)
+    assert_allclose(var, var_expect)
+    assert_allclose(deriv, deriv_expect, atol=1.e-7, rtol=1.e-7)
+
+    # check that with 1D input data can handle 1D prediction data correctly
+
+    gp = GaussianProcessGPU(y, y, nugget=0.)
+
+    gp.fit(np.ones(gp.n_params))
+
+    n_predict = 51
+    mu, var, deriv = gp.predict(np.linspace(0., 1., n_predict))
+
+    assert mu.shape == (n_predict,)
+    assert var.shape == (n_predict,)
+    assert deriv.shape == (n_predict, 1)
+
+    # check unc and deriv flags work
+
+    _, var, deriv = gp.predict(x_test, unc=False, deriv=False)
+
+    assert var is None
+    assert deriv is None
+
+    # check that the returned PredictResult works correctly
+    gp = GaussianProcessGPU(x, y, nugget=0.)
+    theta = np.ones(gp.n_params)
+    gp.fit(theta)
+    x_test = np.array([[2., 3., 4.]])
+
+    pr = gp.predict(x_test)
+
+    assert_allclose(pr.mean, mu_expect)
+    assert_allclose(pr.unc, var_expect)
+    assert_allclose(pr.deriv, deriv_expect, atol=1.e-7, rtol=1.e-7)
+
+    assert_allclose(pr['mean'], mu_expect)
+    assert_allclose(pr['unc'], var_expect)
+    assert_allclose(pr['deriv'], deriv_expect, atol=1.e-7, rtol=1.e-7)
+
+    assert_allclose(pr[0], mu_expect)
+    assert_allclose(pr[1], var_expect)
+    assert_allclose(pr[2], deriv_expect, atol=1.e-7, rtol=1.e-7)
+
+    # check that calling gp is equivalent to predicting
+
+    assert_allclose(gp(x_test), mu_expect)
+
+
 def test_GaussianProcess_predict_nugget(x, y):
     "test that the nugget works correctly when making predictions"
 
@@ -522,11 +786,38 @@ def test_GaussianProcess_predict_nugget(x, y):
 
     assert_allclose(preds.unc, var_expect, atol=1.e-7)
 
+
+@pytest.mark.skipif(not gpu_usable(), reason=GPU_NOT_FOUND_MSG)
+def test_GaussianProcessGPU_predict_nugget(x, y):
+    "test that the nugget works correctly when making predictions"
+
+    nugget = 1.e0
+
+    gp = GaussianProcessGPU(x, y, nugget=nugget)
+    theta = np.ones(gp.n_params)
+
+    gp.fit(theta)
+
+    preds = gp.predict(x)
+
+    K = gp.kernel.kernel_f(x, x, theta[:-1])
+
+    var_expect = np.exp(theta[-2]) + nugget - np.diag(np.dot(K, np.linalg.solve(K + np.eye(gp.n)*nugget, K)))
+
+    assert_allclose(preds.unc, var_expect, atol=1.e-7)
+
+    preds = gp.predict(x, include_nugget=False)
+
+    var_expect = np.exp(theta[-2]) - np.diag(np.dot(K, np.linalg.solve(K + np.eye(gp.n)*nugget, K)))
+
+    assert_allclose(preds.unc, var_expect, atol=1.e-7)
+
+
 def test_GaussianProcess_predict_pivot():
     "test that pivoting gives same predictions as standard version"
 
     # input arrays are re-ordered such that pivoting should re-order the second to match the first
-    
+
     x1 = np.array([1., 4., 2.])
     x2 = np.array([1., 2., 4.])
     y1 = np.array([1., 1., 2.])
@@ -546,7 +837,8 @@ def test_GaussianProcess_predict_pivot():
     assert_allclose(mean1, mean2)
     assert_allclose(var1, var2)
     assert_allclose(deriv1, deriv2)
-    
+
+
 def test_GaussianProcess_predict_variance():
     "confirm that caching factorized matrix produces stable variance predictions"
 
@@ -566,6 +858,7 @@ def test_GaussianProcess_predict_variance():
 
     assert_allclose(np.zeros(101), var, atol = 1.e-3)
 
+
 def test_GaussianProcess_predict_failures(x, y):
     "test situations where predict method of GaussianProcess should fail"
 
@@ -583,8 +876,36 @@ def test_GaussianProcess_predict_failures(x, y):
     with pytest.raises(AssertionError):
         gp.predict(np.array([[2., 4.]]))
 
+
+@pytest.mark.skipif(not gpu_usable(), reason=GPU_NOT_FOUND_MSG)
+def test_GaussianProcessGPU_predict_failures(x, y):
+    "test situations where predict method of GaussianProcessGPU should fail"
+
+    gp = GaussianProcessGPU(x, y)
+
+    with pytest.raises(ValueError):
+        gp.predict(np.array([2., 3., 4.]))
+
+    theta = np.ones(gp.n_params)
+    gp.fit(theta)
+
+    with pytest.raises(AssertionError):
+        gp.predict(np.ones((2, 2, 2)))
+
+    with pytest.raises(AssertionError):
+        gp.predict(np.array([[2., 4.]]))
+
+
 def test_GaussianProcess_str(x, y):
     "Test function for string method"
 
     gp = GaussianProcess(x, y)
     assert (str(gp) == "Gaussian Process with {} training examples and {} input variables".format(x.shape[0], x.shape[1]))
+
+
+@pytest.mark.skipif(not gpu_usable(), reason=GPU_NOT_FOUND_MSG)
+def test_GaussianProcessGPU_str(x, y):
+    "Test function for string method"
+
+    gp = GaussianProcessGPU(x, y)
+    assert (str(gp) == "Gaussian Process with {} training examples and {} input variables".format(x.shape[0], x.shape[1]))
diff --git a/mogp_emulator/tests/test_HistoryMatching.py b/mogp_emulator/tests/test_HistoryMatching.py
index 7b4dd062..877e5759 100644
--- a/mogp_emulator/tests/test_HistoryMatching.py
+++ b/mogp_emulator/tests/test_HistoryMatching.py
@@ -4,6 +4,10 @@
 from ..HistoryMatching import HistoryMatching
 from ..fitting import fit_GP_MAP
 from ..GaussianProcess import GaussianProcess, PredictResult
+from ..LibGPGPU import gpu_usable
+from ..GaussianProcessGPU import GaussianProcessGPU
+
+GPU_NOT_FOUND_MSG = "A compatible GPU could not be found or the GPU library (libgpgpu) could not be loaded"
 
 def get_y_simulated_1D(x):
     n_points = len(x)
@@ -131,6 +135,129 @@ def test_sanity_checks():
     hm = HistoryMatching(obs=obs, gp=gp, coords=coords)
     I = hm.get_implausibility(7.)
 
+
+@pytest.mark.skipif(not gpu_usable(), reason=GPU_NOT_FOUND_MSG)
+def test_sanity_checks_GPU():
+    "test basic functioning of HistoryMatching"
+
+    # Create a gaussian process
+    x_training = np.array([
+        [0.],
+        [10.],
+        [20.],
+        [30.],
+        [43.],
+        [50.]
+    ])
+
+    y_training = get_y_simulated_1D(x_training)
+
+    np.random.seed(47)
+    gp = GaussianProcessGPU(x_training, y_training)
+    gp = fit_GP_MAP(gp)
+
+    # Define observation and implausibility threshold
+    obs = [-0.8, 0.0004]
+
+    # Coords to predict
+    n_rand = 2000
+    x_predict_min = -3
+    x_predict_max = 53
+    x_predict = np.random.rand(n_rand)
+    x_predict = np.sort(x_predict, axis=0)
+    x_predict *= (x_predict_max - x_predict_min)
+    x_predict += x_predict_min
+    x_predict = x_predict[:, None]
+
+    coords = x_predict
+
+    expectations = gp.predict(coords)
+
+    # Create History Matching Instance
+    print("---TEST INPUTS---")
+    print("No Args")
+    hm = HistoryMatching()
+    hm.status()
+
+    print("Obs Only a - list")
+    hm = HistoryMatching(obs=obs)
+    hm.status()
+
+    print("Obs only b - single-element list")
+    hm = HistoryMatching(obs=[3.])
+    hm.status()
+
+    print("Obs only c - single-value")
+    hm = HistoryMatching(obs=3.)
+    hm.status()
+
+    print("gp Only")
+    hm = HistoryMatching(gp=gp)
+    hm.status()
+
+    print("Coords only a - 2d ndarray")
+    hm = HistoryMatching(coords=coords)
+    hm.status()
+
+    print("Coords only b - 1d ndarray")
+    hm = HistoryMatching(coords=np.random.rand(n_rand))
+    hm.status()
+
+    print("Coords only c - list")
+    hm = HistoryMatching(coords=[a for a in range(n_rand)])
+    hm.status()
+
+    print("Expectation only")
+    hm = HistoryMatching(expectations=expectations)
+    hm.status()
+
+    print("Threshold Only")
+    hm = HistoryMatching(threshold=3.)
+    hm.status()
+
+
+    print("---TEST ASSIGNMENT---")
+    print("Assign gp")
+    hm = HistoryMatching(obs)
+    hm.status()
+    hm.set_gp(gp)
+    hm.status()
+
+    print("Assign Obs")
+    hm = HistoryMatching(gp)
+    hm.status()
+    hm.set_obs(obs)
+    hm.status()
+
+    print("Assign Coords")
+    hm = HistoryMatching()
+    hm.status()
+    hm.set_coords(coords)
+    hm.status()
+
+    print("Assign Expectations")
+    hm = HistoryMatching()
+    hm.status()
+    hm.set_expectations(expectations)
+    hm.status()
+
+    print("Assign Threshold")
+    hm = HistoryMatching()
+    hm.status()
+    hm.set_threshold(3.)
+    hm.status()
+
+
+    print("---TEST IMPLAUSABILIY---")
+    print("implausibility test a - no vars")
+    hm = HistoryMatching(obs=obs, gp=gp, coords=coords)
+    I = hm.get_implausibility()
+
+    print("implausibility test b - single value")
+    hm = HistoryMatching(obs=obs, gp=gp, coords=coords)
+    I = hm.get_implausibility(7.)
+
+
 def test_HistoryMatching_init():
     "test the init method of HistoryMatching"
 
@@ -335,6 +462,25 @@ def test_HistoryMatching_set_gp():
     with pytest.raises(TypeError):
         hm.set_gp(gp=1.)
 
+@pytest.mark.skipif(not gpu_usable(), reason=GPU_NOT_FOUND_MSG)
+def test_HistoryMatching_set_gpGPU():
+    'test the set_gp method of HistoryMatching'
+
+    # correct functionality
+
+    gp = GaussianProcessGPU(np.reshape(np.linspace(0., 1.), (-1, 1)), np.linspace(0., 1.))
+
+    hm = HistoryMatching()
+    hm.set_gp(gp)
+
+    assert hm.gp == gp
+
+    # bad type for GP
+
+    with pytest.raises(TypeError):
+        hm.set_gp(gp=1.)
+
+
 def test_HistoryMatching_set_obs():
     'test the set_obs method of HistoryMatching'
 
diff --git a/mogp_emulator/tests/test_MultiOutputGP.py b/mogp_emulator/tests/test_MultiOutputGP.py
index 21c3e910..0c54d995 100644
--- a/mogp_emulator/tests/test_MultiOutputGP.py
+++ b/mogp_emulator/tests/test_MultiOutputGP.py
@@ -3,10 +3,13 @@
 from numpy.testing import assert_allclose
 from ..GaussianProcess import GaussianProcess, PredictResult
 from ..MultiOutputGP import MultiOutputGP
+from ..LibGPGPU import gpu_usable
 from ..MeanFunction import ConstantMean, LinearMean, MeanFunction
 from ..Kernel import Matern52
 from scipy import linalg
 
+GPU_NOT_FOUND_MSG = "A compatible GPU could not be found or the GPU library (libgpgpu) could not be loaded"
+
 @pytest.fixture
 def x():
     return np.array([[1., 2., 3.], [4., 5., 6.]])
@@ -27,6 +30,20 @@ def test_MultiOutputGP_init(x, y):
     assert gp.n == 2
     assert gp.n_emulators == 2
 
+
+@pytest.mark.skipif(not gpu_usable(), reason=GPU_NOT_FOUND_MSG)
+def test_MultiOutputGP_GPU_init(x, y):
+    "Test function for correct functioning of the init method of GaussianProcessGPU"
+
+    gp = MultiOutputGP(x, y, use_gpu=True)
+    assert_allclose(x, gp.emulators[0].inputs)
+    assert_allclose(x, gp.emulators[1].inputs)
+    assert_allclose(y[0], gp.emulators[0].targets)
+    assert_allclose(y[1], gp.emulators[1].targets)
+    assert gp.D == 3
+    assert gp.n == 2
+    assert gp.n_emulators == 2
+
 @pytest.fixture
 def dx():
     return 1.e-6
@@ -97,7 +114,7 @@ def test_MultiOutputGP_predict(x, y, dx):
         assert_allclose(var[i], var_expect)
 
     # test behavior if not all emulators are fit
-    
+
     gp.emulators[1].theta = None
 
     mu, var, deriv = gp.predict(x_test, allow_not_fit=True)
@@ -115,6 +132,88 @@ def test_MultiOutputGP_predict(x, y, dx):
         gp.predict(x_test)
 
 
+@pytest.mark.skipif(not gpu_usable(), reason=GPU_NOT_FOUND_MSG)
+def test_MultiOutputGP_GPU_predict(x, y, dx):
+    "test the predict method of GaussianProcess"
+
+    gp = MultiOutputGP(x, y, nugget=0., use_gpu=True)
+    theta = np.ones(gp.emulators[0].n_params)
+
+    gp.emulators[0].fit(theta)
+    gp.emulators[1].fit(theta)
+
+    x_test = np.array([[2., 3., 4.]])
+
+    mu, var, deriv = gp.predict(x_test)
+
+    for i in range(2):
+
+        K = gp.emulators[i].kernel.kernel_f(x, x, theta[:-1])
+        Ktest = gp.emulators[i].kernel.kernel_f(x_test, x, theta[:-1])
+
+        mu_expect = np.dot(Ktest, gp.emulators[i].invQt)
+        var_expect = np.exp(theta[-2]) - np.diag(np.dot(Ktest, np.linalg.solve(K, Ktest.T)))
+
+        D = gp.D
+
+        deriv_expect = np.zeros((1, D))
+
+        for j in range(D):
+            dx_array = np.zeros(D)
+            dx_array[j] = dx
+            deriv_expect[0, j] = (gp.emulators[i].predict(x_test)[0] - gp.emulators[i].predict(x_test - dx_array)[0])/dx
+
+        assert_allclose(mu[i], mu_expect)
+        assert_allclose(var[i], var_expect)
+        assert_allclose(deriv[i], deriv_expect, atol=1.e-6, rtol=1.e-6)
+
+    nugget = 1.
+    gp = MultiOutputGP(x, y, nugget=nugget, use_gpu=True)
+    theta = np.ones(gp.emulators[0].n_params)
+
+    gp.emulators[0].fit(theta)
+    gp.emulators[1].fit(theta)
+
+    x_test = np.array([[2., 3., 4.]])
+
+    mu, var, deriv = gp.predict(x_test)
+
+    for i in range(2):
+
+        K = gp.emulators[i].kernel.kernel_f(x, x, theta[:-1]) + np.eye(gp.emulators[i].n)*nugget
+        Ktest = gp.emulators[i].kernel.kernel_f(x_test, x, theta[:-1])
+
+        var_expect = np.exp(theta[-2]) + nugget - np.diag(np.dot(Ktest, np.linalg.solve(K, Ktest.T)))
+
+        assert_allclose(var[i], var_expect)
+
+    mu, var, deriv = gp.predict(x_test, include_nugget=False)
+
+    for i in range(2):
+
+        K = gp.emulators[i].kernel.kernel_f(x, x, theta[:-1]) + np.eye(gp.emulators[i].n)*nugget
+        Ktest = gp.emulators[i].kernel.kernel_f(x_test, x, theta[:-1])
+
+        var_expect = np.exp(theta[-2]) - np.diag(np.dot(Ktest, np.linalg.solve(K, Ktest.T)))
+
+        assert_allclose(var[i], var_expect)
+
+    gp.emulators[1].theta = None
+
+    mu, var, deriv = gp.predict(x_test, allow_not_fit=True)
+    assert np.all(np.isnan(mu[1]))
+    assert np.all(np.isnan(var[1]))
+    assert np.all(np.isnan(deriv[1]))
+
+    mu, var, deriv = gp.predict(x_test, unc=False, deriv=False, allow_not_fit=True)
+    assert np.all(np.isnan(mu[1]))
+    assert var is None
+    assert deriv is None
+
+    with pytest.raises(ValueError):
+        gp.predict(x_test, allow_not_fit=False)
+
+
 def test_MultiOutputGP_check(x, y):
     """test the methods of MultiOutputGP that extracts GPs that have or
     have not been fit (or their indices)
@@ -127,7 +226,36 @@ def test_MultiOutputGP_check(x, y):
     assert gp.get_indices_not_fit() == [0, 1]
     assert len(gp.get_emulators_fit()) == 0
     assert len(gp.get_emulators_not_fit()) == 2
-    
+
+    gp.emulators[0].fit(theta)
+
+    assert gp.get_indices_fit() == [0]
+    assert gp.get_indices_not_fit() == [1]
+    assert len(gp.get_emulators_fit()) == 1
+    assert len(gp.get_emulators_not_fit()) == 1
+
+    gp.emulators[1].fit(theta)
+
+    assert gp.get_indices_fit() == [0, 1]
+    assert gp.get_indices_not_fit() == []
+    assert len(gp.get_emulators_fit()) == 2
+    assert len(gp.get_emulators_not_fit()) == 0
+
+
+@pytest.mark.skipif(not gpu_usable(), reason=GPU_NOT_FOUND_MSG)
+def test_MultiOutputGP_GPU_check(x, y):
+    """test the methods of MultiOutputGP that extracts GPs that have or
+    have not been fit (or their indices)
+    """
+
+    gp = MultiOutputGP(x, y, nugget=0., use_gpu=True)
+    theta = np.ones(gp.emulators[0].n_params)
+
+    assert gp.get_indices_fit() == []
+    assert gp.get_indices_not_fit() == [0, 1]
+    assert len(gp.get_emulators_fit()) == 0
+    assert len(gp.get_emulators_not_fit()) == 2
+
     gp.emulators[0].fit(theta)
 
     assert gp.get_indices_fit() == [0]
diff --git a/mogp_emulator/tests/test_SequentialDesign.py b/mogp_emulator/tests/test_SequentialDesign.py
index 1950a1cd..fc1dd8e7 100644
--- a/mogp_emulator/tests/test_SequentialDesign.py
+++ b/mogp_emulator/tests/test_SequentialDesign.py
@@ -6,8 +6,12 @@
 from ..ExperimentalDesign import LatinHypercubeDesign
 from ..SequentialDesign import SequentialDesign, MICEDesign, MICEFastGP
 from ..GaussianProcess import GaussianProcess
+from ..GaussianProcessGPU import GaussianProcessGPU
 from ..fitting import fit_GP_MAP
 from tempfile import TemporaryFile
+from ..LibGPGPU import gpu_usable
+
+GPU_NOT_FOUND_MSG = "A compatible GPU could not be found or the GPU library (libgpgpu) could not be loaded"
 
 def test_SequentialDesign_init():
     "test the init method of ExperimentalDesign"
@@ -859,6 +863,7 @@ def test_MICEDesign_estimate_next_target():
     with pytest.raises(AssertionError):
         md._estimate_next_target(np.zeros(2))
 
+
 def test_MICEDesign_MICE_criterion():
     "test the function to compute the MICE criterion"
 
@@ -927,4 +932,4 @@ def test_MICEFastGP():
         result = gp.fast_predict(-1)
 
     with pytest.raises(AssertionError):
-        result = gp.fast_predict(5)
\ No newline at end of file
+        result = gp.fast_predict(5)
diff --git a/mogp_emulator/tests/test_fitting.py b/mogp_emulator/tests/test_fitting.py
index 7b155ae0..fd6318ec 100644
--- a/mogp_emulator/tests/test_fitting.py
+++ b/mogp_emulator/tests/test_fitting.py
@@ -2,9 +2,13 @@
 import pytest
 from numpy.testing import assert_allclose
 from ..GaussianProcess import GaussianProcess
+from ..GaussianProcessGPU import GaussianProcessGPU
+from ..LibGPGPU import gpu_usable
 from ..MultiOutputGP import MultiOutputGP
 from ..fitting import fit_GP_MAP, _fit_single_GP_MAP, _fit_MOGP_MAP
 
+GPU_NOT_FOUND_MSG = "A compatible GPU could not be found or the GPU library (libgpgpu) could not be loaded"
+
 def minimize_mock(fun, x0, args=(), method=None, jac=None, hess=None, hessp=None, bounds=None, constraints=(), tol=None, callback=None, options=None):
     return {"x": np.array([1.6, -2.1, -0.8]),
             "fun": 0.}
@@ -13,7 +17,7 @@ def test_fit_GP_MAP(monkeypatch):
     "test the fit_GP_MAP function"
 
     monkeypatch.setattr("mogp_emulator.fitting.minimize", minimize_mock)
-    
+
     # test correct basic functioning
 
     x = np.linspace(0., 1.)
@@ -38,6 +42,29 @@ def test_fit_GP_MAP(monkeypatch):
     assert_allclose(gp.current_logpost, logpost_exp)
 
 
+@pytest.mark.skipif(not gpu_usable(), reason=GPU_NOT_FOUND_MSG)
+def test_fit_GP_MAP_GPU(monkeypatch):
+    "test the fit_GP_MAP function"
+
+    monkeypatch.setattr("mogp_emulator.fitting.minimize", minimize_mock)
+
+    # test correct basic functioning
+
+    x = np.linspace(0., 1.)
+    y = x**2
+
+    gp = GaussianProcessGPU(x, y)
+
+    theta_exp = np.array([ 1.6, -2.1 , -0.8])
+    logpost_exp = gp.logposterior(theta_exp)
+
+    gp = fit_GP_MAP(gp)
+
+    assert isinstance(gp, GaussianProcessGPU)
+    assert_allclose(gp.theta, theta_exp)
+    assert_allclose(gp.current_logpost, logpost_exp)
+
+
 def test_fit_GP_MAP_failures():
     "test failures of fit_GP_MAP"
 
@@ -45,7 +72,7 @@ def test_fit_GP_MAP_failures():
     y = x**2
 
     gp = GaussianProcess(x, y)
-    
+
     # minimization fails
 
     with pytest.raises(RuntimeError):
@@ -73,6 +100,42 @@ def test_fit_GP_MAP_failures():
     with pytest.raises(AssertionError):
         fit_GP_MAP(gp, theta0=np.ones(1))
 
+@pytest.mark.skipif(not gpu_usable(), reason=GPU_NOT_FOUND_MSG)
+def test_fit_GP_MAP_GPU_failures():
+    "test failures of fit_GP_MAP using GaussianProcessGPU"
+
+    x = np.linspace(0., 1.)
+    y = x**2
+
+    gp = GaussianProcessGPU(x, y)
+
+    # minimization fails
+
+    with pytest.raises(RuntimeError):
+        fit_GP_MAP(gp, n_tries=1, theta0=-10000.*np.ones(3))
+
+    gp = GaussianProcessGPU(x, y, nugget=0.)
+
+    with pytest.raises(RuntimeError):
+        fit_GP_MAP(gp, n_tries=1)
+
+    with pytest.raises(RuntimeError):
+        fit_GP_MAP(gp, theta0 = np.array([800., 0., 0.]), n_tries=1)
+
+    # bad inputs
+
+    with pytest.raises(TypeError):
+        fit_GP_MAP(x)
+
+    with pytest.raises(TypeError):
+        fit_GP_MAP()
+
+    with pytest.raises(AssertionError):
+        fit_GP_MAP(gp, n_tries=-1)
+
+    with pytest.raises(AssertionError):
+        fit_GP_MAP(gp, theta0=np.ones(1))
+
 
 def test_fit_GP_MAP_MOGP():
     "test the fit_GP_MAP function with multiple outputs"
@@ -91,7 +154,7 @@ def test_fit_GP_MAP_MOGP():
     assert isinstance(gp, MultiOutputGP)
 
     np.random.seed(4335)
-    
+
     # same test, but pass args and kwargs rather than gp
 
     gp = fit_GP_MAP(x, y, mean="0.", use_patsy=False, method="L-BFGS-B", processes=1)
@@ -124,14 +187,14 @@ def test_fit_GP_MAP_MOGP():
 
 def test_fit_GP_MAP_MOGP_failures():
     "test situations where mogp fitting should fail"
-    
+
     x = np.linspace(0., 1.)
     y = np.zeros((2, 50))
     y[0] = x**2
     y[1] = 2. + x**3
 
     gp = MultiOutputGP(x, y)
-    
+
     # minimization fails
 
     gp = fit_GP_MAP(gp, n_tries=1, theta0=-10000.*np.ones(3))
@@ -190,7 +253,7 @@ def test_fit_single_GP_MAP(monkeypatch):
     "test the method to run the minimization algorithm on a GP class"
 
     monkeypatch.setattr("mogp_emulator.fitting.minimize", minimize_mock)
-    
+
     x = np.linspace(0., 1.)
     y = x**2
 
@@ -208,7 +271,7 @@ def test_fit_single_GP_MAP(monkeypatch):
 
 def test_fit_single_GP_MAP_failures():
     "test situation where fitting one emulator should fail"
-    
+
     x = np.linspace(0., 1.)
     y = x**2
 
@@ -280,7 +343,7 @@ def test_fit_MOGP_MAP():
 
     gp = _fit_MOGP_MAP(gp, theta0=np.zeros(3), n_tries=1, refit=True)
     assert not np.allclose(gp.emulators[0].theta, np.ones(3))
-    
+
 
 def test_fit_MOGP_MAP_failures():
     "test situations where fitting should fail"
@@ -291,12 +354,12 @@ def test_fit_MOGP_MAP_failures():
     y[1] = 2. + x**3
 
     gp = MultiOutputGP(x, y)
-    
+
     # minimization fails
 
     gp = _fit_MOGP_MAP(gp, n_tries=1, theta0=-10000.*np.ones(3))
     assert gp.get_indices_not_fit() == [0, 1]
-        
+
     gp = MultiOutputGP(x, y, nugget=0.)
 
     gp = _fit_MOGP_MAP(gp, theta0=np.array([0., 0., 0.]), n_tries=1)
diff --git a/mogp_gpu/test/gpu_basic.cu b/mogp_gpu/test/gpu_basic.cu
new file mode 100644
index 00000000..aebbed8a
--- /dev/null
+++ b/mogp_gpu/test/gpu_basic.cu
@@ -0,0 +1,191 @@
+#include <iostream>
+
+#include <vector>
+#include <algorithm>
+#include <string>
+#include <sstream>
+#include <assert.h>
+#include <stdexcept>
+#include <cuda_runtime.h>
+#include <cublas_v2.h>
+#include <cusolverDn.h>
+#include <thrust/device_ptr.h>
+#include <thrust/device_malloc.h>
+#include <thrust/device_free.h>
+#include <thrust/device_vector.h>
+#include <thrust/functional.h>
+#include <thrust/transform_reduce.h>
+#include <thrust/copy.h>
+
+#include "../src/cov_gpu.hpp"
+#include "../src/gp_gpu.hpp"
+
+#define WARP_SIZE 32
+#define FULL_MASK 0xffffffff
+
+typedef double REAL;
+
+void test_device_vector_copy()
+{
+    const size_t N=5;
+    // std::vector<REAL> a{1.0, 2.0, 3.0, 4.0, 5.0};
+    const double a[] = {1.0, 2.0, 3.0, 4.0, 5.0};
+    thrust::device_vector<REAL> a_d(N);
+
+    std::vector<REAL> result(N);
+    thrust::device_vector<REAL> result_d(N);
+
+    thrust::copy(a, a + N, a_d.begin());
+
+    thrust::fill(result_d.begin(), result_d.end(), 0.0);
+    thrust::transform(a_d.begin(), a_d.end(), a_d.begin(), result_d.begin(), thrust::plus<float>());
+
+    thrust::copy(result_d.begin(), result_d.end(), result.begin());
+
+    for (int i=0; i<N; i++)
+        std::cout << result[i] << " ";
+    std::cout << "\n";
+}
+
+void test_cov()
+{
+  const size_t N=3;
+  const size_t Ninput=1;
+  const size_t Ntheta=Ninput+1;
+  std::vector<REAL> x{1.0, 2.0, 3.0};
+  std::vector<REAL> result(N*N);
+
+  thrust::device_vector<REAL> result_d(N*N, 0.0);
+  thrust::device_vector<REAL> x_d(x);
+  thrust::device_vector<REAL> theta_d(Ntheta, -1.0);
+  /// calculate the covariance matrix between the inputs x
+  cov_batch_gpu(dev_ptr(result_d), N, N, Ninput, dev_ptr(x_d),
+                dev_ptr(x_d), dev_ptr(theta_d));
+
+  thrust::copy(result_d.begin(), result_d.end(), result.begin());
+  std::cout<<"Covariance matrix: "<<std::endl;
+  for (size_t i=0; i<result.size(); i++)
+    std::cout << result[i] << " ";
+  std::cout << "\n";
+}
+
+void test_cov_deriv()
+{
+    const size_t Ninput=3;
+    const size_t Ntheta=Ninput+1;
+
+    std::vector<REAL> x{1.0, 2.0, 3.0};
+    std::vector<REAL> y{4.0, 5.0, 6.0};
+
+    std::vector<REAL> result(Ntheta);
+
+    thrust::device_vector<REAL> result_d(Ntheta, 0.0);
+    thrust::device_vector<REAL> x_d(x);
+    thrust::device_vector<REAL> y_d(y);
+    thrust::device_vector<REAL> theta_d(Ntheta, -1.0);
+
+    // x, y
+    cov_deriv_theta_batch_gpu(dev_ptr(result_d), Ninput, 1, 1, dev_ptr(x_d), dev_ptr(y_d), dev_ptr(theta_d));
+    thrust::copy(result_d.begin(), result_d.end(), result.begin());
+
+    for (size_t i=0; i<Ntheta; i++)
+        std::cout << result[i] << " ";
+    std::cout << "\n";
+
+    // x, x
+    cov_deriv_theta_batch_gpu(dev_ptr(result_d), Ninput, 1, 1, dev_ptr(x_d), dev_ptr(x_d), dev_ptr(theta_d));
+    thrust::copy(result_d.begin(), result_d.end(), result.begin());
+
+    for (size_t i=0; i<Ntheta; i++)
+        std::cout << result[i] << " ";
+    std::cout << "\n";
+}
+
+void test_deriv_x() {
+  int Ninput = 1;
+  int Nx = 1;
+  int Ny = 1;
+  thrust::device_vector<REAL> result_d(Ninput, 0.);
+  thrust::device_vector<REAL> x_d(Nx, 2.);
+  thrust::device_vector<REAL> y_d(Ny, 3.);
+  std::vector<REAL> theta{-1.0, -1.0};
+  thrust::device_vector<REAL> theta_d(theta);
+
+
+    cov_deriv_x_batch_gpu(
+			  dev_ptr(result_d), Ninput, Nx, Ny, dev_ptr(x_d), dev_ptr(y_d), dev_ptr(theta_d)
+			  );
+    std::cout<<"result of test_deriv_x is "<<result_d[0]<<std::endl;
+}
+
+
+void test_sum_log_diag()
+{
+
+  size_t n=3;
+  std::vector<REAL> x{1.0, 2.0, 3.0, 4., 5., 6., 7., 8., 9.};
+  thrust::device_vector<REAL> x_d(x);
+  // determine the size of the buffer for cub::device reduce
+
+
+  // buffer for cub::DeviceReduce::Sum
+  thrust::device_vector<REAL> sum_buffer_d;
+
+    // size of sum_buffer_d
+  size_t sum_buffer_size_bytes;
+  // The following call determines the size of the cub::DeviceReduce::Sum workspace:
+  // sum_buffer_size_bytes is 0 before this call, and the size of sum_buffer_d afterwards.
+  // The end iterators are supplied but are not used.
+  cub::DeviceReduce::Sum(dev_ptr(sum_buffer_d), sum_buffer_size_bytes,
+			 sum_buffer_d.end(), sum_buffer_d.end(), n);
+  sum_buffer_d.resize(sum_buffer_size_bytes);
+
+  double result;
+  thrust::device_vector<double> result_d(1);
+  // call sum_log_diag
+  sum_log_diag(n, dev_ptr(x_d), dev_ptr(result_d), dev_ptr(sum_buffer_d), sum_buffer_size_bytes);
+  thrust::copy(result_d.begin(), result_d.end(), &result);
+
+  std::cout<<"Result of sum_log_diag is "<<result<<std::endl;
+}
+
+void test_trace()
+{
+
+  size_t n=3;
+  std::vector<REAL> x{1.0, 2.0, 3.0, 4., 5., 6., 7., 8., 9.};
+  thrust::device_vector<REAL> x_d(x);
+  // determine the size of the buffer for cub::device reduce
+
+  // buffer for cub::DeviceReduce::Sum
+  thrust::device_vector<REAL> sum_buffer_d;
+
+    // size of sum_buffer_d
+  size_t sum_buffer_size_bytes;
+  // The following call determines the size of the cub::DeviceReduce::Sum workspace:
+  // sum_buffer_size_bytes is 0 before this call, and the size of sum_buffer_d afterwards.
+  // The end iterators are supplied but are not used.
+  cub::DeviceReduce::Sum(dev_ptr(sum_buffer_d), sum_buffer_size_bytes,
+			 sum_buffer_d.end(), sum_buffer_d.end(), n);
+  sum_buffer_d.resize(sum_buffer_size_bytes);
+
+  double result;
+  thrust::device_vector<double> result_d(1);
+  // call trace
+  trace(n, dev_ptr(x_d), dev_ptr(result_d), dev_ptr(sum_buffer_d), sum_buffer_size_bytes);
+  thrust::copy(result_d.begin(), result_d.end(), &result);
+
+  std::cout<<"Result of trace is "<<result<<std::endl;
+}
+
+
+int main(void)
+{
+    test_device_vector_copy();
+    test_cov();
+    test_cov_deriv();
+    test_deriv_x();
+    test_sum_log_diag();
+    test_trace();
+    return 0;
+}
diff --git a/mogp_gpu/test/gpu_cholesky.cu b/mogp_gpu/test/gpu_cholesky.cu
new file mode 100644
index 00000000..5e4edd69
--- /dev/null
+++ b/mogp_gpu/test/gpu_cholesky.cu
@@ -0,0 +1,117 @@
+/// Example of using cuSOLVER Cholesky decomposition, based on StackOverflow
+/// https://stackoverflow.com/questions/29196139/cholesky-decomposition-with-cuda
+
+#include<iostream>
+#include <vector>
+#include <algorithm>
+#include <string>
+#include <sstream>
+#include <assert.h>
+#include <stdexcept>
+
+#include <cuda_runtime.h>
+#include <cublas_v2.h>
+#include <cusolverDn.h>
+#include <thrust/device_ptr.h>
+#include <thrust/device_malloc.h>
+#include <thrust/device_free.h>
+#include <thrust/device_vector.h>
+#include <thrust/functional.h>
+#include <thrust/transform_reduce.h>
+#include <thrust/copy.h>
+
+#define WARP_SIZE 32
+#define FULL_MASK 0xffffffff
+
+typedef double REAL;
+
+template <typename T>
+T *dev_ptr(thrust::device_vector<T>& dv)
+{
+    return dv.data().get();
+}
+
+
+/************************************/
+/*
+   Create test matrix (from Wikipedia)
+      4  12 -16       2 0 0     2 6 -8
+    ( 12 37 -43 ) = ( 6 1 0 ) ( 0 1  5 )
+     -16 -43 98      -8 5 3     0 0 3
+ */
+/************************************/
+
+void setTest3x3Matrix(double * __restrict h_A) {
+  h_A[0] = 4.0; h_A[1] = 12.0; h_A[2] = -16.0;
+  h_A[3] = 12.0; h_A[4] = 37.; h_A[5] = -43.0;
+  h_A[6] = 16.0; h_A[7] = -43.0; h_A[8] = 98.0;
+}
+
+
+void testCholesky()
+{
+
+    // --- CUDA solver initialization
+    cusolverDnHandle_t solver_handle;
+    cusolverDnCreate(&solver_handle);
+
+    // --- CUBLAS initialization
+    cublasHandle_t cublas_handle;
+    cublasCreate(&cublas_handle);
+
+    /**********************
+    Use 3x3 example from wikipedia
+    ************************/
+
+    const int N = 3;
+    double *A_h = (double *)malloc(N * N * sizeof(double));
+    setTest3x3Matrix(A_h);
+
+    // copy to device
+    thrust::device_vector<REAL> A_d(N * N);
+
+    thrust::copy(A_h, A_h + N*N, A_d.begin());
+
+    /****************************************/
+    /* COMPUTING THE CHOLESKY DECOMPOSITION */
+    /****************************************/
+    // --- cuSOLVE input/output parameters/arrays
+    int work_size = 0;
+    int info_h;
+    thrust::device_vector<int> info_d(1);
+    //   int *devInfo;           cudaMalloc(&devInfo, sizeof(int));
+
+    // --- CUDA CHOLESKY initialization
+    cusolverDnDpotrf_bufferSize(solver_handle, CUBLAS_FILL_MODE_LOWER, N, dev_ptr(A_d), N, &work_size);
+
+    std::cout<<"Buffer size is "<<work_size<<std::endl;
+    // make the buffer and resize it to the right size
+    thrust::device_vector<REAL> potrf_buffer_d;
+    potrf_buffer_d.resize(work_size);
+
+    // --- CUDA POTRF execution
+    cusolverStatus_t status = cusolverDnDpotrf(solver_handle, CUBLAS_FILL_MODE_LOWER, N, dev_ptr(A_d), N, dev_ptr(potrf_buffer_d), work_size, dev_ptr(info_d));
+
+    // check the status
+    thrust::copy(info_d.begin(), info_d.end(), &info_h);
+
+    if (status != CUSOLVER_STATUS_SUCCESS || info_h < 0) {
+      std::string msg;
+      std::stringstream smsg(msg);
+      smsg << "Error in potrf: return code " << status << ", info " << info_h;
+      throw std::runtime_error(smsg.str());
+
+    }
+
+    std::vector<double> expected = { 2.0, 6.0, -8.0, 12.0, 1.0, 5.0, 16.0, -43.0, 3.0};
+    for (unsigned int i=0; i< A_d.size(); ++i) {
+      std::cout<<"element "<<i<<" is "<<A_d[i]<<std::endl;
+      assert (A_d[i] == expected[i]);
+    }
+
+}
+
+int main() {
+  testCholesky();
+  return 0;
+}
\ No newline at end of file
diff --git a/mogp_gpu/test/gpu_test_utils.cu b/mogp_gpu/test/gpu_test_utils.cu
new file mode 100644
index 00000000..f7207dcb
--- /dev/null
+++ b/mogp_gpu/test/gpu_test_utils.cu
@@ -0,0 +1,59 @@
+#include <iostream>
+
+#include <vector>
+#include <algorithm>
+#include <string>
+#include <sstream>
+#include <assert.h>
+#include <stdexcept>
+#include <cuda_runtime.h>
+#include <cublas_v2.h>
+#include <cusolverDn.h>
+#include <thrust/device_ptr.h>
+#include <thrust/device_malloc.h>
+#include <thrust/device_free.h>
+#include <thrust/device_vector.h>
+#include <thrust/functional.h>
+#include <thrust/transform_reduce.h>
+#include <thrust/copy.h>
+
+#include "../src/util.hpp"
+
+#define WARP_SIZE 32
+#define FULL_MASK 0xffffffff
+
+typedef double REAL;
+
+template <typename T>
+T *dev_ptr(thrust::device_vector<T>& dv)
+{
+    return dv.data().get();
+}
+
+
+void test_add_diagonal()
+{
+    const size_t N=3;
+    /// create 3x3 matrix and copy to device
+    const double a[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0};
+    thrust::device_vector<REAL> a_d(N*N);
+    thrust::copy(a, a + N*N, a_d.begin());
+    // add b to the diagonal
+    const REAL b = 0.5;
+
+    add_diagonal(N, b, dev_ptr(a_d));
+    std::vector<double> expected = { 1.5, 2.0, 3.0, 4.0, 5.5, 6.0, 7.0, 8.0, 9.5};
+    for (unsigned int i=0; i< a_d.size(); ++i) {
+      std::cout<<"Element "<<i<<" is "<<a_d[i]<<std::endl;
+      assert (a_d[i] == expected[i]);
+    }
+}
+
+
+
+int main(void)
+{
+    test_add_diagonal();
+
+    return 0;
+}