FEAT/TEST: Adding all functions from lapack.h

Author: pavanky

arrayfire/__init__.py

@@ -15,3 +15,4 @@
 from .blas import *
 from .arith import *
 from .statistics import *
+from .lapack import *

arrayfire/blas.py

@@ -16,6 +16,24 @@ def matmul(lhs, rhs, lhs_opts=AF_MAT_NONE, rhs_opts=AF_MAT_NONE):
                              lhs_opts, rhs_opts))
     return out
 
+def matmulTN(lhs, rhs):
+    out = array()
+    safe_call(clib.af_matmul(pointer(out.arr), lhs.arr, rhs.arr,\
+                             AF_MAT_TRANS, AF_MAT_NONE))
+    return out
+
+def matmulNT(lhs, rhs):
+    out = array()
+    safe_call(clib.af_matmul(pointer(out.arr), lhs.arr, rhs.arr,\
+                             AF_MAT_NONE, AF_MAT_TRANS))
+    return out
+
+def matmulTT(lhs, rhs):
+    out = array()
+    safe_call(clib.af_matmul(pointer(out.arr), lhs.arr, rhs.arr,\
+                             AF_MAT_TRANS, AF_MAT_TRANS))
+    return out
+
 def dot(lhs, rhs, lhs_opts=AF_MAT_NONE, rhs_opts=AF_MAT_NONE):
     out = array()
     safe_call(clib.af_dot(pointer(out.arr), lhs.arr, rhs.arr,\

arrayfire/lapack.py

@@ -0,0 +1,80 @@
+#######################################################
+# Copyright (c) 2014, ArrayFire
+# All rights reserved.
+#
+# This file is distributed under 3-clause BSD license.
+# The complete license agreement can be obtained at:
+# http://arrayfire.com/licenses/BSD-3-Clause
+########################################################
+
+from .library import *
+from .array import *
+
+def lu(A):
+    L = array()
+    U = array()
+    P = array()
+    safe_call(clib.af_lu(pointer(L.arr), pointer(U.arr), pointer(P.arr), A.arr))
+    return L,U,P
+
+def lu_inplace(A, pivot="lapack"):
+    P = array()
+    is_pivot_lapack = False if (pivot == "full") else True
+    safe_call(clib.af_lu_inplace(pointer(P.arr), A.arr, is_pivot_lapack))
+    return P
+
+def qr(A):
+    Q = array()
+    R = array()
+    T = array()
+    safe_call(clib.af_lu(pointer(Q.arr), pointer(R.arr), pointer(T.arr), A.arr))
+    return Q,R,T
+
+def qr_inplace(A):
+    T = array()
+    safe_call(clib.af_qr_inplace(pointer(T.arr), A.arr))
+    return T
+
+def cholesky(A, is_upper=True):
+    R = array()
+    info = c_int(0)
+    safe_call(clib.af_cholesky(pointer(R.arr), pointer(info), A.arr, is_upper))
+    return R, info.value
+
+def cholesky_inplace(A, is_upper=True):
+    info = c_int(0)
+    safe_call(clib.af_cholesky_inplace(pointer(info), A.arr, is_upper))
+    return info.value
+
+def solve(A, B, options=AF_MAT_NONE):
+    X = array()
+    safe_call(clib.af_solve(pointer(X.arr), A.arr, B.arr, options))
+    return X
+
+def solve_lu(A, P, B, options=AF_MAT_NONE):
+    X = array()
+    safe_call(clib.af_solve_lu(pointer(X.arr), A.arr, P.arr, B.arr, options))
+    return X
+
+def inverse(A, options=AF_MAT_NONE):
+    I = array()
+    safe_call(clib.af_inverse(pointer(I.arr), A.arr, options))
+    return I
+
+def rank(A, tol=1E-5):
+    r = c_uint(0)
+    safe_call(clib.af_rank(pointer(r), A.arr, c_double(tol)))
+    return r.value
+
+def det(A):
+    re = c_double(0)
+    im = c_double(0)
+    safe_call(clib.af_det(pointer(re), pointer(im), A.arr))
+    re = re.value
+    im = im.value
+    return re if (im == 0) else re + im * 1j
+
+def norm(A, norm_type=AF_NORM_EUCLID, p=1.0, q=1.0):
+    res = c_double(0)
+    safe_call(clib.af_norm(pointer(res), A.arr, norm_type, c_double(p), c_double(q)))
+    return res.value

tests/simple_lapack.py

@@ -0,0 +1,71 @@
+#!/usr/bin/python
+#######################################################
+# Copyright (c) 2014, ArrayFire
+# All rights reserved.
+#
+# This file is distributed under 3-clause BSD license.
+# The complete license agreement can be obtained at:
+# http://arrayfire.com/licenses/BSD-3-Clause
+########################################################
+import arrayfire as af
+
+a = af.randu(5,5)
+
+l,u,p = af.lu(a)
+
+af.print_array(l)
+af.print_array(u)
+af.print_array(p)
+
+p = af.lu_inplace(a, "full")
+
+af.print_array(a)
+af.print_array(p)
+
+a = af.randu(5,3)
+
+q,r,t = af.qr(a)
+
+af.print_array(q)
+af.print_array(r)
+af.print_array(t)
+
+af.qr_inplace(a)
+
+af.print_array(a)
+
+a = af.randu(5, 5)
+a = af.matmulTN(a, a) + 10 * af.identity(5,5)
+
+R,info = af.cholesky(a)
+af.print_array(R)
+print(info)
+
+af.cholesky_inplace(a)
+af.print_array(a)
+
+a = af.randu(5,5)
+ai = af.inverse(a)
+
+af.print_array(a)
+af.print_array(ai)
+
+x0 = af.randu(5, 3)
+b = af.matmul(a, x0)
+x1 = af.solve(a, b)
+
+af.print_array(x0)
+af.print_array(x1)
+
+p = af.lu_inplace(a)
+
+x2 = af.solve_lu(a, p, b)
+
+af.print_array(x2)
+
+print(af.rank(a))
+print(af.det(a))
+print(af.norm(a, af.AF_NORM_EUCLID))
+print(af.norm(a, af.AF_NORM_MATRIX_1))
+print(af.norm(a, af.AF_NORM_MATRIX_INF))
+print(af.norm(a, af.AF_NORM_MATRIX_L_PQ, 1, 1))