diff --git a/mkldiscover.py b/mkldiscover.py
new file mode 100755
index 000000000..6423618fe
--- /dev/null
+++ b/mkldiscover.py
@@ -0,0 +1,86 @@
+#!/usr/bin/env python
+# MIT License
+#
+# Copyright (c) 2017 Anders Steen Christensen
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in all
+# copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+# SOFTWARE.
+
+from __future__ import print_function
+
+import os
+
+def mkl_exists(verbose=False):
+
+    # Get environment variables
+    __MKLROOT__ = os.environ.get('MKLROOT')
+    __LD_LIBRARY_PATH__ = os.environ.get('LD_LIBRARY_PATH')
+
+    # Check if $MKLROOT is set
+
+    if __MKLROOT__ is None:
+
+        if verbose: 
+            print("MKL-discover: MKLROOT was not set")
+
+        return False
+
+    else:
+
+        if verbose: 
+            print("MKL-discover: MKLROOT was set to")
+            print(__MKLROOT__)
+
+    # Check if path exists
+    mklroot_exists = os.path.isdir(__MKLROOT__)
+    
+    if not mklroot_exists:
+        if verbose: 
+            print("MKL-discover: MKLROOT path does not exist")
+
+        return False
+
+    found_libmkl_rt = False
+
+    # Check if libmkl_rt.so exists below $MKLROOT
+    for dirpath, dirnames, filenames in os.walk(__MKLROOT__, followlinks=True):
+
+        if "libmkl_rt.so" in filenames:
+
+            if verbose:
+                print("MKL-discover: Found libmkl_rt.so at ", dirpath)
+
+            # Check that the dirpath where libmkl_rt.so is in $LD_LIBRARY_PATH
+            if dirpath in __LD_LIBRARY_PATH__:
+
+                if verbose:
+                    print("MKL-discover: Found %s in $LD_LIBRARY_PATH" % dirpath)
+                    print("MKL-discover: Concluding that MKL can be used.")
+                found_libmkl_rt = True
+
+
+    return found_libmkl_rt
+
+if __name__ == "__main__":
+
+    mkl_present = mkl_exists(verbose=False)
+
+    if mkl_present:
+        print("MKL could be found")
+    else:
+        print("MKL could NOT be found")
diff --git a/qml/fcho_solve.f90 b/qml/fcho_solve.f90
index fc7c1fee8..4b0add806 100644
--- a/qml/fcho_solve.f90
+++ b/qml/fcho_solve.f90
@@ -46,23 +46,100 @@ subroutine fcho_solve(A,y,x)
 
 end subroutine fcho_solve
 
