In [1]:
import numpy as np
from gth_solve_jit import gth_solve, gth_solve_jit
In [2]:
A = [[0.4, 0.6], [0.2, 0.8]]
In [3]:
gth_solve(A)
Out[3]:
array([ 0.25, 0.75])
In [4]:
gth_solve_jit(A)
Out[4]:
array([ 0.25, 0.75])
In [5]:
B = [[0, 1, 0, 0], [1, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]] # Reducible matrix
In [6]:
gth_solve(B)
Out[6]:
array([ 0.5, 0.5, 0. , 0. ])
In [7]:
gth_solve_jit(B)
Out[7]:
array([ 0.5, 0.5, 0. , 0. ])
In [8]:
def random_probvec(k, m):
"""
Create probability vectors.
Parameters
----------
k : scalar(int)
Dimension of each probability vectors.
m : scalar(int)
Number of probability vectors.
Returns
-------
ndarray(float, ndim=2)
Array of shape (m, k) containing probability vectors as rows.
"""
x = np.empty((m, k+1))
r = np.random.rand(m, k-1)
r.sort(axis=-1)
x[:, 0], x[:, 1:k], x[:, k] = 0, r, 1
return np.diff(x, axis=-1)
In [9]:
random_probvec(2, 2)
Out[9]:
array([[ 0.56655772, 0.43344228],
[ 0.43229485, 0.56770515]])
In [10]:
sizes = [10, 50, 100, 1000]
rand_matrices = []
for n in sizes:
Q = random_probvec(n, n)
rand_matrices.append(Q)
In [11]:
for i, Q in enumerate(rand_matrices):
print 'rand_matrices[{0}] ({1} x {2})'.format(i, Q.shape[0], Q.shape[1])
%timeit gth_solve(Q)
%timeit gth_solve_jit(Q)
rand_matrices[0] (10 x 10)
1000 loops, best of 3: 234 µs per loop
1 loops, best of 3: 20 µs per loop
rand_matrices[1] (50 x 50)
1000 loops, best of 3: 1.44 ms per loop
10000 loops, best of 3: 97.4 µs per loop
rand_matrices[2] (100 x 100)
100 loops, best of 3: 4.24 ms per loop
1000 loops, best of 3: 633 µs per loop
rand_matrices[3] (1000 x 1000)
1 loops, best of 3: 1.64 s per loop
1 loops, best of 3: 617 ms per loop
In [12]:
import platform
print platform.platform()
Darwin-10.8.0-i386-64bit
In [13]:
import sys
print sys.version
2.7.8 (default, Aug 12 2014, 11:32:49)
[GCC 4.2.1 (Apple Inc. build 5666) (dot 3)]
In [14]:
print np.__version__
1.9.2
In [15]:
import numba
print numba.__version__
0.15.1
In [16]:
gth_solve_jit.inspect_types()
gth_solve_jit (array(float64, 2d, C, nonconst),)
--------------------------------------------------------------------------------
# File: gth_solve_jit.py
# --- LINE 97 ---
@jit
# --- LINE 98 ---
def gth_solve_jit(A):
# --- LINE 99 ---
"""
# --- LINE 100 ---
JIT-compiled version of `gth_solve`.
# --- LINE 101 ---
# --- LINE 102 ---
"""
# --- LINE 103 ---
#A1 = np.array(A, dtype=np.float64)
# --- LINE 104 ---
# label 0
# A.1 = A :: pyobject
# del A
# $0.1 = global(np: <module 'numpy' from '/usr/local/lib/python2.7/site-packages/numpy/__init__.pyc'>) :: pyobject
# $0.2 = getattr(attr=array, value=$0.1) :: pyobject
# del $0.1
# $0.5 = global(float: <type 'float'>) :: pyobject
# $0.6 = call $0.2(A.1, dtype=$0.5) :: pyobject
# del A.1
# del $0.5
# del $0.2
# A1 = $0.6 :: pyobject
# del $0.6
A1 = np.array(A, dtype=float)
# --- LINE 105 ---
# --- LINE 106 ---
# $0.7 = global(len: <built-in function len>) :: pyobject
# $0.9 = getattr(attr=shape, value=A1) :: pyobject
# $0.10 = call $0.7($0.9, ) :: pyobject
# del $0.9
# del $0.7
# $const0.11 = const(int, 2) :: pyobject
# $0.12 = $0.10 != $const0.11 :: pyobject
# del $const0.11
# del $0.10
# branch $0.12, 68, 42
# label 42
# $42.2 = getattr(attr=shape, value=A1) :: pyobject
# $const42.3 = const(int, 0) :: pyobject
# $42.4 = getitem(index=$const42.3, value=$42.2) :: pyobject
# del $const42.3
# del $42.2
# $42.6 = getattr(attr=shape, value=A1) :: pyobject
# $const42.7 = const(int, 1) :: pyobject
# $42.8 = getitem(index=$const42.7, value=$42.6) :: pyobject
# del $const42.7
# del $42.6
# $42.9 = $42.4 != $42.8 :: pyobject
# del $42.8
# del $42.4
# branch $42.9, 68, 77
if len(A1.shape) != 2 or A1.shape[0] != A1.shape[1]:
# --- LINE 107 ---
# raise ValueError('matrix must be square') # Not supported
# --- LINE 108 ---
# label 68
# del A1
# del $42.9
# del $0.12
# $68.1 = global(ValueError: <type 'exceptions.ValueError'>) :: pyobject
# del $68.1
# raise <type 'exceptions.ValueError'>
raise ValueError
# --- LINE 109 ---
# --- LINE 110 ---
# label 77
# del $42.9
# del $0.12
# $77.2 = getattr(attr=shape, value=A1) :: pyobject
# $const77.3 = const(int, 0) :: pyobject
# $77.4 = getitem(index=$const77.3, value=$77.2) :: pyobject
# del $const77.3
# del $77.2
# n = $77.4 :: pyobject
# del $77.4
n = A1.shape[0]
# --- LINE 111 ---
# --- LINE 112 ---
#x = np.zeros(n, dtype=np.float64)
# --- LINE 113 ---
# $77.5 = global(np: <module 'numpy' from '/usr/local/lib/python2.7/site-packages/numpy/__init__.pyc'>) :: pyobject
# $77.6 = getattr(attr=zeros, value=$77.5) :: pyobject
# del $77.5
# $77.8 = call $77.6(n, ) :: pyobject
# del $77.6
# x = $77.8 :: pyobject
# del $77.8
x = np.zeros(n)
# --- LINE 114 ---
# --- LINE 115 ---
# === Reduction === #
# --- LINE 116 ---
# $const77.9 = const(int, 1) :: pyobject
# SCALE_NONZERO = $const77.9 :: pyobject
# del $const77.9
SCALE_NONZERO = 1
# --- LINE 117 ---
# jump 114.1
# label 114.1
# $const114.1.1 = const(LiftedLoop, LiftedLoop(<function gth_solve_jit at 0x10915e1b8>)) :: pyobject
# $114.1.5 = call $const114.1.1(A1, SCALE_NONZERO, n, ) :: pyobject
# del $const114.1.1
# $114.1.8 = exhaust_iter(count=2, value=$114.1.5) :: pyobject
# del $114.1.5
# $114.1.6 = static_getitem(index=0, value=$114.1.8) :: pyobject
# $114.1.7 = static_getitem(index=1, value=$114.1.8) :: pyobject
# del $114.1.8
# SCALE_NONZERO = $114.1.6 :: pyobject
# del SCALE_NONZERO
# del $114.1.6
# n = $114.1.7 :: pyobject
# del $114.1.7
# jump 347
for k in range(n-1):
# --- LINE 118 ---
if SCALE_NONZERO == 1:
# --- LINE 119 ---
scale = np.sum(A1[k, k+1:n])
# --- LINE 120 ---
if scale <= 0:
# --- LINE 121 ---
# There is one (and only one) recurrent class contained in
# --- LINE 122 ---
# {0, ..., k};
# --- LINE 123 ---
# compute the solution associated with that recurrent class.
# --- LINE 124 ---
n = k+1
# --- LINE 125 ---
# break
# --- LINE 126 ---
# `break` statement is not supported for auto jitting
# --- LINE 127 ---
# see github.com/numba/numba/pull/819
# --- LINE 128 ---
# so need to use the flag `SCALE_NONZERO`
# --- LINE 129 ---
SCALE_NONZERO = 0
# --- LINE 130 ---
for i in range(k+1, n):
# --- LINE 131 ---
A1[i, k] /= scale
# --- LINE 132 ---
# --- LINE 133 ---
for j in range(k+1, n):
# --- LINE 134 ---
A1[i, j] += A1[i, k] * A1[k, j]
# --- LINE 135 ---
# --- LINE 136 ---
# === Backward substitution === #
# --- LINE 137 ---
# label 347
# $const347.1 = const(int, 1) :: pyobject
# $const347.4 = const(int, 1) :: pyobject
# $347.5 = n - $const347.4 :: pyobject
# del $const347.4
# x[$347.5] = $const347.1 :: pyobject
# del $const347.1
# del $347.5
x[n-1] = 1
# --- LINE 138 ---
# label 364.1
# $const364.1.1 = const(LiftedLoop, LiftedLoop(<function gth_solve_jit at 0x10915e1b8>)) :: pyobject
# $364.1.5 = call $const364.1.1(A1, x, n, ) :: pyobject
# del A1
# del $const364.1.1
# del $364.1.5
# jump 458
# jump 364.1
for k in range(n-2, -1, -1):
# --- LINE 139 ---
for i in range(k+1, n):
# --- LINE 140 ---
x[k] += x[i] * A1[i, k]
# --- LINE 141 ---
# --- LINE 142 ---
# === Normalization === #
# --- LINE 143 ---
# label 458
# $458.1 = global(np: <module 'numpy' from '/usr/local/lib/python2.7/site-packages/numpy/__init__.pyc'>) :: pyobject
# $458.2 = getattr(attr=sum, value=$458.1) :: pyobject
# del $458.1
# $458.4 = call $458.2(x, ) :: pyobject
# del $458.2
# norm = $458.4 :: pyobject
# del $458.4
norm = np.sum(x)
# --- LINE 144 ---
# label 476.1
# $const476.1.1 = const(LiftedLoop, LiftedLoop(<function gth_solve_jit at 0x10915e1b8>)) :: pyobject
# $476.1.5 = call $const476.1.1(x, norm, n, ) :: pyobject
# del norm
# del n
# del $const476.1.1
# del $476.1.5
# jump 512
# jump 476.1
for k in range(n):
# --- LINE 145 ---
x[k] /= norm
# --- LINE 146 ---
# --- LINE 147 ---
# label 512
# $512.2 = cast(value=x) :: pyobject
# del x
# return $512.2
return x
# The function contains lifted loops
# Loop at line 117
# Has 1 overloads
# File: gth_solve_jit.py
# --- LINE 97 ---
@jit
# --- LINE 98 ---
def gth_solve_jit(A):
# --- LINE 99 ---
"""
# --- LINE 100 ---
JIT-compiled version of `gth_solve`.
# --- LINE 101 ---
# --- LINE 102 ---
"""
# --- LINE 103 ---
#A1 = np.array(A, dtype=np.float64)
# --- LINE 104 ---
A1 = np.array(A, dtype=float)
# --- LINE 105 ---
# --- LINE 106 ---
if len(A1.shape) != 2 or A1.shape[0] != A1.shape[1]:
# --- LINE 107 ---
# raise ValueError('matrix must be square') # Not supported
# --- LINE 108 ---
raise ValueError
# --- LINE 109 ---
# --- LINE 110 ---
n = A1.shape[0]
# --- LINE 111 ---
# --- LINE 112 ---
#x = np.zeros(n, dtype=np.float64)
# --- LINE 113 ---
x = np.zeros(n)
# --- LINE 114 ---
# --- LINE 115 ---
# === Reduction === #
# --- LINE 116 ---
SCALE_NONZERO = 1
# --- LINE 117 ---
# label 128
# $128.2 = iternext(value=$phi128.1) :: pair<int64, bool>
# $128.3 = pair_first(value=$128.2) :: int64
# $128.4 = pair_second(value=$128.2) :: bool
# del $128.2
# $phi131.1 = $128.3 :: int64
# del $128.3
# branch $128.4, 131, 346
# label 131
# k = $phi131.1 :: int64
# del $phi131.1
# label 111
# A1.1 = A1 :: array(float64, 2d, C, nonconst)
# del A1
# SCALE_NONZERO.1 = SCALE_NONZERO :: int64
# del SCALE_NONZERO
# n.1 = n :: int64
# del n
# $111.1 = global(range: <built-in function range>) :: range
# $const111.3 = const(int, 1) :: int32
# $111.4 = n.1 - $const111.3 :: int64
# del $const111.3
# $111.5 = call $111.1($111.4, ) :: (int64,) -> range_state64
# del $111.4
# del $111.1
# $111.6 = getiter(value=$111.5) :: range_iter64
# del $111.5
# $phi128.1 = $111.6 :: range_iter64
# del $111.6
# jump 128
for k in range(n-1):
# --- LINE 118 ---
# $const131.3 = const(int, 1) :: int32
# $131.4 = SCALE_NONZERO.1 == $const131.3 :: bool
# del $const131.3
# branch $131.4, 146, 128
if SCALE_NONZERO == 1:
# --- LINE 119 ---
# label 146
# del $131.4
# $146.1 = global(np: <module 'numpy' from '/usr/local/lib/python2.7/site-packages/numpy/__init__.pyc'>) :: Module(<module 'numpy' from '/usr/local/lib/python2.7/site-packages/numpy/__init__.pyc'>)
# $146.2 = getattr(attr=sum, value=$146.1) :: Function(<unbound method Numpy_reduce_<function sum at 0x1027ad848>.sum>)
# del $146.1
# $k146.4 = k :: int64
# $const146.6 = const(int, 1) :: int32
# $146.7 = k + $const146.6 :: int64
# del $const146.6
# $146.9 = global(slice: <type 'slice'>) :: slice
# $146.10 = call $146.9($146.7, n.1, ) :: (int64, int64) -> slice3_type
# del $146.9
# del $146.7
# $146.11 = build_tuple(items=[Var($k146.4, gth_solve_jit.py (119)), Var($146.10, gth_solve_jit.py (119))]) :: (int64, slice3_type)
# del $k146.4
# del $146.10
# $146.12 = getitem(index=$146.11, value=A1.1) :: array(float64, 1d, A, nonconst)
# del $146.11
# $146.13 = call $146.2($146.12, ) :: (array(float64, 1d, A, nonconst),) -> float64
# del $146.2
# del $146.12
# scale = $146.13 :: float64
# del $146.13
scale = np.sum(A1[k, k+1:n])
# --- LINE 120 ---
# $const146.15 = const(int, 0) :: int32
# $146.16 = scale <= $const146.15 :: bool
# del $const146.15
# branch $146.16, 193, 212
if scale <= 0:
# --- LINE 121 ---
# There is one (and only one) recurrent class contained in
# --- LINE 122 ---
# {0, ..., k};
# --- LINE 123 ---
# compute the solution associated with that recurrent class.
# --- LINE 124 ---
# label 193
# $const193.2 = const(int, 1) :: int32
# $193.3 = k + $const193.2 :: int64
# del $const193.2
# n.1 = $193.3 :: int64
# del $193.3
n = k+1
# --- LINE 125 ---
# break
# --- LINE 126 ---
# `break` statement is not supported for auto jitting
# --- LINE 127 ---
# see github.com/numba/numba/pull/819
# --- LINE 128 ---
# so need to use the flag `SCALE_NONZERO`
# --- LINE 129 ---
# $const193.4 = const(int, 0) :: int32
# SCALE_NONZERO.1 = $const193.4 :: int64
# del $const193.4
# jump 212
SCALE_NONZERO = 0
# --- LINE 130 ---
# label 232
# $232.2 = iternext(value=$phi232.1) :: pair<int64, bool>
# $232.3 = pair_first(value=$232.2) :: int64
# $232.4 = pair_second(value=$232.2) :: bool
# del $232.2
# $phi235.1 = $232.3 :: int64
# del $232.3
# branch $232.4, 235, 339
# label 235
# i = $phi235.1 :: int64
# del $phi235.1
# label 212
# del $146.16
# $212.1 = global(range: <built-in function range>) :: range
# $const212.3 = const(int, 1) :: int32
# $212.4 = k + $const212.3 :: int64
# del $const212.3
# $212.6 = call $212.1($212.4, n.1, ) :: (int64, int64) -> range_state64
# del $212.4
# del $212.1
# $212.7 = getiter(value=$212.6) :: range_iter64
# del $212.6
# $phi232.1 = $212.7 :: range_iter64
# del $212.7
# jump 232
for i in range(k+1, n):
# --- LINE 131 ---
# $A1235.2 = A1.1 :: array(float64, 2d, C, nonconst)
# $235.5 = build_tuple(items=[Var(i, gth_solve_jit.py (130)), Var(k, gth_solve_jit.py (117))]) :: (int64 x 2)
# $235.8 = getitem(index=$235.5, value=A1.1) :: float64
# $235.10 = inplace_binop(rhs=scale, lhs=$235.8, fn=/?) :: float64
# del $235.8
# $A1235.2[$235.5] = $235.10 :: (array(float64, 2d, C, nonconst), (int64 x 2), float64) -> none
# del $A1235.2
# del $235.5
# del $235.10
A1[i, k] /= scale
# --- LINE 132 ---
# --- LINE 133 ---
# label 260
# $260.1 = global(range: <built-in function range>) :: range
# $const260.3 = const(int, 1) :: int32
# $260.4 = k + $const260.3 :: int64
# del $const260.3
# $260.6 = call $260.1($260.4, n.1, ) :: (int64, int64) -> range_state64
# del $260.4
# del $260.1
# $260.7 = getiter(value=$260.6) :: range_iter64
# del $260.6
# $phi280.1 = $260.7 :: range_iter64
# del $260.7
# jump 280
# jump 260
# label 280
# $280.2 = iternext(value=$phi280.1) :: pair<int64, bool>
# $280.3 = pair_first(value=$280.2) :: int64
# $280.4 = pair_second(value=$280.2) :: bool
# del $280.2
# $phi283.1 = $280.3 :: int64
# del $280.3
# branch $280.4, 283, 335
# label 283
# j = $phi283.1 :: int64
# del $phi283.1
for j in range(k+1, n):
# --- LINE 134 ---
# label 335
# del i
# del $phi283.1
# del $phi280.1
# del $280.4
# jump 336
# label 336
# jump 232
# label 347
# $347.3 = build_tuple(items=[Var(SCALE_NONZERO.1, gth_solve_jit.py (117)), Var(n.1, gth_solve_jit.py (117))]) :: (int64 x 2)
# del n.1
# del SCALE_NONZERO.1
# $347.4 = cast(value=$347.3) :: (int64 x 2)
# del $347.3
# return $347.4
# label 339
# del scale
# del k
# del $phi235.1
# del $phi232.1
# del $232.4
# jump 128
# label 346
# del k
# del A1.1
# del $phi131.1
# del $phi128.1
# del $131.4
# del $128.4
# jump 347
# $A1283.2 = A1.1 :: array(float64, 2d, C, nonconst)
# $283.5 = build_tuple(items=[Var(i, gth_solve_jit.py (130)), Var(j, gth_solve_jit.py (133))]) :: (int64 x 2)
# $283.8 = getitem(index=$283.5, value=A1.1) :: float64
# $283.12 = build_tuple(items=[Var(i, gth_solve_jit.py (130)), Var(k, gth_solve_jit.py (117))]) :: (int64 x 2)
# $283.13 = getitem(index=$283.12, value=A1.1) :: float64
# del $283.12
# $283.17 = build_tuple(items=[Var(k, gth_solve_jit.py (117)), Var(j, gth_solve_jit.py (133))]) :: (int64 x 2)
# del j
# $283.18 = getitem(index=$283.17, value=A1.1) :: float64
# del $283.17
# $283.19 = $283.13 * $283.18 :: float64
# del $283.18
# del $283.13
# $283.20 = inplace_binop(rhs=$283.19, lhs=$283.8, fn=+) :: float64
# del $283.8
# del $283.19
# $A1283.2[$283.5] = $283.20 :: (array(float64, 2d, C, nonconst), (int64 x 2), float64) -> none
# del $A1283.2
# del $283.5
# del $283.20
# jump 280
A1[i, j] += A1[i, k] * A1[k, j]
# --- LINE 135 ---
# --- LINE 136 ---
# === Backward substitution === #
# --- LINE 137 ---
x[n-1] = 1
# --- LINE 138 ---
for k in range(n-2, -1, -1):
# --- LINE 139 ---
for i in range(k+1, n):
# --- LINE 140 ---
x[k] += x[i] * A1[i, k]
# --- LINE 141 ---
# --- LINE 142 ---
# === Normalization === #
# --- LINE 143 ---
norm = np.sum(x)
# --- LINE 144 ---
for k in range(n):
# --- LINE 145 ---
x[k] /= norm
# --- LINE 146 ---
# --- LINE 147 ---
return x
# Loop at line 138
# Has 1 overloads
# File: gth_solve_jit.py
# --- LINE 97 ---
@jit
# --- LINE 98 ---
def gth_solve_jit(A):
# --- LINE 99 ---
"""
# --- LINE 100 ---
JIT-compiled version of `gth_solve`.
# --- LINE 101 ---
# --- LINE 102 ---
"""
# --- LINE 103 ---
#A1 = np.array(A, dtype=np.float64)
# --- LINE 104 ---
A1 = np.array(A, dtype=float)
# --- LINE 105 ---
# --- LINE 106 ---
if len(A1.shape) != 2 or A1.shape[0] != A1.shape[1]:
# --- LINE 107 ---
# raise ValueError('matrix must be square') # Not supported
# --- LINE 108 ---
raise ValueError
# --- LINE 109 ---
# --- LINE 110 ---
n = A1.shape[0]
# --- LINE 111 ---
# --- LINE 112 ---
#x = np.zeros(n, dtype=np.float64)
# --- LINE 113 ---
x = np.zeros(n)
# --- LINE 114 ---
# --- LINE 115 ---
# === Reduction === #
# --- LINE 116 ---
SCALE_NONZERO = 1
# --- LINE 117 ---
for k in range(n-1):
# --- LINE 118 ---
if SCALE_NONZERO == 1:
# --- LINE 119 ---
scale = np.sum(A1[k, k+1:n])
# --- LINE 120 ---
if scale <= 0:
# --- LINE 121 ---
# There is one (and only one) recurrent class contained in
# --- LINE 122 ---
# {0, ..., k};
# --- LINE 123 ---
# compute the solution associated with that recurrent class.
# --- LINE 124 ---
n = k+1
# --- LINE 125 ---
# break
# --- LINE 126 ---
# `break` statement is not supported for auto jitting
# --- LINE 127 ---
# see github.com/numba/numba/pull/819
# --- LINE 128 ---
# so need to use the flag `SCALE_NONZERO`
# --- LINE 129 ---
SCALE_NONZERO = 0
# --- LINE 130 ---
for i in range(k+1, n):
# --- LINE 131 ---
A1[i, k] /= scale
# --- LINE 132 ---
# --- LINE 133 ---
for j in range(k+1, n):
# --- LINE 134 ---
A1[i, j] += A1[i, k] * A1[k, j]
# --- LINE 135 ---
# --- LINE 136 ---
# === Backward substitution === #
# --- LINE 137 ---
x[n-1] = 1
# --- LINE 138 ---
# label 384
# $384.2 = iternext(value=$phi384.1) :: pair<int64, bool>
# $384.3 = pair_first(value=$384.2) :: int64
# $384.4 = pair_second(value=$384.2) :: bool
# del $384.2
# $phi387.1 = $384.3 :: int64
# del $384.3
# branch $384.4, 387, 457
# label 387
# k = $phi387.1 :: int64
# del $phi387.1
# label 361
# A1.1 = A1 :: array(float64, 2d, C, nonconst)
# del A1
# x.1 = x :: array(float64, 1d, C, nonconst)
# del x
# n.1 = n :: int64
# del n
# $361.1 = global(range: <built-in function range>) :: range
# $const361.3 = const(int, 2) :: int32
# $361.4 = n.1 - $const361.3 :: int64
# del $const361.3
# $const361.5 = const(int, -1) :: int32
# $const361.6 = const(int, -1) :: int32
# $361.7 = call $361.1($361.4, $const361.5, $const361.6, ) :: (int64, int64, int64) -> range_state64
# del $const361.6
# del $const361.5
# del $361.4
# del $361.1
# $361.8 = getiter(value=$361.7) :: range_iter64
# del $361.7
# $phi384.1 = $361.8 :: range_iter64
# del $361.8
# jump 384
for k in range(n-2, -1, -1):
# --- LINE 139 ---
# jump 390
# label 390
# $390.1 = global(range: <built-in function range>) :: range
# $const390.3 = const(int, 1) :: int32
# $390.4 = k + $const390.3 :: int64
# del $const390.3
# $390.6 = call $390.1($390.4, n.1, ) :: (int64, int64) -> range_state64
# del $390.4
# del $390.1
# $390.7 = getiter(value=$390.6) :: range_iter64
# del $390.6
# $phi410.1 = $390.7 :: range_iter64
# del $390.7
# jump 410
# label 410
# $410.2 = iternext(value=$phi410.1) :: pair<int64, bool>
# $410.3 = pair_first(value=$410.2) :: int64
# $410.4 = pair_second(value=$410.2) :: bool
# del $410.2
# $phi413.1 = $410.3 :: int64
# del $410.3
# branch $410.4, 413, 453
# label 413
# i = $phi413.1 :: int64
# del $phi413.1
for i in range(k+1, n):
# --- LINE 140 ---
# label 453
# del k
# del $phi413.1
# del $phi410.1
# del $410.4
# jump 454
# label 454
# jump 384
# label 458
# $const458.1 = const(NoneType, None) :: none
# $458.2 = cast(value=$const458.1) :: none
# del $const458.1
# return $458.2
# label 457
# del x.1
# del n.1
# del A1.1
# del $phi387.1
# del $phi384.1
# del $384.4
# jump 458
# $x413.2 = x.1 :: array(float64, 1d, C, nonconst)
# $k413.3 = k :: int64
# $413.6 = getitem(index=k, value=x.1) :: float64
# $413.9 = getitem(index=i, value=x.1) :: float64
# $413.13 = build_tuple(items=[Var(i, gth_solve_jit.py (139)), Var(k, gth_solve_jit.py (138))]) :: (int64 x 2)
# del i
# $413.14 = getitem(index=$413.13, value=A1.1) :: float64
# del $413.13
# $413.15 = $413.9 * $413.14 :: float64
# del $413.9
# del $413.14
# $413.16 = inplace_binop(rhs=$413.15, lhs=$413.6, fn=+) :: float64
# del $413.6
# del $413.15
# $x413.2[$k413.3] = $413.16 :: (array(float64, 1d, C, nonconst), int64, float64) -> none
# del $x413.2
# del $k413.3
# del $413.16
# jump 410
x[k] += x[i] * A1[i, k]
# --- LINE 141 ---
# --- LINE 142 ---
# === Normalization === #
# --- LINE 143 ---
norm = np.sum(x)
# --- LINE 144 ---
for k in range(n):
# --- LINE 145 ---
x[k] /= norm
# --- LINE 146 ---
# --- LINE 147 ---
return x
# Loop at line 144
# Has 1 overloads
# File: gth_solve_jit.py
# --- LINE 97 ---
@jit
# --- LINE 98 ---
def gth_solve_jit(A):
# --- LINE 99 ---
"""
# --- LINE 100 ---
JIT-compiled version of `gth_solve`.
# --- LINE 101 ---
# --- LINE 102 ---
"""
# --- LINE 103 ---
#A1 = np.array(A, dtype=np.float64)
# --- LINE 104 ---
A1 = np.array(A, dtype=float)
# --- LINE 105 ---
# --- LINE 106 ---
if len(A1.shape) != 2 or A1.shape[0] != A1.shape[1]:
# --- LINE 107 ---
# raise ValueError('matrix must be square') # Not supported
# --- LINE 108 ---
raise ValueError
# --- LINE 109 ---
# --- LINE 110 ---
n = A1.shape[0]
# --- LINE 111 ---
# --- LINE 112 ---
#x = np.zeros(n, dtype=np.float64)
# --- LINE 113 ---
x = np.zeros(n)
# --- LINE 114 ---
# --- LINE 115 ---
# === Reduction === #
# --- LINE 116 ---
SCALE_NONZERO = 1
# --- LINE 117 ---
for k in range(n-1):
# --- LINE 118 ---
if SCALE_NONZERO == 1:
# --- LINE 119 ---
scale = np.sum(A1[k, k+1:n])
# --- LINE 120 ---
if scale <= 0:
# --- LINE 121 ---
# There is one (and only one) recurrent class contained in
# --- LINE 122 ---
# {0, ..., k};
# --- LINE 123 ---
# compute the solution associated with that recurrent class.
# --- LINE 124 ---
n = k+1
# --- LINE 125 ---
# break
# --- LINE 126 ---
# `break` statement is not supported for auto jitting
# --- LINE 127 ---
# see github.com/numba/numba/pull/819
# --- LINE 128 ---
# so need to use the flag `SCALE_NONZERO`
# --- LINE 129 ---
SCALE_NONZERO = 0
# --- LINE 130 ---
for i in range(k+1, n):
# --- LINE 131 ---
A1[i, k] /= scale
# --- LINE 132 ---
# --- LINE 133 ---
for j in range(k+1, n):
# --- LINE 134 ---
A1[i, j] += A1[i, k] * A1[k, j]
# --- LINE 135 ---
# --- LINE 136 ---
# === Backward substitution === #
# --- LINE 137 ---
x[n-1] = 1
# --- LINE 138 ---
for k in range(n-2, -1, -1):
# --- LINE 139 ---
for i in range(k+1, n):
# --- LINE 140 ---
x[k] += x[i] * A1[i, k]
# --- LINE 141 ---
# --- LINE 142 ---
# === Normalization === #
# --- LINE 143 ---
norm = np.sum(x)
# --- LINE 144 ---
# label 473
# x.1 = x :: array(float64, 1d, C, nonconst)
# del x
# norm.1 = norm :: float64
# del norm
# n.1 = n :: int64
# del n
# $473.1 = global(range: <built-in function range>) :: range
# $473.3 = call $473.1(n.1, ) :: (int64,) -> range_state64
# del n.1
# del $473.1
# $473.4 = getiter(value=$473.3) :: range_iter64
# del $473.3
# $phi486.1 = $473.4 :: range_iter64
# del $473.4
# jump 486
# label 486
# $486.2 = iternext(value=$phi486.1) :: pair<int64, bool>
# $486.3 = pair_first(value=$486.2) :: int64
# $486.4 = pair_second(value=$486.2) :: bool
# del $486.2
# $phi489.1 = $486.3 :: int64
# del $486.3
# branch $486.4, 489, 511
# label 489
# k = $phi489.1 :: int64
# del $phi489.1
for k in range(n):
# --- LINE 145 ---
# label 512
# $const512.1 = const(NoneType, None) :: none
# $512.2 = cast(value=$const512.1) :: none
# del $const512.1
# return $512.2
# label 511
# del x.1
# del norm.1
# del $phi489.1
# del $phi486.1
# del $486.4
# jump 512
# $x489.2 = x.1 :: array(float64, 1d, C, nonconst)
# $k489.3 = k :: int64
# $489.6 = getitem(index=k, value=x.1) :: float64
# del k
# $489.8 = inplace_binop(rhs=norm.1, lhs=$489.6, fn=/?) :: float64
# del $489.6
# $x489.2[$k489.3] = $489.8 :: (array(float64, 1d, C, nonconst), int64, float64) -> none
# del $x489.2
# del $k489.3
# del $489.8
# jump 486
x[k] /= norm
# --- LINE 146 ---
# --- LINE 147 ---
return x
================================================================================
gth_solve_jit (pyobject,)
--------------------------------------------------------------------------------
# File: gth_solve_jit.py
# --- LINE 97 ---
@jit
# --- LINE 98 ---
def gth_solve_jit(A):
# --- LINE 99 ---
"""
# --- LINE 100 ---
JIT-compiled version of `gth_solve`.
# --- LINE 101 ---
# --- LINE 102 ---
"""
# --- LINE 103 ---
#A1 = np.array(A, dtype=np.float64)
# --- LINE 104 ---
# label 0
# A.1 = A :: pyobject
# del A
# $0.1 = global(np: <module 'numpy' from '/usr/local/lib/python2.7/site-packages/numpy/__init__.pyc'>) :: pyobject
# $0.2 = getattr(attr=array, value=$0.1) :: pyobject
# del $0.1
# $0.5 = global(float: <type 'float'>) :: pyobject
# $0.6 = call $0.2(A.1, dtype=$0.5) :: pyobject
# del A.1
# del $0.5
# del $0.2
# A1 = $0.6 :: pyobject
# del $0.6
A1 = np.array(A, dtype=float)
# --- LINE 105 ---
# --- LINE 106 ---
# $0.7 = global(len: <built-in function len>) :: pyobject
# $0.9 = getattr(attr=shape, value=A1) :: pyobject
# $0.10 = call $0.7($0.9, ) :: pyobject
# del $0.9
# del $0.7
# $const0.11 = const(int, 2) :: pyobject
# $0.12 = $0.10 != $const0.11 :: pyobject
# del $const0.11
# del $0.10
# branch $0.12, 68, 42
# label 42
# $42.2 = getattr(attr=shape, value=A1) :: pyobject
# $const42.3 = const(int, 0) :: pyobject
# $42.4 = getitem(index=$const42.3, value=$42.2) :: pyobject
# del $const42.3
# del $42.2
# $42.6 = getattr(attr=shape, value=A1) :: pyobject
# $const42.7 = const(int, 1) :: pyobject
# $42.8 = getitem(index=$const42.7, value=$42.6) :: pyobject
# del $const42.7
# del $42.6
# $42.9 = $42.4 != $42.8 :: pyobject
# del $42.8
# del $42.4
# branch $42.9, 68, 77
if len(A1.shape) != 2 or A1.shape[0] != A1.shape[1]:
# --- LINE 107 ---
# raise ValueError('matrix must be square') # Not supported
# --- LINE 108 ---
# label 68
# del A1
# del $42.9
# del $0.12
# $68.1 = global(ValueError: <type 'exceptions.ValueError'>) :: pyobject
# del $68.1
# raise <type 'exceptions.ValueError'>
raise ValueError
# --- LINE 109 ---
# --- LINE 110 ---
# label 77
# del $42.9
# del $0.12
# $77.2 = getattr(attr=shape, value=A1) :: pyobject
# $const77.3 = const(int, 0) :: pyobject
# $77.4 = getitem(index=$const77.3, value=$77.2) :: pyobject
# del $const77.3
# del $77.2
# n = $77.4 :: pyobject
# del $77.4
n = A1.shape[0]
# --- LINE 111 ---
# --- LINE 112 ---
#x = np.zeros(n, dtype=np.float64)
# --- LINE 113 ---
# $77.5 = global(np: <module 'numpy' from '/usr/local/lib/python2.7/site-packages/numpy/__init__.pyc'>) :: pyobject
# $77.6 = getattr(attr=zeros, value=$77.5) :: pyobject
# del $77.5
# $77.8 = call $77.6(n, ) :: pyobject
# del $77.6
# x = $77.8 :: pyobject
# del $77.8
x = np.zeros(n)
# --- LINE 114 ---
# --- LINE 115 ---
# === Reduction === #
# --- LINE 116 ---
# $const77.9 = const(int, 1) :: pyobject
# SCALE_NONZERO = $const77.9 :: pyobject
# del $const77.9
SCALE_NONZERO = 1
# --- LINE 117 ---
# jump 114.1
# label 114.1
# $const114.1.1 = const(LiftedLoop, LiftedLoop(<function gth_solve_jit at 0x10915e1b8>)) :: pyobject
# $114.1.5 = call $const114.1.1(A1, SCALE_NONZERO, n, ) :: pyobject
# del $const114.1.1
# $114.1.8 = exhaust_iter(count=2, value=$114.1.5) :: pyobject
# del $114.1.5
# $114.1.6 = static_getitem(index=0, value=$114.1.8) :: pyobject
# $114.1.7 = static_getitem(index=1, value=$114.1.8) :: pyobject
# del $114.1.8
# SCALE_NONZERO = $114.1.6 :: pyobject
# del SCALE_NONZERO
# del $114.1.6
# n = $114.1.7 :: pyobject
# del $114.1.7
# jump 347
for k in range(n-1):
# --- LINE 118 ---
if SCALE_NONZERO == 1:
# --- LINE 119 ---
scale = np.sum(A1[k, k+1:n])
# --- LINE 120 ---
if scale <= 0:
# --- LINE 121 ---
# There is one (and only one) recurrent class contained in
# --- LINE 122 ---
# {0, ..., k};
# --- LINE 123 ---
# compute the solution associated with that recurrent class.
# --- LINE 124 ---
n = k+1
# --- LINE 125 ---
# break
# --- LINE 126 ---
# `break` statement is not supported for auto jitting
# --- LINE 127 ---
# see github.com/numba/numba/pull/819
# --- LINE 128 ---
# so need to use the flag `SCALE_NONZERO`
# --- LINE 129 ---
SCALE_NONZERO = 0
# --- LINE 130 ---
for i in range(k+1, n):
# --- LINE 131 ---
A1[i, k] /= scale
# --- LINE 132 ---
# --- LINE 133 ---
for j in range(k+1, n):
# --- LINE 134 ---
A1[i, j] += A1[i, k] * A1[k, j]
# --- LINE 135 ---
# --- LINE 136 ---
# === Backward substitution === #
# --- LINE 137 ---
# label 347
# $const347.1 = const(int, 1) :: pyobject
# $const347.4 = const(int, 1) :: pyobject
# $347.5 = n - $const347.4 :: pyobject
# del $const347.4
# x[$347.5] = $const347.1 :: pyobject
# del $const347.1
# del $347.5
x[n-1] = 1
# --- LINE 138 ---
# label 364.1
# $const364.1.1 = const(LiftedLoop, LiftedLoop(<function gth_solve_jit at 0x10915e1b8>)) :: pyobject
# $364.1.5 = call $const364.1.1(A1, x, n, ) :: pyobject
# del A1
# del $const364.1.1
# del $364.1.5
# jump 458
# jump 364.1
for k in range(n-2, -1, -1):
# --- LINE 139 ---
for i in range(k+1, n):
# --- LINE 140 ---
x[k] += x[i] * A1[i, k]
# --- LINE 141 ---
# --- LINE 142 ---
# === Normalization === #
# --- LINE 143 ---
# label 458
# $458.1 = global(np: <module 'numpy' from '/usr/local/lib/python2.7/site-packages/numpy/__init__.pyc'>) :: pyobject
# $458.2 = getattr(attr=sum, value=$458.1) :: pyobject
# del $458.1
# $458.4 = call $458.2(x, ) :: pyobject
# del $458.2
# norm = $458.4 :: pyobject
# del $458.4
norm = np.sum(x)
# --- LINE 144 ---
# label 476.1
# $const476.1.1 = const(LiftedLoop, LiftedLoop(<function gth_solve_jit at 0x10915e1b8>)) :: pyobject
# $476.1.5 = call $const476.1.1(x, norm, n, ) :: pyobject
# del norm
# del n
# del $const476.1.1
# del $476.1.5
# jump 512
# jump 476.1
for k in range(n):
# --- LINE 145 ---
x[k] /= norm
# --- LINE 146 ---
# --- LINE 147 ---
# label 512
# $512.2 = cast(value=x) :: pyobject
# del x
# return $512.2
return x
# The function contains lifted loops
# Loop at line 117
# Has 1 overloads
# File: gth_solve_jit.py
# --- LINE 97 ---
@jit
# --- LINE 98 ---
def gth_solve_jit(A):
# --- LINE 99 ---
"""
# --- LINE 100 ---
JIT-compiled version of `gth_solve`.
# --- LINE 101 ---
# --- LINE 102 ---
"""
# --- LINE 103 ---
#A1 = np.array(A, dtype=np.float64)
# --- LINE 104 ---
A1 = np.array(A, dtype=float)
# --- LINE 105 ---
# --- LINE 106 ---
if len(A1.shape) != 2 or A1.shape[0] != A1.shape[1]:
# --- LINE 107 ---
# raise ValueError('matrix must be square') # Not supported
# --- LINE 108 ---
raise ValueError
# --- LINE 109 ---
# --- LINE 110 ---
n = A1.shape[0]
# --- LINE 111 ---
# --- LINE 112 ---
#x = np.zeros(n, dtype=np.float64)
# --- LINE 113 ---
x = np.zeros(n)
# --- LINE 114 ---
# --- LINE 115 ---
# === Reduction === #
# --- LINE 116 ---
SCALE_NONZERO = 1
# --- LINE 117 ---
# label 128
# $128.2 = iternext(value=$phi128.1) :: pair<int64, bool>
# $128.3 = pair_first(value=$128.2) :: int64
# $128.4 = pair_second(value=$128.2) :: bool
# del $128.2
# $phi131.1 = $128.3 :: int64
# del $128.3
# branch $128.4, 131, 346
# label 131
# k = $phi131.1 :: int64
# del $phi131.1
# label 111
# A1.1 = A1 :: array(float64, 2d, C, nonconst)
# del A1
# SCALE_NONZERO.1 = SCALE_NONZERO :: int64
# del SCALE_NONZERO
# n.1 = n :: int64
# del n
# $111.1 = global(range: <built-in function range>) :: range
# $const111.3 = const(int, 1) :: int32
# $111.4 = n.1 - $const111.3 :: int64
# del $const111.3
# $111.5 = call $111.1($111.4, ) :: (int64,) -> range_state64
# del $111.4
# del $111.1
# $111.6 = getiter(value=$111.5) :: range_iter64
# del $111.5
# $phi128.1 = $111.6 :: range_iter64
# del $111.6
# jump 128
for k in range(n-1):
# --- LINE 118 ---
# $const131.3 = const(int, 1) :: int32
# $131.4 = SCALE_NONZERO.1 == $const131.3 :: bool
# del $const131.3
# branch $131.4, 146, 128
if SCALE_NONZERO == 1:
# --- LINE 119 ---
# label 146
# del $131.4
# $146.1 = global(np: <module 'numpy' from '/usr/local/lib/python2.7/site-packages/numpy/__init__.pyc'>) :: Module(<module 'numpy' from '/usr/local/lib/python2.7/site-packages/numpy/__init__.pyc'>)
# $146.2 = getattr(attr=sum, value=$146.1) :: Function(<unbound method Numpy_reduce_<function sum at 0x1027ad848>.sum>)
# del $146.1
# $k146.4 = k :: int64
# $const146.6 = const(int, 1) :: int32
# $146.7 = k + $const146.6 :: int64
# del $const146.6
# $146.9 = global(slice: <type 'slice'>) :: slice
# $146.10 = call $146.9($146.7, n.1, ) :: (int64, int64) -> slice3_type
# del $146.9
# del $146.7
# $146.11 = build_tuple(items=[Var($k146.4, gth_solve_jit.py (119)), Var($146.10, gth_solve_jit.py (119))]) :: (int64, slice3_type)
# del $k146.4
# del $146.10
# $146.12 = getitem(index=$146.11, value=A1.1) :: array(float64, 1d, A, nonconst)
# del $146.11
# $146.13 = call $146.2($146.12, ) :: (array(float64, 1d, A, nonconst),) -> float64
# del $146.2
# del $146.12
# scale = $146.13 :: float64
# del $146.13
scale = np.sum(A1[k, k+1:n])
# --- LINE 120 ---
# $const146.15 = const(int, 0) :: int32
# $146.16 = scale <= $const146.15 :: bool
# del $const146.15
# branch $146.16, 193, 212
if scale <= 0:
# --- LINE 121 ---
# There is one (and only one) recurrent class contained in
# --- LINE 122 ---
# {0, ..., k};
# --- LINE 123 ---
# compute the solution associated with that recurrent class.
# --- LINE 124 ---
# label 193
# $const193.2 = const(int, 1) :: int32
# $193.3 = k + $const193.2 :: int64
# del $const193.2
# n.1 = $193.3 :: int64
# del $193.3
n = k+1
# --- LINE 125 ---
# break
# --- LINE 126 ---
# `break` statement is not supported for auto jitting
# --- LINE 127 ---
# see github.com/numba/numba/pull/819
# --- LINE 128 ---
# so need to use the flag `SCALE_NONZERO`
# --- LINE 129 ---
# $const193.4 = const(int, 0) :: int32
# SCALE_NONZERO.1 = $const193.4 :: int64
# del $const193.4
# jump 212
SCALE_NONZERO = 0
# --- LINE 130 ---
# label 232
# $232.2 = iternext(value=$phi232.1) :: pair<int64, bool>
# $232.3 = pair_first(value=$232.2) :: int64
# $232.4 = pair_second(value=$232.2) :: bool
# del $232.2
# $phi235.1 = $232.3 :: int64
# del $232.3
# branch $232.4, 235, 339
# label 235
# i = $phi235.1 :: int64
# del $phi235.1
# label 212
# del $146.16
# $212.1 = global(range: <built-in function range>) :: range
# $const212.3 = const(int, 1) :: int32
# $212.4 = k + $const212.3 :: int64
# del $const212.3
# $212.6 = call $212.1($212.4, n.1, ) :: (int64, int64) -> range_state64
# del $212.4
# del $212.1
# $212.7 = getiter(value=$212.6) :: range_iter64
# del $212.6
# $phi232.1 = $212.7 :: range_iter64
# del $212.7
# jump 232
for i in range(k+1, n):
# --- LINE 131 ---
# $A1235.2 = A1.1 :: array(float64, 2d, C, nonconst)
# $235.5 = build_tuple(items=[Var(i, gth_solve_jit.py (130)), Var(k, gth_solve_jit.py (117))]) :: (int64 x 2)
# $235.8 = getitem(index=$235.5, value=A1.1) :: float64
# $235.10 = inplace_binop(rhs=scale, lhs=$235.8, fn=/?) :: float64
# del $235.8
# $A1235.2[$235.5] = $235.10 :: (array(float64, 2d, C, nonconst), (int64 x 2), float64) -> none
# del $A1235.2
# del $235.5
# del $235.10
A1[i, k] /= scale
# --- LINE 132 ---
# --- LINE 133 ---
# label 260
# $260.1 = global(range: <built-in function range>) :: range
# $const260.3 = const(int, 1) :: int32
# $260.4 = k + $const260.3 :: int64
# del $const260.3
# $260.6 = call $260.1($260.4, n.1, ) :: (int64, int64) -> range_state64
# del $260.4
# del $260.1
# $260.7 = getiter(value=$260.6) :: range_iter64
# del $260.6
# $phi280.1 = $260.7 :: range_iter64
# del $260.7
# jump 280
# jump 260
# label 280
# $280.2 = iternext(value=$phi280.1) :: pair<int64, bool>
# $280.3 = pair_first(value=$280.2) :: int64
# $280.4 = pair_second(value=$280.2) :: bool
# del $280.2
# $phi283.1 = $280.3 :: int64
# del $280.3
# branch $280.4, 283, 335
# label 283
# j = $phi283.1 :: int64
# del $phi283.1
for j in range(k+1, n):
# --- LINE 134 ---
# label 335
# del i
# del $phi283.1
# del $phi280.1
# del $280.4
# jump 336
# label 336
# jump 232
# label 347
# $347.3 = build_tuple(items=[Var(SCALE_NONZERO.1, gth_solve_jit.py (117)), Var(n.1, gth_solve_jit.py (117))]) :: (int64 x 2)
# del n.1
# del SCALE_NONZERO.1
# $347.4 = cast(value=$347.3) :: (int64 x 2)
# del $347.3
# return $347.4
# label 339
# del scale
# del k
# del $phi235.1
# del $phi232.1
# del $232.4
# jump 128
# label 346
# del k
# del A1.1
# del $phi131.1
# del $phi128.1
# del $131.4
# del $128.4
# jump 347
# $A1283.2 = A1.1 :: array(float64, 2d, C, nonconst)
# $283.5 = build_tuple(items=[Var(i, gth_solve_jit.py (130)), Var(j, gth_solve_jit.py (133))]) :: (int64 x 2)
# $283.8 = getitem(index=$283.5, value=A1.1) :: float64
# $283.12 = build_tuple(items=[Var(i, gth_solve_jit.py (130)), Var(k, gth_solve_jit.py (117))]) :: (int64 x 2)
# $283.13 = getitem(index=$283.12, value=A1.1) :: float64
# del $283.12
# $283.17 = build_tuple(items=[Var(k, gth_solve_jit.py (117)), Var(j, gth_solve_jit.py (133))]) :: (int64 x 2)
# del j
# $283.18 = getitem(index=$283.17, value=A1.1) :: float64
# del $283.17
# $283.19 = $283.13 * $283.18 :: float64
# del $283.18
# del $283.13
# $283.20 = inplace_binop(rhs=$283.19, lhs=$283.8, fn=+) :: float64
# del $283.8
# del $283.19
# $A1283.2[$283.5] = $283.20 :: (array(float64, 2d, C, nonconst), (int64 x 2), float64) -> none
# del $A1283.2
# del $283.5
# del $283.20
# jump 280
A1[i, j] += A1[i, k] * A1[k, j]
# --- LINE 135 ---
# --- LINE 136 ---
# === Backward substitution === #
# --- LINE 137 ---
x[n-1] = 1
# --- LINE 138 ---
for k in range(n-2, -1, -1):
# --- LINE 139 ---
for i in range(k+1, n):
# --- LINE 140 ---
x[k] += x[i] * A1[i, k]
# --- LINE 141 ---
# --- LINE 142 ---
# === Normalization === #
# --- LINE 143 ---
norm = np.sum(x)
# --- LINE 144 ---
for k in range(n):
# --- LINE 145 ---
x[k] /= norm
# --- LINE 146 ---
# --- LINE 147 ---
return x
# Loop at line 138
# Has 1 overloads
# File: gth_solve_jit.py
# --- LINE 97 ---
@jit
# --- LINE 98 ---
def gth_solve_jit(A):
# --- LINE 99 ---
"""
# --- LINE 100 ---
JIT-compiled version of `gth_solve`.
# --- LINE 101 ---
# --- LINE 102 ---
"""
# --- LINE 103 ---
#A1 = np.array(A, dtype=np.float64)
# --- LINE 104 ---
A1 = np.array(A, dtype=float)
# --- LINE 105 ---
# --- LINE 106 ---
if len(A1.shape) != 2 or A1.shape[0] != A1.shape[1]:
# --- LINE 107 ---
# raise ValueError('matrix must be square') # Not supported
# --- LINE 108 ---
raise ValueError
# --- LINE 109 ---
# --- LINE 110 ---
n = A1.shape[0]
# --- LINE 111 ---
# --- LINE 112 ---
#x = np.zeros(n, dtype=np.float64)
# --- LINE 113 ---
x = np.zeros(n)
# --- LINE 114 ---
# --- LINE 115 ---
# === Reduction === #
# --- LINE 116 ---
SCALE_NONZERO = 1
# --- LINE 117 ---
for k in range(n-1):
# --- LINE 118 ---
if SCALE_NONZERO == 1:
# --- LINE 119 ---
scale = np.sum(A1[k, k+1:n])
# --- LINE 120 ---
if scale <= 0:
# --- LINE 121 ---
# There is one (and only one) recurrent class contained in
# --- LINE 122 ---
# {0, ..., k};
# --- LINE 123 ---
# compute the solution associated with that recurrent class.
# --- LINE 124 ---
n = k+1
# --- LINE 125 ---
# break
# --- LINE 126 ---
# `break` statement is not supported for auto jitting
# --- LINE 127 ---
# see github.com/numba/numba/pull/819
# --- LINE 128 ---
# so need to use the flag `SCALE_NONZERO`
# --- LINE 129 ---
SCALE_NONZERO = 0
# --- LINE 130 ---
for i in range(k+1, n):
# --- LINE 131 ---
A1[i, k] /= scale
# --- LINE 132 ---
# --- LINE 133 ---
for j in range(k+1, n):
# --- LINE 134 ---
A1[i, j] += A1[i, k] * A1[k, j]
# --- LINE 135 ---
# --- LINE 136 ---
# === Backward substitution === #
# --- LINE 137 ---
x[n-1] = 1
# --- LINE 138 ---
# label 384
# $384.2 = iternext(value=$phi384.1) :: pair<int64, bool>
# $384.3 = pair_first(value=$384.2) :: int64
# $384.4 = pair_second(value=$384.2) :: bool
# del $384.2
# $phi387.1 = $384.3 :: int64
# del $384.3
# branch $384.4, 387, 457
# label 387
# k = $phi387.1 :: int64
# del $phi387.1
# label 361
# A1.1 = A1 :: array(float64, 2d, C, nonconst)
# del A1
# x.1 = x :: array(float64, 1d, C, nonconst)
# del x
# n.1 = n :: int64
# del n
# $361.1 = global(range: <built-in function range>) :: range
# $const361.3 = const(int, 2) :: int32
# $361.4 = n.1 - $const361.3 :: int64
# del $const361.3
# $const361.5 = const(int, -1) :: int32
# $const361.6 = const(int, -1) :: int32
# $361.7 = call $361.1($361.4, $const361.5, $const361.6, ) :: (int64, int64, int64) -> range_state64
# del $const361.6
# del $const361.5
# del $361.4
# del $361.1
# $361.8 = getiter(value=$361.7) :: range_iter64
# del $361.7
# $phi384.1 = $361.8 :: range_iter64
# del $361.8
# jump 384
for k in range(n-2, -1, -1):
# --- LINE 139 ---
# jump 390
# label 390
# $390.1 = global(range: <built-in function range>) :: range
# $const390.3 = const(int, 1) :: int32
# $390.4 = k + $const390.3 :: int64
# del $const390.3
# $390.6 = call $390.1($390.4, n.1, ) :: (int64, int64) -> range_state64
# del $390.4
# del $390.1
# $390.7 = getiter(value=$390.6) :: range_iter64
# del $390.6
# $phi410.1 = $390.7 :: range_iter64
# del $390.7
# jump 410
# label 410
# $410.2 = iternext(value=$phi410.1) :: pair<int64, bool>
# $410.3 = pair_first(value=$410.2) :: int64
# $410.4 = pair_second(value=$410.2) :: bool
# del $410.2
# $phi413.1 = $410.3 :: int64
# del $410.3
# branch $410.4, 413, 453
# label 413
# i = $phi413.1 :: int64
# del $phi413.1
for i in range(k+1, n):
# --- LINE 140 ---
# label 453
# del k
# del $phi413.1
# del $phi410.1
# del $410.4
# jump 454
# label 454
# jump 384
# label 458
# $const458.1 = const(NoneType, None) :: none
# $458.2 = cast(value=$const458.1) :: none
# del $const458.1
# return $458.2
# label 457
# del x.1
# del n.1
# del A1.1
# del $phi387.1
# del $phi384.1
# del $384.4
# jump 458
# $x413.2 = x.1 :: array(float64, 1d, C, nonconst)
# $k413.3 = k :: int64
# $413.6 = getitem(index=k, value=x.1) :: float64
# $413.9 = getitem(index=i, value=x.1) :: float64
# $413.13 = build_tuple(items=[Var(i, gth_solve_jit.py (139)), Var(k, gth_solve_jit.py (138))]) :: (int64 x 2)
# del i
# $413.14 = getitem(index=$413.13, value=A1.1) :: float64
# del $413.13
# $413.15 = $413.9 * $413.14 :: float64
# del $413.9
# del $413.14
# $413.16 = inplace_binop(rhs=$413.15, lhs=$413.6, fn=+) :: float64
# del $413.6
# del $413.15
# $x413.2[$k413.3] = $413.16 :: (array(float64, 1d, C, nonconst), int64, float64) -> none
# del $x413.2
# del $k413.3
# del $413.16
# jump 410
x[k] += x[i] * A1[i, k]
# --- LINE 141 ---
# --- LINE 142 ---
# === Normalization === #
# --- LINE 143 ---
norm = np.sum(x)
# --- LINE 144 ---
for k in range(n):
# --- LINE 145 ---
x[k] /= norm
# --- LINE 146 ---
# --- LINE 147 ---
return x
# Loop at line 144
# Has 1 overloads
# File: gth_solve_jit.py
# --- LINE 97 ---
@jit
# --- LINE 98 ---
def gth_solve_jit(A):
# --- LINE 99 ---
"""
# --- LINE 100 ---
JIT-compiled version of `gth_solve`.
# --- LINE 101 ---
# --- LINE 102 ---
"""
# --- LINE 103 ---
#A1 = np.array(A, dtype=np.float64)
# --- LINE 104 ---
A1 = np.array(A, dtype=float)
# --- LINE 105 ---
# --- LINE 106 ---
if len(A1.shape) != 2 or A1.shape[0] != A1.shape[1]:
# --- LINE 107 ---
# raise ValueError('matrix must be square') # Not supported
# --- LINE 108 ---
raise ValueError
# --- LINE 109 ---
# --- LINE 110 ---
n = A1.shape[0]
# --- LINE 111 ---
# --- LINE 112 ---
#x = np.zeros(n, dtype=np.float64)
# --- LINE 113 ---
x = np.zeros(n)
# --- LINE 114 ---
# --- LINE 115 ---
# === Reduction === #
# --- LINE 116 ---
SCALE_NONZERO = 1
# --- LINE 117 ---
for k in range(n-1):
# --- LINE 118 ---
if SCALE_NONZERO == 1:
# --- LINE 119 ---
scale = np.sum(A1[k, k+1:n])
# --- LINE 120 ---
if scale <= 0:
# --- LINE 121 ---
# There is one (and only one) recurrent class contained in
# --- LINE 122 ---
# {0, ..., k};
# --- LINE 123 ---
# compute the solution associated with that recurrent class.
# --- LINE 124 ---
n = k+1
# --- LINE 125 ---
# break
# --- LINE 126 ---
# `break` statement is not supported for auto jitting
# --- LINE 127 ---
# see github.com/numba/numba/pull/819
# --- LINE 128 ---
# so need to use the flag `SCALE_NONZERO`
# --- LINE 129 ---
SCALE_NONZERO = 0
# --- LINE 130 ---
for i in range(k+1, n):
# --- LINE 131 ---
A1[i, k] /= scale
# --- LINE 132 ---
# --- LINE 133 ---
for j in range(k+1, n):
# --- LINE 134 ---
A1[i, j] += A1[i, k] * A1[k, j]
# --- LINE 135 ---
# --- LINE 136 ---
# === Backward substitution === #
# --- LINE 137 ---
x[n-1] = 1
# --- LINE 138 ---
for k in range(n-2, -1, -1):
# --- LINE 139 ---
for i in range(k+1, n):
# --- LINE 140 ---
x[k] += x[i] * A1[i, k]
# --- LINE 141 ---
# --- LINE 142 ---
# === Normalization === #
# --- LINE 143 ---
norm = np.sum(x)
# --- LINE 144 ---
# label 473
# x.1 = x :: array(float64, 1d, C, nonconst)
# del x
# norm.1 = norm :: float64
# del norm
# n.1 = n :: int64
# del n
# $473.1 = global(range: <built-in function range>) :: range
# $473.3 = call $473.1(n.1, ) :: (int64,) -> range_state64
# del n.1
# del $473.1
# $473.4 = getiter(value=$473.3) :: range_iter64
# del $473.3
# $phi486.1 = $473.4 :: range_iter64
# del $473.4
# jump 486
# label 486
# $486.2 = iternext(value=$phi486.1) :: pair<int64, bool>
# $486.3 = pair_first(value=$486.2) :: int64
# $486.4 = pair_second(value=$486.2) :: bool
# del $486.2
# $phi489.1 = $486.3 :: int64
# del $486.3
# branch $486.4, 489, 511
# label 489
# k = $phi489.1 :: int64
# del $phi489.1
for k in range(n):
# --- LINE 145 ---
# label 512
# $const512.1 = const(NoneType, None) :: none
# $512.2 = cast(value=$const512.1) :: none
# del $const512.1
# return $512.2
# label 511
# del x.1
# del norm.1
# del $phi489.1
# del $phi486.1
# del $486.4
# jump 512
# $x489.2 = x.1 :: array(float64, 1d, C, nonconst)
# $k489.3 = k :: int64
# $489.6 = getitem(index=k, value=x.1) :: float64
# del k
# $489.8 = inplace_binop(rhs=norm.1, lhs=$489.6, fn=/?) :: float64
# del $489.6
# $x489.2[$k489.3] = $489.8 :: (array(float64, 1d, C, nonconst), int64, float64) -> none
# del $x489.2
# del $k489.3
# del $489.8
# jump 486
x[k] /= norm
# --- LINE 146 ---
# --- LINE 147 ---
return x
================================================================================
In [16]:
Content source: oyamad/gth_solve_numba
Similar notebooks: