#######################################################
# Copyright (c) 2015, 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
########################################################
"""
Dense Linear Algebra functions (solve, inverse, etc).
"""
from .library import *
from .array import *
[docs]def lu(A):
"""
LU decomposition.
Parameters
----------
A: af.Array
A 2 dimensional arrayfire array.
Returns
-------
(L,U,P): tuple of af.Arrays
- L - Lower triangular matrix.
- U - Upper triangular matrix.
- P - Permutation array.
Note
----
The original matrix `A` can be reconstructed using the outputs in the following manner.
>>> A[P, :] = af.matmul(L, U)
"""
L = Array()
U = Array()
P = Array()
safe_call(backend.get().af_lu(ct.pointer(L.arr), ct.pointer(U.arr), ct.pointer(P.arr), A.arr))
return L,U,P
[docs]def lu_inplace(A, pivot="lapack"):
"""
In place LU decomposition.
Parameters
----------
A: af.Array
- a 2 dimensional arrayfire array on entry.
- Contains L in the lower triangle on exit.
- Contains U in the upper triangle on exit.
Returns
-------
P: af.Array
- Permutation array.
Note
----
This function is primarily used with `af.solve_lu` to reduce computations.
"""
P = Array()
is_pivot_lapack = False if (pivot == "full") else True
safe_call(backend.get().af_lu_inplace(ct.pointer(P.arr), A.arr, is_pivot_lapack))
return P
[docs]def qr(A):
"""
QR decomposition.
Parameters
----------
A: af.Array
A 2 dimensional arrayfire array.
Returns
-------
(Q,R,T): tuple of af.Arrays
- Q - Orthogonal matrix.
- R - Upper triangular matrix.
- T - Vector containing additional information to solve a least squares problem.
Note
----
The outputs of this funciton have the following properties.
>>> A = af.matmul(Q, R)
>>> I = af.matmulNT(Q, Q) # Identity matrix
"""
Q = Array()
R = Array()
T = Array()
safe_call(backend.get().af_lu(ct.pointer(Q.arr), ct.pointer(R.arr), ct.pointer(T.arr), A.arr))
return Q,R,T
[docs]def qr_inplace(A):
"""
In place QR decomposition.
Parameters
----------
A: af.Array
- a 2 dimensional arrayfire array on entry.
- Packed Q and R matrices on exit.
Returns
-------
T: af.Array
- Vector containing additional information to solve a least squares problem.
Note
----
This function is used to save space only when `R` is required.
"""
T = Array()
safe_call(backend.get().af_qr_inplace(ct.pointer(T.arr), A.arr))
return T
[docs]def cholesky(A, is_upper=True):
"""
Cholesky decomposition
Parameters
----------
A: af.Array
A 2 dimensional, symmetric, positive definite matrix.
is_upper: optional: bool. default: True
Specifies if output `R` is upper triangular (if True) or lower triangular (if False).
Returns
-------
(R,info): tuple of af.Array, int.
- R - triangular matrix.
- info - 0 if decomposition sucessful.
Note
----
The original matrix `A` can be reconstructed using the outputs in the following manner.
>>> A = af.matmulNT(R, R) #if R is upper triangular
"""
R = Array()
info = ct.c_int(0)
safe_call(backend.get().af_cholesky(ct.pointer(R.arr), ct.pointer(info), A.arr, is_upper))
return R, info.value
[docs]def cholesky_inplace(A, is_upper=True):
"""
In place Cholesky decomposition.
Parameters
----------
A: af.Array
- a 2 dimensional, symmetric, positive definite matrix.
- Trinangular matrix on exit.
is_upper: optional: bool. default: True.
Specifies if output `R` is upper triangular (if True) or lower triangular (if False).
Returns
-------
info : int.
0 if decomposition sucessful.
"""
info = ct.c_int(0)
safe_call(backend.get().af_cholesky_inplace(ct.pointer(info), A.arr, is_upper))
return info.value
[docs]def solve(A, B, options=MATPROP.NONE):
"""
Solve a system of linear equations.
Parameters
----------
A: af.Array
A 2 dimensional arrayfire array representing the coefficients of the system.
B: af.Array
A 1 or 2 dimensional arrayfire array representing the constants of the system.
options: optional: af.MATPROP. default: af.MATPROP.NONE.
- Additional options to speed up computations.
- Currently needs to be one of `af.MATPROP.NONE`, `af.MATPROP.LOWER`, `af.MATPROP.UPPER`.
Returns
-------
X: af.Array
A 1 or 2 dimensional arrayfire array representing the unknowns in the system.
"""
X = Array()
safe_call(backend.get().af_solve(ct.pointer(X.arr), A.arr, B.arr, options.value))
return X
[docs]def solve_lu(A, P, B, options=MATPROP.NONE):
"""
Solve a system of linear equations, using LU decomposition.
Parameters
----------
A: af.Array
- A 2 dimensional arrayfire array representing the coefficients of the system.
- This matrix should be decomposed previously using `lu_inplace(A)`.
P: af.Array
- Permutation array.
- This array is the output of an earlier call to `lu_inplace(A)`
B: af.Array
A 1 or 2 dimensional arrayfire array representing the constants of the system.
Returns
-------
X: af.Array
A 1 or 2 dimensional arrayfire array representing the unknowns in the system.
"""
X = Array()
safe_call(backend.get().af_solve_lu(ct.pointer(X.arr), A.arr, P.arr, B.arr, options.value))
return X
[docs]def inverse(A, options=MATPROP.NONE):
"""
Invert a matrix.
Parameters
----------
A: af.Array
- A 2 dimensional arrayfire array
options: optional: af.MATPROP. default: af.MATPROP.NONE.
- Additional options to speed up computations.
- Currently needs to be one of `af.MATPROP.NONE`.
Returns
-------
AI: af.Array
- A 2 dimensional array that is the inverse of `A`
Note
----
`A` needs to be a square matrix.
"""
AI = Array()
safe_call(backend.get().af_inverse(ct.pointer(AI.arr), A.arr, options.value))
return AI
[docs]def rank(A, tol=1E-5):
"""
Rank of a matrix.
Parameters
----------
A: af.Array
- A 2 dimensional arrayfire array
tol: optional: scalar. default: 1E-5.
- Tolerance for calculating rank
Returns
-------
r: int
- Rank of `A` within the given tolerance
"""
r = ct.c_uint(0)
safe_call(backend.get().af_rank(ct.pointer(r), A.arr, ct.c_double(tol)))
return r.value
[docs]def det(A):
"""
Determinant of a matrix.
Parameters
----------
A: af.Array
- A 2 dimensional arrayfire array
Returns
-------
res: scalar
- Determinant of the matrix.
"""
re = ct.c_double(0)
im = ct.c_double(0)
safe_call(backend.get().af_det(ct.pointer(re), ct.pointer(im), A.arr))
re = re.value
im = im.value
return re if (im == 0) else re + im * 1j
[docs]def norm(A, norm_type=NORM.EUCLID, p=1.0, q=1.0):
"""
Norm of an array or a matrix.
Parameters
----------
A: af.Array
- A 1 or 2 dimensional arrayfire array
norm_type: optional: af.NORM. default: af.NORM.EUCLID.
- Type of norm to be calculated.
p: scalar. default 1.0.
- Used only if `norm_type` is one of `af.NORM.VECTOR_P`, `af.NORM_MATRIX_L_PQ`
q: scalar. default 1.0.
- Used only if `norm_type` is `af.NORM_MATRIX_L_PQ`
Returns
-------
res: scalar
- norm of the input
"""
res = ct.c_double(0)
safe_call(backend.get().af_norm(ct.pointer(res), A.arr, norm_type.value,
ct.c_double(p), ct.c_double(q)))
return res.value
[docs]def svd(A):
"""
Singular Value Decomposition
Parameters
----------
A: af.Array
A 2 dimensional arrayfire array.
Returns
-------
(U,S,Vt): tuple of af.Arrays
- U - A unitary matrix
- S - An array containing the elements of diagonal matrix
- Vt - A unitary matrix
Note
----
- The original matrix `A` is preserved and additional storage space is required for decomposition.
- If the original matrix `A` need not be preserved, use `svd_inplace` instead.
- The original matrix `A` can be reconstructed using the outputs in the following manner.
>>> Smat = af.diag(S, 0, False)
>>> A_recon = af.matmul(af.matmul(U, Smat), Vt)
"""
U = Array()
S = Array()
Vt = Array()
safe_call(backend.get().af_svd(ct.pointer(U.arr), ct.pointer(S.arr), ct.pointer(Vt.arr), A.arr))
return U, S, Vt
[docs]def svd_inplace(A):
"""
Singular Value Decomposition
Parameters
----------
A: af.Array
A 2 dimensional arrayfire array.
Returns
-------
(U,S,Vt): tuple of af.Arrays
- U - A unitary matrix
- S - An array containing the elements of diagonal matrix
- Vt - A unitary matrix
Note
----
- The original matrix `A` is not preserved.
- If the original matrix `A` needs to be preserved, use `svd` instead.
- The original matrix `A` can be reconstructed using the outputs in the following manner.
>>> Smat = af.diag(S, 0, False)
>>> A_recon = af.matmul(af.matmul(U, Smat), Vt)
"""
U = Array()
S = Array()
Vt = Array()
safe_call(backend.get().af_svd_inplace(ct.pointer(U.arr), ct.pointer(S.arr), ct.pointer(Vt.arr),
A.arr))
return U, S, Vt
[docs]def is_lapack_available():
"""
Function to check if the arrayfire library was built with lapack support.
"""
res = ct.c_bool(False)
safe_call(backend.get().af_is_lapack_available(ct.pointer(res)))
return res.value