-! subroutine fcho_invert(A)
-! 
-!     implicit none
-! 
-!     double precision, dimension(:,:), intent(inout) :: A
-!     integer :: info, na
-! 
-!     na = size(A, dim=1)
-! 
-!     call dpotrf("L", na, A , na, info)
-!     if (info > 0) then
-!         write (*,*) "WARNING: Cholesky decomposition DPOTRF() exited with error code:", info
-!     endif
-! 
-!     call dpotri("L", na, A , na, info )
-!     if (info > 0) then
-!         write (*,*) "WARNING: Cholesky inversion DPOTRI() exited with error code:", info
-!     endif
-! 
-! end subroutine fcho_invert
+subroutine fcho_invert(A)
+
+    implicit none
+
+    double precision, dimension(:,:), intent(inout) :: A
+    integer :: info, na
+
+    na = size(A, dim=1)
+
+    call dpotrf("L", na, A , na, info)
+    if (info > 0) then
+        write (*,*) "WARNING: Cholesky decomposition DPOTRF() exited with error code:", info
+    endif
+
+    call dpotri("L", na, A , na, info )
+    if (info > 0) then
+        write (*,*) "WARNING: Cholesky inversion DPOTRI() exited with error code:", info
+    endif
+
+end subroutine fcho_invert
+
+
+subroutine fbkf_invert(A)
+
+    implicit none
+
+    double precision, dimension(:,:), intent(inout) :: A
+    integer :: info, na, nb
+
+    integer, dimension(size(A,1)) :: ipiv   ! pivot indices
+    integer :: ilaenv
+    
+    integer :: lwork
+
+    double precision, allocatable, dimension(:) :: work
+
+    na = size(A, dim=1)
+    
+    nb = ilaenv( 1, 'DSYTRF', "L", na, -1, -1, -1 )
+
+    lwork = na*nb
+
+    allocate(work(lwork))
+
+    ! call dpotrf("L", na, A , na, info)
+    call dsytrf("L", na, A, na, ipiv, work, lwork, info)
+    if (info > 0) then
+        write (*,*) "WARNING: Bunch-Kaufman factorization DSYTRI() exited with error code:", info
+    endif
+
+    ! call dpotri("L", na, A , na, info )
+    call dsytri( "L", na, a, na, ipiv, work, info )
+    if (info > 0) then
+        write (*,*) "WARNING: BKF inversion DPOTRI() exited with error code:", info
+    endif
+
+    deallocate(work)
+
+end subroutine fbkf_invert
+
+
+subroutine fbkf_solve(A,y,x)
+
+    implicit none
+
+    double precision, dimension(:,:), intent(in) :: A
+    double precision, dimension(:), intent(in) :: y
+    double precision, dimension(:), intent(inout) :: x
+
+    double precision, allocatable, dimension(:) :: work
+    integer :: ilaenv
+
+    integer, dimension(size(A,1)) :: ipiv   ! pivot indices
+    integer :: info, na, nb, lwork
+
+    na = size(A, dim=1)
+
+    nb = ilaenv( 1, 'DSYTRF', "L", na, -1, -1, -1 )
+
+    lwork = na*nb
+    allocate(work(lwork))
+
+    call dsytrf("L", na, A, na, ipiv, work, lwork, info)
+    if (info > 0) then
+        write (*,*) "WARNING: Bunch-Kaufman factorization DSYTRI() exited with error code:", info
+    endif
+
+    x(:na) = y(:na)
+
+    call dsytrs("L", na, 1, A, na, ipiv, x, na, info )
+
+    if (info > 0) then
+        write (*,*) "WARNING: Bunch-Kaufman solver DSYTRS() exited with error code:", info
+    endif
+
+    deallocate(work)
+end subroutine fbkf_solve
diff --git a/qml/math.py b/qml/math.py
index 825c67c95..500f9358e 100644
--- a/qml/math.py
+++ b/qml/math.py
@@ -22,45 +22,56 @@
 
 import numpy as np
 
+from copy import deepcopy
+
 from .fcho_solve import fcho_solve
-# from .fcho_solve import fcho_invert
-
-
-# Disabled due to bug.
-# def cho_invert(A):
-#     """ Solves [A x = y] for x using a Cholesky decomposition
-#         via calls to LAPACK dpotrf and dpotri in the F2PY module.
-# 
-#         Arguments:
-#         ==============
-#         A -- the A-matrix (symmetric and positive definite).
-# 
-#         Returns:
-#         ==============
-#         A -- the inverted A-matrix
-#     """
-# 
-#     B = np.asfortranarray(A)
-#     fcho_invert(B)
-# 
-#     return B
+from .fcho_solve import fcho_invert
+from .fcho_solve import fbkf_solve
+from .fcho_solve import fbkf_invert
+
+
+def cho_invert(A):
+    """ Returns the inverse of a positive definite matrix, using a Cholesky decomposition
+        via calls to LAPACK dpotrf and dpotri in the F2PY module.
+
+        :param A: Matrix (symmetric and positive definite, left-hand side).
+        :type A: numpy array
+
+        :return: The inverse matrix
+        :rtype: numpy array
+    """
+
+    if len(A.shape) != 2 or A.shape[0] != A.shape[1]:
+        raise ValueError('expected square matrix')
+
+    I = np.asfortranarray(A)
+
+    fcho_invert(I)
+
+    # Matrix to store the inverse
+    i_lower = np.tril_indices_from(A)
+
+    # Copy lower triangle to upper
+    I.T[i_lower] = I[i_lower]
+
+    return I
 
 
 def cho_solve(A, y):
     """ Solves the equation
 
-            :math:`A x = y` 
+            :math:`A x = y`
 
-        for x using a Cholesky decomposition  via calls to LAPACK dpotrf and dpotrs in the F2PY module.
+        for x using a Cholesky decomposition  via calls to LAPACK dpotrf and dpotrs in the F2PY module. Preserves the input matrix A.
 
         :param A: Matrix (symmetric and positive definite, left-hand side).
-        :type A: numpy array 
+        :type A: numpy array
         :param y: Vector (right-hand side of the equation).
-        :type y: numpy array 
+        :type y: numpy array
 
         :return: The solution vector.
-        :rtype: numpy array 
-    """
+        :rtype: numpy array
+        """
 
     if len(A.shape) != 2 or A.shape[0] != A.shape[1]:
         raise ValueError('expected square matrix')
@@ -70,7 +81,84 @@ def cho_solve(A, y):
 
     n = A.shape[0]
 
+    # Backup diagonal before Cholesky-decomposition
+    A_diag = A[np.diag_indices_from(A)]
+
     x = np.zeros(n)
     fcho_solve(A, y, x)
 
+    # Reset diagonal after Cholesky-decomposition
+    A[np.diag_indices_from(A)] = A_diag
+
+    # Copy lower triangle to upper
+    i_lower = np.tril_indices_from(A)
+    A.T[i_lower] = A[i_lower]
+
+    return x
+
+
+def bkf_invert(A):
+    """ Returns the inverse of a positive definite matrix, using a Cholesky decomposition
+        via calls to LAPACK dpotrf and dpotri in the F2PY module.
+
+        :param A: Matrix (symmetric and positive definite, left-hand side).
+        :type A: numpy array
+
+        :return: The inverse matrix
+        :rtype: numpy array
+    """
+
+    if len(A.shape) != 2 or A.shape[0] != A.shape[1]:
+        raise ValueError('expected square matrix')
+
+    I = np.asfortranarray(A)
+
+    fbkf_invert(I)
+
+    # Matrix to store the inverse
+    i_lower = np.tril_indices_from(A)
+
+    # Copy lower triangle to upper
+    I.T[i_lower] = I[i_lower]
+
+    return I
+
+
+def bkf_solve(A, y):
+    """ Solves the equation
+
+            :math:`A x = y`
+
+        for x using a Cholesky decomposition  via calls to LAPACK dpotrf and dpotrs in the F2PY module. Preserves the input matrix A.
+
+        :param A: Matrix (symmetric and positive definite, left-hand side).
+        :type A: numpy array
+        :param y: Vector (right-hand side of the equation).
+        :type y: numpy array
+
+        :return: The solution vector.
+        :rtype: numpy array
+        """
+
+    if len(A.shape) != 2 or A.shape[0] != A.shape[1]:
+        raise ValueError('expected square matrix')
+
+    if len(y.shape) != 1 or y.shape[0] != A.shape[1]:
+        raise ValueError('expected matrix and vector of same stride size')
+
+    n = A.shape[0]
+
+    # Backup diagonal before Cholesky-decomposition
+    A_diag = A[np.diag_indices_from(A)]
+
+    x = np.zeros(n)
+    fbkf_solve(A, y, x)
+
+    # Reset diagonal after Cholesky-decomposition
+    A[np.diag_indices_from(A)] = A_diag
+
+    # Copy lower triangle to upper
+    i_lower = np.tril_indices_from(A)
+    A.T[i_lower] = A[i_lower]
+
     return x
diff --git a/setup.py b/setup.py
index f677bbf1e..b99f40468 100755
--- a/setup.py
+++ b/setup.py
@@ -1,6 +1,8 @@
 import sys
 from numpy.distutils.core import Extension, setup
 
+from mkldiscover import mkl_exists
+
 __author__ = "Anders S. Christensen"
 __copyright__ = "Copyright 2016"
 __credits__ = ["Anders S. Christensen (2016) https://github.com/qmlcode/qml"]
@@ -21,6 +23,11 @@
 LINKER_FLAGS = ["-lgomp"]
 MATH_LINKER_FLAGS = ["-lblas", "-llapack"]
 
+# UNCOMMENT TO FORCE LINKING TO MKL with GNU compilers:
+if mkl_exists(verbose=True):
+    LINKER_FLAGS = ["-lgomp", " -lpthread", "-lm", "-ldl"]
+    MATH_LINKER_FLAGS = ["-L${MKLROOT}/lib/intel64", "-lmkl_rt"]
+
 # For clang without OpenMP: (i.e. most Apple/mac system)
 if sys.platform == "darwin" and all(["gnu" not in arg for arg in sys.argv]):
     COMPILER_FLAGS = ["-O3", "-m64", "-march=native", "-fPIC"]
@@ -35,9 +42,6 @@
     MATH_LINKER_FLAGS = ["-L${MKLROOT}/lib/intel64", "-lmkl_rt"]
 
 
-# UNCOMMENT TO FORCE LINKING TO MKL with GNU compilers:
-# LINKER_FLAGS = ["-lgomp", " -lpthread", "-lm", "-ldl"]
-# MATH_LINKER_FLAGS = ["-L${MKLROOT}/lib/intel64", "-lmkl_rt"]
 
 
 ext_farad_kernels = Extension(name = 'farad_kernels',
diff --git a/tests/data/y_cho_solve.txt b/tests/data/y_cho_solve.txt
new file mode 100644
index 000000000..3e692b133
--- /dev/null
+++ b/tests/data/y_cho_solve.txt
@@ -0,0 +1,200 @@
+-1.620680000000000064e+03
+-1.623029999999999973e+03
+-1.314670000000000073e+03
+-1.757589999999999918e+03
+-1.423589999999999918e+03
+-1.503670000000000073e+03
+-1.606660000000000082e+03
+-1.374029999999999973e+03
+-1.551160000000000082e+03
+-1.565059999999999945e+03
+-1.362329999999999927e+03
+-1.751089999999999918e+03
+-1.274170000000000073e+03
+-1.484069999999999936e+03
+-1.619380000000000109e+03
+-1.354339999999999918e+03
+-1.122490000000000009e+03
+-1.618779999999999973e+03
+-1.015860000000000014e+03
+-1.597410000000000082e+03
+-1.631990000000000009e+03
+-1.215059999999999945e+03
+-1.475160000000000082e+03
+-1.223019999999999982e+03
+-1.481029999999999973e+03
+-1.468859999999999900e+03
+-1.419259999999999991e+03
+-1.427509999999999991e+03
+-1.552740000000000009e+03
+-8.690960000000000036e+02
+-1.319079999999999927e+03
+-1.550259999999999991e+03
+-1.312400000000000091e+03
+-1.201980000000000018e+03
+-1.749150000000000091e+03
+-1.274930000000000064e+03
+-1.246619999999999891e+03
+-1.018419999999999959e+03
+-1.408779999999999973e+03
+-1.892740000000000009e+03
+-1.339890000000000100e+03
+-1.618609999999999900e+03
+-1.615450000000000045e+03
+-8.299379999999999882e+02
+-1.427359999999999900e+03
+-1.218000000000000000e+03
+-1.387460000000000036e+03
+-1.617269999999999982e+03
+-1.405730000000000018e+03
+-1.485869999999999891e+03
+-1.412180000000000064e+03
+-1.542200000000000045e+03
+-1.292559999999999945e+03
+-1.160720000000000027e+03
+-1.214240000000000009e+03
+-1.604859999999999900e+03
+-1.156960000000000036e+03
+-1.427509999999999991e+03
+-1.190119999999999891e+03
+-1.593500000000000000e+03
+-1.193900000000000091e+03
+-1.483720000000000027e+03
+-1.567079999999999927e+03
+-1.232490000000000009e+03
+-9.582459999999999809e+02
+-1.313950000000000045e+03
+-1.324789999999999964e+03
+-1.236809999999999945e+03
+-1.304039999999999964e+03
+-1.564019999999999982e+03
+-1.366839999999999918e+03
+-1.246589999999999918e+03
+-1.358829999999999927e+03
+-1.089440000000000055e+03
+-1.094210000000000036e+03
+-1.262619999999999891e+03
+-1.254670000000000073e+03
+-1.073349999999999909e+03
+-1.259480000000000018e+03
+-1.460960000000000036e+03
+-1.548589999999999918e+03
+-1.065210000000000036e+03
+-1.489920000000000073e+03
+-1.502990000000000009e+03
+-1.613589999999999918e+03
+-1.247680000000000064e+03
+-1.063880000000000109e+03
+-7.240489999999999782e+02
+-1.487309999999999945e+03
+-1.354900000000000091e+03
+-1.223349999999999909e+03
+-1.314420000000000073e+03
+-1.208299999999999955e+03
+-1.316880000000000109e+03
+-1.075410000000000082e+03
+-1.357829999999999927e+03
+-1.242680000000000064e+03
+-1.692319999999999936e+03
+-1.621019999999999982e+03
+-9.469809999999999945e+02
+-1.614200000000000045e+03
+-1.294980000000000018e+03
+-1.239869999999999891e+03
+-1.491369999999999891e+03
+-1.221779999999999973e+03
+-1.274450000000000045e+03
+-1.192150000000000091e+03
+-1.102049999999999955e+03
+-1.626779999999999973e+03
+-1.472829999999999927e+03
+-1.746630000000000109e+03
+-1.621529999999999973e+03
+-1.622410000000000082e+03
+-1.417619999999999891e+03
+-8.701760000000000446e+02
+-1.268380000000000109e+03
+-1.603259999999999991e+03
+-1.540599999999999909e+03
+-1.225049999999999955e+03
+-1.766819999999999936e+03
+-1.617650000000000091e+03
+-1.472019999999999982e+03
+-1.174049999999999955e+03
+-1.455650000000000091e+03
+-1.230190000000000055e+03
+-1.419059999999999945e+03
+-1.307369999999999891e+03
+-1.316980000000000018e+03
+-1.213740000000000009e+03
+-1.226049999999999955e+03
+-1.341470000000000027e+03
+-1.332970000000000027e+03
+-1.228190000000000055e+03
+-1.384039999999999964e+03
+-1.228829999999999927e+03
+-1.628349999999999909e+03
+-1.489329999999999927e+03
+-1.451309999999999945e+03
+-1.506309999999999945e+03
+-1.422140000000000100e+03
+-1.007370000000000005e+03
+-1.319789999999999964e+03
+-1.024740000000000009e+03
+-1.429380000000000109e+03
+-1.153009999999999991e+03
+-1.232880000000000109e+03
+-1.269279999999999973e+03
+-1.057710000000000036e+03
+-1.182099999999999909e+03
+-1.194319999999999936e+03
+-1.613680000000000064e+03
+-1.167599999999999909e+03
+-1.348190000000000055e+03
+-1.374549999999999955e+03
+-1.224140000000000100e+03
+-9.713200000000000500e+02
+-1.158740000000000009e+03
+-1.629089999999999918e+03
+-1.758170000000000073e+03
+-1.104319999999999936e+03
+-1.121480000000000018e+03
+-1.697569999999999936e+03
+-1.265269999999999982e+03
+-1.364519999999999982e+03
+-1.689690000000000055e+03
+-9.284790000000000418e+02
+-1.612400000000000091e+03
+-1.124779999999999973e+03
+-1.098869999999999891e+03
+-1.462680000000000064e+03
+-9.756960000000000264e+02
+-1.348420000000000073e+03
+-1.541660000000000082e+03
+-1.242730000000000018e+03
+-1.255970000000000027e+03
+-1.408680000000000064e+03
+-1.487529999999999973e+03
+-1.480559999999999945e+03
+-8.064950000000000045e+02
+-1.106400000000000091e+03
+-1.549269999999999982e+03
+-1.369900000000000091e+03
+-1.510900000000000091e+03
+-1.603650000000000091e+03
+-1.378970000000000027e+03
+-1.325579999999999927e+03
+-1.355849999999999909e+03
+-1.684599999999999909e+03
+-1.243000000000000000e+03
+-1.546710000000000036e+03
+-1.185190000000000055e+03
+-1.757579999999999927e+03
+-1.611410000000000082e+03
+-1.471180000000000064e+03
+-1.616240000000000009e+03
+-1.295740000000000009e+03
+-1.748529999999999973e+03
+-9.805030000000000427e+02
+-1.454150000000000091e+03
+-1.121940000000000055e+03
diff --git a/tests/test_math.py b/tests/test_math.py
index f51d55fc9..5cb461ddc 100644
--- a/tests/test_math.py
+++ b/tests/test_math.py
@@ -20,23 +20,100 @@
 # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 # SOFTWARE.
 
+import os
+
 import numpy as np
 
+from copy import deepcopy
+
 import qml
 import qml.math
 
-def test_cholesky():
+def test_cho_solve():
 
-    A = np.asarray([[ 2.0, -1.0,  0.0],
-                  [-1.0,  2.0, -1.0],
-                  [ 0.0, -1.0,  2.0]])
+    test_dir = os.path.dirname(os.path.realpath(__file__))
 
-    y = np.asarray([1.0, 1.0, 1.0])
+    A_ref = np.loadtxt(test_dir + "/data/K_local_gaussian.txt")
+    y_ref = np.loadtxt(test_dir + "/data/y_cho_solve.txt")
 
+    A = deepcopy(A_ref) 
+    y = deepcopy(y_ref) 
     x_qml   = qml.math.cho_solve(A,y)
+
+    # Check arrays are unchanged
+    assert np.allclose(y, y_ref)
+    assert np.allclose(A, A_ref)
+
+    A = deepcopy(A_ref) 
     x_scipy = np.linalg.solve(A, y)
 
+    # Check for correct solution
     assert np.allclose(x_qml, x_scipy)
 
+
+def test_cho_invert():
+
+    test_dir = os.path.dirname(os.path.realpath(__file__))
+
+    A_ref = np.loadtxt(test_dir + "/data/K_local_gaussian.txt")
+
+    A = deepcopy(A_ref) 
+    Ai_qml = qml.math.cho_invert(A)
+
+    # Check A is unchanged
+    assert np.allclose(A, A_ref)
+
+    A = deepcopy(A_ref) 
+    one = np.eye(A.shape[0])
+    
+    # Check that it is a true inverse
+    assert np.allclose(np.matmul(A, Ai_qml), one, atol=1e-7)
+
+
+def test_bkf_invert():
+
+    test_dir = os.path.dirname(os.path.realpath(__file__))
+
+    A_ref = np.loadtxt(test_dir + "/data/K_local_gaussian.txt")
+
+    A = deepcopy(A_ref) 
+    Ai_qml = qml.math.bkf_invert(A)
+
+    # Check A is unchanged
+    assert np.allclose(A, A_ref)
+
+    A = deepcopy(A_ref) 
+    one = np.eye(A.shape[0])
+   
+    np.set_printoptions(linewidth=20000)
+    assert np.allclose(np.matmul(A, Ai_qml), one, atol=1e-7)
+
+
+def test_bkf_solve():
+
+    test_dir = os.path.dirname(os.path.realpath(__file__))
+
+    A_ref = np.loadtxt(test_dir + "/data/K_local_gaussian.txt")
+    y_ref = np.loadtxt(test_dir + "/data/y_cho_solve.txt")
+
+    A = deepcopy(A_ref) 
+    y = deepcopy(y_ref)
+    x_qml   = qml.math.bkf_solve(A,y)
+
+    # Check arrays are unchanged
+    assert np.allclose(y, y_ref)
+    assert np.allclose(A, A_ref)
+
+    A = deepcopy(A_ref) 
+    y = deepcopy(y_ref) 
+    x_scipy = np.linalg.solve(A, y)
+
+    # Check for correct solution
+    assert np.allclose(x_qml, x_scipy)
+
+
 if __name__ == "__main__":
-    test_cholesky()
+    test_cho_solve()
+    test_cho_invert()
+    test_bkf_invert()
+    test_bkf_solve()