Linalg

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")

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

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