Initial commit

Author: dwfmarchant

.gitignore

@@ -0,0 +1,64 @@
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*$py.class
+
+# C extensions
+*.so
+
+# Distribution / packaging
+.Python
+env/
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+*.egg-info/
+.installed.cfg
+*.egg
+
+# PyInstaller
+#  Usually these files are written by a python script from a template
+#  before PyInstaller builds the exe, so as to inject date/other infos into it.
+*.manifest
+*.spec
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.coverage
+.coverage.*
+.cache
+nosetests.xml
+coverage.xml
+*,cover
+.hypothesis/
+
+# Translations
+*.mo
+*.pot
+
+# Django stuff:
+*.log
+
+# Sphinx documentation
+docs/_build/
+
+# PyBuilder
+target/
+
+#Ipython Notebook
+.ipynb_checkpoints
+*.sublime-project
+*.sublime-workspace

LICENSE

@@ -0,0 +1,21 @@
+The MIT License (MIT)
+
+Copyright (c) 2016 David Marchant
+
+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.

README.md

@@ -0,0 +1,2 @@
+# pyMKLwrap
+Python wrapper of Intel MKL routines

pyMKL/__init__.py

@@ -0,0 +1,10 @@
+import numpy as np
+import scipy.sparse as sp
+from ctypes import CDLL, cdll, RTLD_GLOBAL
+from ctypes import POINTER, byref, c_int, c_longlong
+
+path = 'libmkl_intel_lp64.dylib'
+MKLlib = CDLL(path, RTLD_GLOBAL)
+
+from pardisoInterface import pardisoinit, pardiso
+from pardisoSolver import pardisoSolver

pyMKL/pardisoInterface.py

@@ -0,0 +1,28 @@
+from pyMKL import MKLlib
+from ctypes import POINTER, c_int, c_longlong
+
+pardisoinit = MKLlib.pardisoinit
+
+pardisoinit.argtypes = [POINTER(c_longlong), 
+                        POINTER(c_int), 
+                        POINTER(c_int)]
+pardisoinit.restype = None
+
+pardiso = MKLlib.pardiso
+pardiso.argtypes = [POINTER(c_longlong), # pt
+                    POINTER(c_int),      # maxfct
+                    POINTER(c_int),      # mnum
+                    POINTER(c_int),      # mtype
+                    POINTER(c_int),      # phase
+                    POINTER(c_int),      # n
+                    POINTER(None),       # a
+                    POINTER(c_int),      # ia
+                    POINTER(c_int),      # ja
+                    POINTER(c_int),      # perm
+                    POINTER(c_int),      # nrhs
+                    POINTER(c_int),      # iparm
+                    POINTER(c_int),      # msglvl
+                    POINTER(None),       # b
+                    POINTER(None),       # x
+                    POINTER(c_int)]      # error)
+pardiso.restype = None

pyMKL/pardisoSolver.py

@@ -0,0 +1,127 @@
+from pyMKL import pardisoinit, pardiso
+from ctypes import POINTER, byref, c_longlong, c_int, c_float
+import numpy as np
+
+"""
+mtype options
+1 -> real and structurally symmetric
+2 -> real and symmetric positive definite
+-2 -> real and symmetric indefinite
+3 -> complex and structurally symmetric
+4 -> complex and Hermitian positive definite
+-4 -> complex and Hermitian indefinite
+6 -> complex and symmetric
+11 -> real and nonsymmetric
+13 -> complex and nonsymmetric
+
+
+phase options
+11 -> Analysis
+12 -> Analysis, numerical factorization
+13 -> Analysis, numerical factorization, solve, iterative refinement
+22 -> Numerical factorization
+23 -> Numerical factorization, solve, iterative refinement
+33 -> Solve, iterative refinement
+331 -> like phase=33, but only forward substitution
+332 -> like phase=33, but only diagonal substitution (if available)
+333 -> like phase=33, but only backward substitution
+0 -> Release internal memory for L and U matrix number mnum
+-1 -> Release all internal memory for all matrices
+"""
+
+class pardisoSolver(object):
+    """docstring for pardisoSolver"""
+    def __init__(self, A, mtype=11):
+        
+        # Prep the matrix
+        self.A = A.tocsr()
+        self.a = A.data
+        self.ia = A.indptr
+        self.ja = A.indices
+
+        self._MKL_a = self.a.ctypes.data_as(POINTER(c_float))
+        self._MKL_ia = self.ia.ctypes.data_as(POINTER(c_int))
+        self._MKL_ja = self.ja.ctypes.data_as(POINTER(c_int))
+
+        self.mtype = mtype
+        self.n = A.shape[0]
+
+        # Hardcode some parameters for now...
+        self.maxfct = 1
+        self.mnum = 1
+        self.perm = 0
+        self.msglvl = 1
+
+        # Initialize handle to data structure
+        self.pt = np.zeros(64, np.int64)
+        self._MKL_pt = self.pt.ctypes.data_as(POINTER(c_longlong))
+
+        # Initialize parameters
+        self.iparm = np.zeros(64, dtype=np.int32)
+        self._MKL_iparm = self.iparm.ctypes.data_as(POINTER(c_int))
+
+        # Initialize pardiso
+        pardisoinit(self._MKL_pt, byref(c_int(self.mtype)), self._MKL_iparm)
+
+        # Set iparm
+        self.iparm[1] = 3 # Use parallel nested dissection for reordering
+        self.iparm[23] = 1 # Use parallel factorization
+        self.iparm[34] = 1 # Zero base indexing
+
+
+    def run_pardiso(self, phase, rhs=None):
+        
+        if rhs is None:
+            nrhs = 0
+            x = np.zeros(1)
+            rhs = np.zeros(1)
+        else:
+            nrhs = 1
+            x = np.zeros(self.n)
+
+        MKL_rhs = rhs.ctypes.data_as(POINTER(c_float))
+        MKL_x = x.ctypes.data_as(POINTER(c_float))
+        ERR = 0
+
+        pardiso(self._MKL_pt,               # pt
+                byref(c_int(self.maxfct)),  # maxfct
+                byref(c_int(self.mnum)),    # mnum
+                byref(c_int(self.mtype)),   # mtype
+                byref(c_int(phase)),        # phase
+                byref(c_int(self.n)),       # n
+                self._MKL_a,                # a
+                self._MKL_ia,               # ia
+                self._MKL_ja,               # ja
+                byref(c_int(self.perm)),    # perm
+                byref(c_int(nrhs)),         # nrhs
+                self._MKL_iparm,            # iparm
+                byref(c_int(self.msglvl)),  # msglvl
+                MKL_rhs,                    # b
+                MKL_x,                      # x
+                byref(c_int(ERR)))          # error
+
+        return x
+    
+
+
+if __name__ == '__main__':
+    
+    import scipy.sparse as sp
+
+    nSize = 20
+    A = sp.rand(nSize, nSize, 0.05, format='csr', random_state=10)
+    A = A.T.dot(A) + sp.spdiags(np.ones(nSize), 0, nSize, nSize)
+    A = A.tocsr()
+
+    np.random.seed(1)
+    xTrue = np.random.rand(nSize)
+    rhs = A.dot(xTrue)
+
+    pSolve = pardisoSolver(A)
+
+    pSolve.run_pardiso(12)
+    x = pSolve.run_pardiso(33, rhs)
+
+    print np.linalg.norm(x-xTrue)/np.linalg.norm(xTrue)
+
+