In [1]:
import numpy as np
from numba import jit
from gth_solve_jit import gth_solve, gth_solve_jit
In [2]:
@jit('float64[:](float64[:,:])')
def gth_solve_jit2(A):
A1 = np.array(A, dtype=np.float64)
if len(A1.shape) != 2 or A1.shape[0] != A1.shape[1]:
# raise ValueError('matrix must be square') # Not supported
raise ValueError
n = A1.shape[0]
x = np.zeros(n, dtype=np.float64)
# === Reduction === #
for k in range(n-1):
scale = np.sum(A1[k, k+1:n])
#if scale <= 0:
# # There is one (and only one) recurrent class contained in
# # {0, ..., k};
# # compute the solution associated with that recurrent class.
# n = k+1
# break
for i in range(k+1, n):
A1[i, k] /= scale
for j in range(k+1, n):
A1[i, j] += A1[i, k] * A1[k, j]
# === Backward substitution === #
x[n-1] = 1
for k in range(n-2, -1, -1):
for i in range(k+1, n):
x[k] += x[i] * A1[i, k]
# === Normalization === #
norm = np.sum(x)
for k in range(n):
x[k] /= norm
return x
In [3]:
gth_solve_jit2([[0.4, 0.6], [0.2, 0.8]])
Out[3]:
array([ 0.25, 0.75])
In [4]:
gth_solve_jit2.inspect_types()
gth_solve_jit2 (array(float64, 2d, A, nonconst),)
--------------------------------------------------------------------------------
# File: <ipython-input-2-e178d06ec04e>
# --- LINE 1 ---
@jit('float64[:](float64[:,:])')
# --- LINE 2 ---
def gth_solve_jit2(A):
# --- LINE 3 ---
# 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(np: <module 'numpy' from '/usr/local/lib/python2.7/site-packages/numpy/__init__.pyc'>) :: pyobject
# $0.6 = getattr(attr=float64, value=$0.5) :: pyobject
# del $0.5
# $0.7 = call $0.2(A.1, dtype=$0.6) :: pyobject
# del A.1
# del $0.6
# del $0.2
# A1 = $0.7 :: pyobject
# del $0.7
A1 = np.array(A, dtype=np.float64)
# --- LINE 4 ---
# --- LINE 5 ---
# $0.8 = global(len: <built-in function len>) :: pyobject
# $0.10 = getattr(attr=shape, value=A1) :: pyobject
# $0.11 = call $0.8($0.10, ) :: pyobject
# del $0.8
# del $0.10
# $const0.12 = const(int, 2) :: pyobject
# $0.13 = $0.11 != $const0.12 :: pyobject
# del $const0.12
# del $0.11
# branch $0.13, 71, 45
# label 45
# $45.2 = getattr(attr=shape, value=A1) :: pyobject
# $const45.3 = const(int, 0) :: pyobject
# $45.4 = getitem(index=$const45.3, value=$45.2) :: pyobject
# del $const45.3
# del $45.2
# $45.6 = getattr(attr=shape, value=A1) :: pyobject
# $const45.7 = const(int, 1) :: pyobject
# $45.8 = getitem(index=$const45.7, value=$45.6) :: pyobject
# del $const45.7
# del $45.6
# $45.9 = $45.4 != $45.8 :: pyobject
# del $45.8
# del $45.4
# branch $45.9, 71, 80
if len(A1.shape) != 2 or A1.shape[0] != A1.shape[1]:
# --- LINE 6 ---
# raise ValueError('matrix must be square') # Not supported
# --- LINE 7 ---
# label 71
# del A1
# del $45.9
# del $0.13
# $71.1 = global(ValueError: <type 'exceptions.ValueError'>) :: pyobject
# del $71.1
# raise <type 'exceptions.ValueError'>
raise ValueError
# --- LINE 8 ---
# --- LINE 9 ---
# label 80
# del $45.9
# del $0.13
# $80.2 = getattr(attr=shape, value=A1) :: pyobject
# $const80.3 = const(int, 0) :: pyobject
# $80.4 = getitem(index=$const80.3, value=$80.2) :: pyobject
# del $const80.3
# del $80.2
# n = $80.4 :: pyobject
# del $80.4
n = A1.shape[0]
# --- LINE 10 ---
# --- LINE 11 ---
# $80.5 = global(np: <module 'numpy' from '/usr/local/lib/python2.7/site-packages/numpy/__init__.pyc'>) :: pyobject
# $80.6 = getattr(attr=zeros, value=$80.5) :: pyobject
# del $80.5
# $80.9 = global(np: <module 'numpy' from '/usr/local/lib/python2.7/site-packages/numpy/__init__.pyc'>) :: pyobject
# $80.10 = getattr(attr=float64, value=$80.9) :: pyobject
# del $80.9
# $80.11 = call $80.6(n, dtype=$80.10) :: pyobject
# del $80.6
# del $80.10
# x = $80.11 :: pyobject
# del $80.11
x = np.zeros(n, dtype=np.float64)
# --- LINE 12 ---
# --- LINE 13 ---
# === Reduction === #
# --- LINE 14 ---
# jump 120.1
# label 120.1
# $const120.1.1 = const(LiftedLoop, LiftedLoop(<function gth_solve_jit2 at 0x109b51758>)) :: pyobject
# $120.1.4 = call $const120.1.1(A1, n, ) :: pyobject
# del $const120.1.1
# del $120.1.4
# jump 307
for k in range(n-1):
# --- LINE 15 ---
scale = np.sum(A1[k, k+1:n])
# --- LINE 16 ---
#if scale <= 0:
# --- LINE 17 ---
# # There is one (and only one) recurrent class contained in
# --- LINE 18 ---
# # {0, ..., k};
# --- LINE 19 ---
# # compute the solution associated with that recurrent class.
# --- LINE 20 ---
# n = k+1
# --- LINE 21 ---
# break
# --- LINE 22 ---
for i in range(k+1, n):
# --- LINE 23 ---
A1[i, k] /= scale
# --- LINE 24 ---
# --- LINE 25 ---
for j in range(k+1, n):
# --- LINE 26 ---
A1[i, j] += A1[i, k] * A1[k, j]
# --- LINE 27 ---
# --- LINE 28 ---
# === Backward substitution === #
# --- LINE 29 ---
# label 307
# $const307.1 = const(int, 1) :: pyobject
# $const307.4 = const(int, 1) :: pyobject
# $307.5 = n - $const307.4 :: pyobject
# del $const307.4
# x[$307.5] = $const307.1 :: pyobject
# del $const307.1
# del $307.5
x[n-1] = 1
# --- LINE 30 ---
# jump 324.1
# label 324.1
# $const324.1.1 = const(LiftedLoop, LiftedLoop(<function gth_solve_jit2 at 0x109b51758>)) :: pyobject
# $324.1.5 = call $const324.1.1(A1, x, n, ) :: pyobject
# del A1
# del $const324.1.1
# del $324.1.5
# jump 418
for k in range(n-2, -1, -1):
# --- LINE 31 ---
for i in range(k+1, n):
# --- LINE 32 ---
x[k] += x[i] * A1[i, k]
# --- LINE 33 ---
# --- LINE 34 ---
# === Normalization === #
# --- LINE 35 ---
# label 418
# $418.1 = global(np: <module 'numpy' from '/usr/local/lib/python2.7/site-packages/numpy/__init__.pyc'>) :: pyobject
# $418.2 = getattr(attr=sum, value=$418.1) :: pyobject
# del $418.1
# $418.4 = call $418.2(x, ) :: pyobject
# del $418.2
# norm = $418.4 :: pyobject
# del $418.4
norm = np.sum(x)
# --- LINE 36 ---
# jump 436.1
# label 436.1
# $const436.1.1 = const(LiftedLoop, LiftedLoop(<function gth_solve_jit2 at 0x109b51758>)) :: pyobject
# $436.1.5 = call $const436.1.1(x, norm, n, ) :: pyobject
# del norm
# del n
# del $const436.1.1
# del $436.1.5
# jump 472
for k in range(n):
# --- LINE 37 ---
x[k] /= norm
# --- LINE 38 ---
# --- LINE 39 ---
# label 472
# $472.2 = cast(value=x) :: pyobject
# del x
# return $472.2
return x
# The function contains lifted loops
# Loop at line 14
# Has 1 overloads
# File: <ipython-input-2-e178d06ec04e>
# --- LINE 1 ---
@jit('float64[:](float64[:,:])')
# --- LINE 2 ---
def gth_solve_jit2(A):
# --- LINE 3 ---
A1 = np.array(A, dtype=np.float64)
# --- LINE 4 ---
# --- LINE 5 ---
if len(A1.shape) != 2 or A1.shape[0] != A1.shape[1]:
# --- LINE 6 ---
# raise ValueError('matrix must be square') # Not supported
# --- LINE 7 ---
raise ValueError
# --- LINE 8 ---
# --- LINE 9 ---
n = A1.shape[0]
# --- LINE 10 ---
# --- LINE 11 ---
x = np.zeros(n, dtype=np.float64)
# --- LINE 12 ---
# --- LINE 13 ---
# === Reduction === #
# --- LINE 14 ---
# label 134
# $134.2 = iternext(value=$phi134.1) :: pair<int64, bool>
# $134.3 = pair_first(value=$134.2) :: int64
# $134.4 = pair_second(value=$134.2) :: bool
# del $134.2
# $phi137.1 = $134.3 :: int64
# del $134.3
# branch $134.4, 137, 306
# label 137
# k = $phi137.1 :: int64
# del $phi137.1
# label 117
# A1.1 = A1 :: array(float64, 2d, C, nonconst)
# del A1
# n.1 = n :: int64
# del n
# $117.1 = global(range: <built-in function range>) :: range
# $const117.3 = const(int, 1) :: int32
# $117.4 = n.1 - $const117.3 :: int64
# del $const117.3
# $117.5 = call $117.1($117.4, ) :: (int64,) -> range_state64
# del $117.4
# del $117.1
# $117.6 = getiter(value=$117.5) :: range_iter64
# del $117.5
# $phi134.1 = $117.6 :: range_iter64
# del $117.6
# jump 134
for k in range(n-1):
# --- LINE 15 ---
# $137.2 = 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'>)
# $137.3 = getattr(attr=sum, value=$137.2) :: Function(<unbound method Numpy_reduce_<function sum at 0x1027c47d0>.sum>)
# del $137.2
# $k137.5 = k :: int64
# $const137.7 = const(int, 1) :: int32
# $137.8 = k + $const137.7 :: int64
# del $const137.7
# $137.10 = global(slice: <type 'slice'>) :: slice
# $137.11 = call $137.10($137.8, n.1, ) :: (int64, int64) -> slice3_type
# del $137.8
# del $137.10
# $137.12 = build_tuple(items=[Var($k137.5, <ipython-input-2-e178d06ec04e> (15)), Var($137.11, <ipython-input-2-e178d06ec04e> (15))]) :: (int64, slice3_type)
# del $k137.5
# del $137.11
# $137.13 = getitem(index=$137.12, value=A1.1) :: array(float64, 1d, A, nonconst)
# del $137.12
# $137.14 = call $137.3($137.13, ) :: (array(float64, 1d, A, nonconst),) -> float64
# del $137.3
# del $137.13
# scale = $137.14 :: float64
# del $137.14
scale = np.sum(A1[k, k+1:n])
# --- LINE 16 ---
#if scale <= 0:
# --- LINE 17 ---
# # There is one (and only one) recurrent class contained in
# --- LINE 18 ---
# # {0, ..., k};
# --- LINE 19 ---
# # compute the solution associated with that recurrent class.
# --- LINE 20 ---
# n = k+1
# --- LINE 21 ---
# break
# --- LINE 22 ---
# label 195
# $195.2 = iternext(value=$phi195.1) :: pair<int64, bool>
# $195.3 = pair_first(value=$195.2) :: int64
# $195.4 = pair_second(value=$195.2) :: bool
# del $195.2
# $phi198.1 = $195.3 :: int64
# del $195.3
# branch $195.4, 198, 302
# label 198
# i = $phi198.1 :: int64
# del $phi198.1
# jump 175
# label 175
# $175.1 = global(range: <built-in function range>) :: range
# $const175.3 = const(int, 1) :: int32
# $175.4 = k + $const175.3 :: int64
# del $const175.3
# $175.6 = call $175.1($175.4, n.1, ) :: (int64, int64) -> range_state64
# del $175.4
# del $175.1
# $175.7 = getiter(value=$175.6) :: range_iter64
# del $175.6
# $phi195.1 = $175.7 :: range_iter64
# del $175.7
# jump 195
for i in range(k+1, n):
# --- LINE 23 ---
# $A1198.2 = A1.1 :: array(float64, 2d, C, nonconst)
# $198.5 = build_tuple(items=[Var(i, <ipython-input-2-e178d06ec04e> (22)), Var(k, <ipython-input-2-e178d06ec04e> (14))]) :: (int64 x 2)
# $198.8 = getitem(index=$198.5, value=A1.1) :: float64
# $198.10 = inplace_binop(rhs=scale, lhs=$198.8, fn=/?) :: float64
# del $198.8
# $A1198.2[$198.5] = $198.10 :: (array(float64, 2d, C, nonconst), (int64 x 2), float64) -> none
# del $A1198.2
# del $198.5
# del $198.10
A1[i, k] /= scale
# --- LINE 24 ---
# --- LINE 25 ---
# jump 223
# label 243
# $243.2 = iternext(value=$phi243.1) :: pair<int64, bool>
# $243.3 = pair_first(value=$243.2) :: int64
# $243.4 = pair_second(value=$243.2) :: bool
# del $243.2
# $phi246.1 = $243.3 :: int64
# del $243.3
# branch $243.4, 246, 298
# label 246
# j = $phi246.1 :: int64
# del $phi246.1
# label 223
# $223.1 = global(range: <built-in function range>) :: range
# $const223.3 = const(int, 1) :: int32
# $223.4 = k + $const223.3 :: int64
# del $const223.3
# $223.6 = call $223.1($223.4, n.1, ) :: (int64, int64) -> range_state64
# del $223.4
# del $223.1
# $223.7 = getiter(value=$223.6) :: range_iter64
# del $223.6
# $phi243.1 = $223.7 :: range_iter64
# del $223.7
# jump 243
for j in range(k+1, n):
# --- LINE 26 ---
# label 307
# $const307.1 = const(NoneType, None) :: none
# $307.2 = cast(value=$const307.1) :: none
# del $const307.1
# return $307.2
# label 298
# del i
# del $phi246.1
# del $phi243.1
# del $243.4
# jump 299
# label 299
# jump 195
# label 302
# del scale
# del k
# del $phi198.1
# del $phi195.1
# del $195.4
# jump 303
# label 306
# del n.1
# del A1.1
# del $phi137.1
# del $phi134.1
# del $134.4
# jump 307
# $A1246.2 = A1.1 :: array(float64, 2d, C, nonconst)
# $246.5 = build_tuple(items=[Var(i, <ipython-input-2-e178d06ec04e> (22)), Var(j, <ipython-input-2-e178d06ec04e> (25))]) :: (int64 x 2)
# $246.8 = getitem(index=$246.5, value=A1.1) :: float64
# $246.12 = build_tuple(items=[Var(i, <ipython-input-2-e178d06ec04e> (22)), Var(k, <ipython-input-2-e178d06ec04e> (14))]) :: (int64 x 2)
# $246.13 = getitem(index=$246.12, value=A1.1) :: float64
# del $246.12
# $246.17 = build_tuple(items=[Var(k, <ipython-input-2-e178d06ec04e> (14)), Var(j, <ipython-input-2-e178d06ec04e> (25))]) :: (int64 x 2)
# del j
# $246.18 = getitem(index=$246.17, value=A1.1) :: float64
# del $246.17
# $246.19 = $246.13 * $246.18 :: float64
# del $246.18
# del $246.13
# $246.20 = inplace_binop(rhs=$246.19, lhs=$246.8, fn=+) :: float64
# del $246.8
# del $246.19
# $A1246.2[$246.5] = $246.20 :: (array(float64, 2d, C, nonconst), (int64 x 2), float64) -> none
# del $A1246.2
# del $246.5
# del $246.20
# jump 243
# label 303
# jump 134
A1[i, j] += A1[i, k] * A1[k, j]
# --- LINE 27 ---
# --- LINE 28 ---
# === Backward substitution === #
# --- LINE 29 ---
x[n-1] = 1
# --- LINE 30 ---
for k in range(n-2, -1, -1):
# --- LINE 31 ---
for i in range(k+1, n):
# --- LINE 32 ---
x[k] += x[i] * A1[i, k]
# --- LINE 33 ---
# --- LINE 34 ---
# === Normalization === #
# --- LINE 35 ---
norm = np.sum(x)
# --- LINE 36 ---
for k in range(n):
# --- LINE 37 ---
x[k] /= norm
# --- LINE 38 ---
# --- LINE 39 ---
return x
# Loop at line 30
# Has 1 overloads
# File: <ipython-input-2-e178d06ec04e>
# --- LINE 1 ---
@jit('float64[:](float64[:,:])')
# --- LINE 2 ---
def gth_solve_jit2(A):
# --- LINE 3 ---
A1 = np.array(A, dtype=np.float64)
# --- LINE 4 ---
# --- LINE 5 ---
if len(A1.shape) != 2 or A1.shape[0] != A1.shape[1]:
# --- LINE 6 ---
# raise ValueError('matrix must be square') # Not supported
# --- LINE 7 ---
raise ValueError
# --- LINE 8 ---
# --- LINE 9 ---
n = A1.shape[0]
# --- LINE 10 ---
# --- LINE 11 ---
x = np.zeros(n, dtype=np.float64)
# --- LINE 12 ---
# --- LINE 13 ---
# === Reduction === #
# --- LINE 14 ---
for k in range(n-1):
# --- LINE 15 ---
scale = np.sum(A1[k, k+1:n])
# --- LINE 16 ---
#if scale <= 0:
# --- LINE 17 ---
# # There is one (and only one) recurrent class contained in
# --- LINE 18 ---
# # {0, ..., k};
# --- LINE 19 ---
# # compute the solution associated with that recurrent class.
# --- LINE 20 ---
# n = k+1
# --- LINE 21 ---
# break
# --- LINE 22 ---
for i in range(k+1, n):
# --- LINE 23 ---
A1[i, k] /= scale
# --- LINE 24 ---
# --- LINE 25 ---
for j in range(k+1, n):
# --- LINE 26 ---
A1[i, j] += A1[i, k] * A1[k, j]
# --- LINE 27 ---
# --- LINE 28 ---
# === Backward substitution === #
# --- LINE 29 ---
x[n-1] = 1
# --- LINE 30 ---
# label 321
# 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
# $321.1 = global(range: <built-in function range>) :: range
# $const321.3 = const(int, 2) :: int32
# $321.4 = n.1 - $const321.3 :: int64
# del $const321.3
# $const321.5 = const(int, -1) :: int32
# $const321.6 = const(int, -1) :: int32
# $321.7 = call $321.1($321.4, $const321.5, $const321.6, ) :: (int64, int64, int64) -> range_state64
# del $const321.6
# del $const321.5
# del $321.4
# del $321.1
# $321.8 = getiter(value=$321.7) :: range_iter64
# del $321.7
# $phi344.1 = $321.8 :: range_iter64
# del $321.8
# jump 344
# label 344
# $344.2 = iternext(value=$phi344.1) :: pair<int64, bool>
# $344.3 = pair_first(value=$344.2) :: int64
# $344.4 = pair_second(value=$344.2) :: bool
# del $344.2
# $phi347.1 = $344.3 :: int64
# del $344.3
# branch $344.4, 347, 417
# label 347
# k = $phi347.1 :: int64
# del $phi347.1
for k in range(n-2, -1, -1):
# --- LINE 31 ---
# label 370
# $370.2 = iternext(value=$phi370.1) :: pair<int64, bool>
# $370.3 = pair_first(value=$370.2) :: int64
# $370.4 = pair_second(value=$370.2) :: bool
# del $370.2
# $phi373.1 = $370.3 :: int64
# del $370.3
# branch $370.4, 373, 413
# label 373
# i = $phi373.1 :: int64
# del $phi373.1
# jump 350
# label 350
# $350.1 = global(range: <built-in function range>) :: range
# $const350.3 = const(int, 1) :: int32
# $350.4 = k + $const350.3 :: int64
# del $const350.3
# $350.6 = call $350.1($350.4, n.1, ) :: (int64, int64) -> range_state64
# del $350.4
# del $350.1
# $350.7 = getiter(value=$350.6) :: range_iter64
# del $350.6
# $phi370.1 = $350.7 :: range_iter64
# del $350.7
# jump 370
for i in range(k+1, n):
# --- LINE 32 ---
# label 418
# $const418.1 = const(NoneType, None) :: none
# $418.2 = cast(value=$const418.1) :: none
# del $const418.1
# return $418.2
# label 417
# del x.1
# del n.1
# del A1.1
# del $phi347.1
# del $phi344.1
# del $344.4
# jump 418
# $x373.2 = x.1 :: array(float64, 1d, C, nonconst)
# $k373.3 = k :: int64
# $373.6 = getitem(index=k, value=x.1) :: float64
# $373.9 = getitem(index=i, value=x.1) :: float64
# $373.13 = build_tuple(items=[Var(i, <ipython-input-2-e178d06ec04e> (31)), Var(k, <ipython-input-2-e178d06ec04e> (30))]) :: (int64 x 2)
# del i
# $373.14 = getitem(index=$373.13, value=A1.1) :: float64
# del $373.13
# $373.15 = $373.9 * $373.14 :: float64
# del $373.9
# del $373.14
# $373.16 = inplace_binop(rhs=$373.15, lhs=$373.6, fn=+) :: float64
# del $373.6
# del $373.15
# $x373.2[$k373.3] = $373.16 :: (array(float64, 1d, C, nonconst), int64, float64) -> none
# del $x373.2
# del $k373.3
# del $373.16
# jump 370
# label 414
# jump 344
# label 413
# del k
# del $phi373.1
# del $phi370.1
# del $370.4
# jump 414
x[k] += x[i] * A1[i, k]
# --- LINE 33 ---
# --- LINE 34 ---
# === Normalization === #
# --- LINE 35 ---
norm = np.sum(x)
# --- LINE 36 ---
for k in range(n):
# --- LINE 37 ---
x[k] /= norm
# --- LINE 38 ---
# --- LINE 39 ---
return x
# Loop at line 36
# Has 1 overloads
# File: <ipython-input-2-e178d06ec04e>
# --- LINE 1 ---
@jit('float64[:](float64[:,:])')
# --- LINE 2 ---
def gth_solve_jit2(A):
# --- LINE 3 ---
A1 = np.array(A, dtype=np.float64)
# --- LINE 4 ---
# --- LINE 5 ---
if len(A1.shape) != 2 or A1.shape[0] != A1.shape[1]:
# --- LINE 6 ---
# raise ValueError('matrix must be square') # Not supported
# --- LINE 7 ---
raise ValueError
# --- LINE 8 ---
# --- LINE 9 ---
n = A1.shape[0]
# --- LINE 10 ---
# --- LINE 11 ---
x = np.zeros(n, dtype=np.float64)
# --- LINE 12 ---
# --- LINE 13 ---
# === Reduction === #
# --- LINE 14 ---
for k in range(n-1):
# --- LINE 15 ---
scale = np.sum(A1[k, k+1:n])
# --- LINE 16 ---
#if scale <= 0:
# --- LINE 17 ---
# # There is one (and only one) recurrent class contained in
# --- LINE 18 ---
# # {0, ..., k};
# --- LINE 19 ---
# # compute the solution associated with that recurrent class.
# --- LINE 20 ---
# n = k+1
# --- LINE 21 ---
# break
# --- LINE 22 ---
for i in range(k+1, n):
# --- LINE 23 ---
A1[i, k] /= scale
# --- LINE 24 ---
# --- LINE 25 ---
for j in range(k+1, n):
# --- LINE 26 ---
A1[i, j] += A1[i, k] * A1[k, j]
# --- LINE 27 ---
# --- LINE 28 ---
# === Backward substitution === #
# --- LINE 29 ---
x[n-1] = 1
# --- LINE 30 ---
for k in range(n-2, -1, -1):
# --- LINE 31 ---
for i in range(k+1, n):
# --- LINE 32 ---
x[k] += x[i] * A1[i, k]
# --- LINE 33 ---
# --- LINE 34 ---
# === Normalization === #
# --- LINE 35 ---
norm = np.sum(x)
# --- LINE 36 ---
# label 433
# x.1 = x :: array(float64, 1d, C, nonconst)
# del x
# norm.1 = norm :: float64
# del norm
# n.1 = n :: int64
# del n
# $433.1 = global(range: <built-in function range>) :: range
# $433.3 = call $433.1(n.1, ) :: (int64,) -> range_state64
# del n.1
# del $433.1
# $433.4 = getiter(value=$433.3) :: range_iter64
# del $433.3
# $phi446.1 = $433.4 :: range_iter64
# del $433.4
# jump 446
# label 446
# $446.2 = iternext(value=$phi446.1) :: pair<int64, bool>
# $446.3 = pair_first(value=$446.2) :: int64
# $446.4 = pair_second(value=$446.2) :: bool
# del $446.2
# $phi449.1 = $446.3 :: int64
# del $446.3
# branch $446.4, 449, 471
# label 449
# k = $phi449.1 :: int64
# del $phi449.1
for k in range(n):
# --- LINE 37 ---
# label 472
# $const472.1 = const(NoneType, None) :: none
# $472.2 = cast(value=$const472.1) :: none
# del $const472.1
# return $472.2
# label 471
# del x.1
# del norm.1
# del $phi449.1
# del $phi446.1
# del $446.4
# jump 472
# $x449.2 = x.1 :: array(float64, 1d, C, nonconst)
# $k449.3 = k :: int64
# $449.6 = getitem(index=k, value=x.1) :: float64
# del k
# $449.8 = inplace_binop(rhs=norm.1, lhs=$449.6, fn=/?) :: float64
# del $449.6
# $x449.2[$k449.3] = $449.8 :: (array(float64, 1d, C, nonconst), int64, float64) -> none
# del $x449.2
# del $k449.3
# del $449.8
# jump 446
x[k] /= norm
# --- LINE 38 ---
# --- LINE 39 ---
return x
================================================================================
In [5]:
sizes = [10, 50, 100, 1000]
rand_matrices = []
for n in sizes:
Q = np.random.rand(n, n)
Q /= np.sum(Q, axis=1, keepdims=True)
rand_matrices.append(Q)
In [6]:
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)
%timeit gth_solve_jit2(Q)
rand_matrices[0] (10 x 10)
1000 loops, best of 3: 240 µs per loop
100000 loops, best of 3: 7.01 µs per loop
100000 loops, best of 3: 17.6 µs per loop
rand_matrices[1] (50 x 50)
1000 loops, best of 3: 1.54 ms per loop
10000 loops, best of 3: 78.2 µs per loop
10000 loops, best of 3: 99.8 µs per loop
rand_matrices[2] (100 x 100)
100 loops, best of 3: 4.42 ms per loop
1000 loops, best of 3: 605 µs per loop
1000 loops, best of 3: 650 µs per loop
rand_matrices[3] (1000 x 1000)
1 loops, best of 3: 1.81 s per loop
1 loops, best of 3: 678 ms per loop
1 loops, best of 3: 634 ms per loop
In [7]:
import platform
print platform.platform()
Darwin-10.8.0-i386-64bit
In [8]:
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 [9]:
print np.__version__
1.9.0
In [10]:
import numba
print numba.__version__
0.15.1
In [11]:
import llvm
print llvm.__version__
0.12.7-6-gfad9ce0
In [12]:
@jit
def gth_solve_jit3(A):
A1 = np.array(A, dtype=np.float64)
if len(A1.shape) != 2 or A1.shape[0] != A1.shape[1]:
# raise ValueError('matrix must be square') # Not supported
raise ValueError
n = A1.shape[0]
x = np.zeros(n, dtype=np.float64)
# === Reduction === #
for k in range(n-1):
scale = np.sum(A1[k, k+1:n])
#if scale <= 0:
# # There is one (and only one) recurrent class contained in
# # {0, ..., k};
# # compute the solution associated with that recurrent class.
# n = k+1
# break
for i in range(k+1, n):
A1[i, k] /= scale
for j in range(k+1, n):
A1[i, j] += A1[i, k] * A1[k, j]
# === Backward substitution === #
x[n-1] = 1
for k in range(n-2, -1, -1):
for i in range(k+1, n):
x[k] += x[i] * A1[i, k]
# === Normalization === #
norm = np.sum(x)
for k in range(n):
x[k] /= norm
return x
In [13]:
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_jit(Q)
%timeit gth_solve_jit2(Q)
%timeit gth_solve_jit3(Q)
rand_matrices[0] (10 x 10)
100000 loops, best of 3: 6.96 µs per loop
100000 loops, best of 3: 17.6 µs per loop
1 loops, best of 3: 24.1 µs per loop
rand_matrices[1] (50 x 50)
10000 loops, best of 3: 78.2 µs per loop
10000 loops, best of 3: 100 µs per loop
10000 loops, best of 3: 99.6 µs per loop
rand_matrices[2] (100 x 100)
1000 loops, best of 3: 603 µs per loop
1000 loops, best of 3: 651 µs per loop
1000 loops, best of 3: 650 µs per loop
rand_matrices[3] (1000 x 1000)
1 loops, best of 3: 633 ms per loop
1 loops, best of 3: 624 ms per loop
1 loops, best of 3: 636 ms per loop
In [14]:
gth_solve_jit3.inspect_types()
gth_solve_jit3 (array(float64, 2d, C, nonconst),)
--------------------------------------------------------------------------------
# File: <ipython-input-12-f7e6a02b961c>
# --- LINE 1 ---
@jit
# --- LINE 2 ---
def gth_solve_jit3(A):
# --- LINE 3 ---
# 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(np: <module 'numpy' from '/usr/local/lib/python2.7/site-packages/numpy/__init__.pyc'>) :: pyobject
# $0.6 = getattr(attr=float64, value=$0.5) :: pyobject
# del $0.5
# $0.7 = call $0.2(A.1, dtype=$0.6) :: pyobject
# del A.1
# del $0.6
# del $0.2
# A1 = $0.7 :: pyobject
# del $0.7
A1 = np.array(A, dtype=np.float64)
# --- LINE 4 ---
# --- LINE 5 ---
# $0.8 = global(len: <built-in function len>) :: pyobject
# $0.10 = getattr(attr=shape, value=A1) :: pyobject
# $0.11 = call $0.8($0.10, ) :: pyobject
# del $0.8
# del $0.10
# $const0.12 = const(int, 2) :: pyobject
# $0.13 = $0.11 != $const0.12 :: pyobject
# del $const0.12
# del $0.11
# branch $0.13, 71, 45
# label 45
# $45.2 = getattr(attr=shape, value=A1) :: pyobject
# $const45.3 = const(int, 0) :: pyobject
# $45.4 = getitem(index=$const45.3, value=$45.2) :: pyobject
# del $const45.3
# del $45.2
# $45.6 = getattr(attr=shape, value=A1) :: pyobject
# $const45.7 = const(int, 1) :: pyobject
# $45.8 = getitem(index=$const45.7, value=$45.6) :: pyobject
# del $const45.7
# del $45.6
# $45.9 = $45.4 != $45.8 :: pyobject
# del $45.8
# del $45.4
# branch $45.9, 71, 80
if len(A1.shape) != 2 or A1.shape[0] != A1.shape[1]:
# --- LINE 6 ---
# raise ValueError('matrix must be square') # Not supported
# --- LINE 7 ---
# label 71
# del A1
# del $45.9
# del $0.13
# $71.1 = global(ValueError: <type 'exceptions.ValueError'>) :: pyobject
# del $71.1
# raise <type 'exceptions.ValueError'>
raise ValueError
# --- LINE 8 ---
# --- LINE 9 ---
# label 80
# del $45.9
# del $0.13
# $80.2 = getattr(attr=shape, value=A1) :: pyobject
# $const80.3 = const(int, 0) :: pyobject
# $80.4 = getitem(index=$const80.3, value=$80.2) :: pyobject
# del $const80.3
# del $80.2
# n = $80.4 :: pyobject
# del $80.4
n = A1.shape[0]
# --- LINE 10 ---
# --- LINE 11 ---
# $80.5 = global(np: <module 'numpy' from '/usr/local/lib/python2.7/site-packages/numpy/__init__.pyc'>) :: pyobject
# $80.6 = getattr(attr=zeros, value=$80.5) :: pyobject
# del $80.5
# $80.9 = global(np: <module 'numpy' from '/usr/local/lib/python2.7/site-packages/numpy/__init__.pyc'>) :: pyobject
# $80.10 = getattr(attr=float64, value=$80.9) :: pyobject
# del $80.9
# $80.11 = call $80.6(n, dtype=$80.10) :: pyobject
# del $80.6
# del $80.10
# x = $80.11 :: pyobject
# del $80.11
x = np.zeros(n, dtype=np.float64)
# --- LINE 12 ---
# --- LINE 13 ---
# === Reduction === #
# --- LINE 14 ---
# jump 120.1
# label 120.1
# $const120.1.1 = const(LiftedLoop, LiftedLoop(<function gth_solve_jit3 at 0x109b908c0>)) :: pyobject
# $120.1.4 = call $const120.1.1(A1, n, ) :: pyobject
# del $const120.1.1
# del $120.1.4
# jump 307
for k in range(n-1):
# --- LINE 15 ---
scale = np.sum(A1[k, k+1:n])
# --- LINE 16 ---
#if scale <= 0:
# --- LINE 17 ---
# # There is one (and only one) recurrent class contained in
# --- LINE 18 ---
# # {0, ..., k};
# --- LINE 19 ---
# # compute the solution associated with that recurrent class.
# --- LINE 20 ---
# n = k+1
# --- LINE 21 ---
# break
# --- LINE 22 ---
for i in range(k+1, n):
# --- LINE 23 ---
A1[i, k] /= scale
# --- LINE 24 ---
# --- LINE 25 ---
for j in range(k+1, n):
# --- LINE 26 ---
A1[i, j] += A1[i, k] * A1[k, j]
# --- LINE 27 ---
# --- LINE 28 ---
# === Backward substitution === #
# --- LINE 29 ---
# label 307
# $const307.1 = const(int, 1) :: pyobject
# $const307.4 = const(int, 1) :: pyobject
# $307.5 = n - $const307.4 :: pyobject
# del $const307.4
# x[$307.5] = $const307.1 :: pyobject
# del $const307.1
# del $307.5
x[n-1] = 1
# --- LINE 30 ---
# jump 324.1
# label 324.1
# $const324.1.1 = const(LiftedLoop, LiftedLoop(<function gth_solve_jit3 at 0x109b908c0>)) :: pyobject
# $324.1.5 = call $const324.1.1(A1, x, n, ) :: pyobject
# del A1
# del $const324.1.1
# del $324.1.5
# jump 418
for k in range(n-2, -1, -1):
# --- LINE 31 ---
for i in range(k+1, n):
# --- LINE 32 ---
x[k] += x[i] * A1[i, k]
# --- LINE 33 ---
# --- LINE 34 ---
# === Normalization === #
# --- LINE 35 ---
# label 418
# $418.1 = global(np: <module 'numpy' from '/usr/local/lib/python2.7/site-packages/numpy/__init__.pyc'>) :: pyobject
# $418.2 = getattr(attr=sum, value=$418.1) :: pyobject
# del $418.1
# $418.4 = call $418.2(x, ) :: pyobject
# del $418.2
# norm = $418.4 :: pyobject
# del $418.4
norm = np.sum(x)
# --- LINE 36 ---
# jump 436.1
# label 436.1
# $const436.1.1 = const(LiftedLoop, LiftedLoop(<function gth_solve_jit3 at 0x109b908c0>)) :: pyobject
# $436.1.5 = call $const436.1.1(x, norm, n, ) :: pyobject
# del norm
# del n
# del $const436.1.1
# del $436.1.5
# jump 472
for k in range(n):
# --- LINE 37 ---
x[k] /= norm
# --- LINE 38 ---
# --- LINE 39 ---
# label 472
# $472.2 = cast(value=x) :: pyobject
# del x
# return $472.2
return x
# The function contains lifted loops
# Loop at line 14
# Has 1 overloads
# File: <ipython-input-12-f7e6a02b961c>
# --- LINE 1 ---
@jit
# --- LINE 2 ---
def gth_solve_jit3(A):
# --- LINE 3 ---
A1 = np.array(A, dtype=np.float64)
# --- LINE 4 ---
# --- LINE 5 ---
if len(A1.shape) != 2 or A1.shape[0] != A1.shape[1]:
# --- LINE 6 ---
# raise ValueError('matrix must be square') # Not supported
# --- LINE 7 ---
raise ValueError
# --- LINE 8 ---
# --- LINE 9 ---
n = A1.shape[0]
# --- LINE 10 ---
# --- LINE 11 ---
x = np.zeros(n, dtype=np.float64)
# --- LINE 12 ---
# --- LINE 13 ---
# === Reduction === #
# --- LINE 14 ---
# label 134
# $134.2 = iternext(value=$phi134.1) :: pair<int64, bool>
# $134.3 = pair_first(value=$134.2) :: int64
# $134.4 = pair_second(value=$134.2) :: bool
# del $134.2
# $phi137.1 = $134.3 :: int64
# del $134.3
# branch $134.4, 137, 306
# label 137
# k = $phi137.1 :: int64
# del $phi137.1
# label 117
# A1.1 = A1 :: array(float64, 2d, C, nonconst)
# del A1
# n.1 = n :: int64
# del n
# $117.1 = global(range: <built-in function range>) :: range
# $const117.3 = const(int, 1) :: int32
# $117.4 = n.1 - $const117.3 :: int64
# del $const117.3
# $117.5 = call $117.1($117.4, ) :: (int64,) -> range_state64
# del $117.4
# del $117.1
# $117.6 = getiter(value=$117.5) :: range_iter64
# del $117.5
# $phi134.1 = $117.6 :: range_iter64
# del $117.6
# jump 134
for k in range(n-1):
# --- LINE 15 ---
# $137.2 = 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'>)
# $137.3 = getattr(attr=sum, value=$137.2) :: Function(<unbound method Numpy_reduce_<function sum at 0x1027c47d0>.sum>)
# del $137.2
# $k137.5 = k :: int64
# $const137.7 = const(int, 1) :: int32
# $137.8 = k + $const137.7 :: int64
# del $const137.7
# $137.10 = global(slice: <type 'slice'>) :: slice
# $137.11 = call $137.10($137.8, n.1, ) :: (int64, int64) -> slice3_type
# del $137.8
# del $137.10
# $137.12 = build_tuple(items=[Var($k137.5, <ipython-input-12-f7e6a02b961c> (15)), Var($137.11, <ipython-input-12-f7e6a02b961c> (15))]) :: (int64, slice3_type)
# del $k137.5
# del $137.11
# $137.13 = getitem(index=$137.12, value=A1.1) :: array(float64, 1d, A, nonconst)
# del $137.12
# $137.14 = call $137.3($137.13, ) :: (array(float64, 1d, A, nonconst),) -> float64
# del $137.3
# del $137.13
# scale = $137.14 :: float64
# del $137.14
scale = np.sum(A1[k, k+1:n])
# --- LINE 16 ---
#if scale <= 0:
# --- LINE 17 ---
# # There is one (and only one) recurrent class contained in
# --- LINE 18 ---
# # {0, ..., k};
# --- LINE 19 ---
# # compute the solution associated with that recurrent class.
# --- LINE 20 ---
# n = k+1
# --- LINE 21 ---
# break
# --- LINE 22 ---
# label 195
# $195.2 = iternext(value=$phi195.1) :: pair<int64, bool>
# $195.3 = pair_first(value=$195.2) :: int64
# $195.4 = pair_second(value=$195.2) :: bool
# del $195.2
# $phi198.1 = $195.3 :: int64
# del $195.3
# branch $195.4, 198, 302
# label 198
# i = $phi198.1 :: int64
# del $phi198.1
# jump 175
# label 175
# $175.1 = global(range: <built-in function range>) :: range
# $const175.3 = const(int, 1) :: int32
# $175.4 = k + $const175.3 :: int64
# del $const175.3
# $175.6 = call $175.1($175.4, n.1, ) :: (int64, int64) -> range_state64
# del $175.4
# del $175.1
# $175.7 = getiter(value=$175.6) :: range_iter64
# del $175.6
# $phi195.1 = $175.7 :: range_iter64
# del $175.7
# jump 195
for i in range(k+1, n):
# --- LINE 23 ---
# $A1198.2 = A1.1 :: array(float64, 2d, C, nonconst)
# $198.5 = build_tuple(items=[Var(i, <ipython-input-12-f7e6a02b961c> (22)), Var(k, <ipython-input-12-f7e6a02b961c> (14))]) :: (int64 x 2)
# $198.8 = getitem(index=$198.5, value=A1.1) :: float64
# $198.10 = inplace_binop(rhs=scale, lhs=$198.8, fn=/?) :: float64
# del $198.8
# $A1198.2[$198.5] = $198.10 :: (array(float64, 2d, C, nonconst), (int64 x 2), float64) -> none
# del $A1198.2
# del $198.5
# del $198.10
A1[i, k] /= scale
# --- LINE 24 ---
# --- LINE 25 ---
# jump 223
# label 243
# $243.2 = iternext(value=$phi243.1) :: pair<int64, bool>
# $243.3 = pair_first(value=$243.2) :: int64
# $243.4 = pair_second(value=$243.2) :: bool
# del $243.2
# $phi246.1 = $243.3 :: int64
# del $243.3
# branch $243.4, 246, 298
# label 246
# j = $phi246.1 :: int64
# del $phi246.1
# label 223
# $223.1 = global(range: <built-in function range>) :: range
# $const223.3 = const(int, 1) :: int32
# $223.4 = k + $const223.3 :: int64
# del $const223.3
# $223.6 = call $223.1($223.4, n.1, ) :: (int64, int64) -> range_state64
# del $223.4
# del $223.1
# $223.7 = getiter(value=$223.6) :: range_iter64
# del $223.6
# $phi243.1 = $223.7 :: range_iter64
# del $223.7
# jump 243
for j in range(k+1, n):
# --- LINE 26 ---
# label 307
# $const307.1 = const(NoneType, None) :: none
# $307.2 = cast(value=$const307.1) :: none
# del $const307.1
# return $307.2
# label 298
# del i
# del $phi246.1
# del $phi243.1
# del $243.4
# jump 299
# label 299
# jump 195
# label 302
# del scale
# del k
# del $phi198.1
# del $phi195.1
# del $195.4
# jump 303
# label 306
# del n.1
# del A1.1
# del $phi137.1
# del $phi134.1
# del $134.4
# jump 307
# $A1246.2 = A1.1 :: array(float64, 2d, C, nonconst)
# $246.5 = build_tuple(items=[Var(i, <ipython-input-12-f7e6a02b961c> (22)), Var(j, <ipython-input-12-f7e6a02b961c> (25))]) :: (int64 x 2)
# $246.8 = getitem(index=$246.5, value=A1.1) :: float64
# $246.12 = build_tuple(items=[Var(i, <ipython-input-12-f7e6a02b961c> (22)), Var(k, <ipython-input-12-f7e6a02b961c> (14))]) :: (int64 x 2)
# $246.13 = getitem(index=$246.12, value=A1.1) :: float64
# del $246.12
# $246.17 = build_tuple(items=[Var(k, <ipython-input-12-f7e6a02b961c> (14)), Var(j, <ipython-input-12-f7e6a02b961c> (25))]) :: (int64 x 2)
# del j
# $246.18 = getitem(index=$246.17, value=A1.1) :: float64
# del $246.17
# $246.19 = $246.13 * $246.18 :: float64
# del $246.18
# del $246.13
# $246.20 = inplace_binop(rhs=$246.19, lhs=$246.8, fn=+) :: float64
# del $246.8
# del $246.19
# $A1246.2[$246.5] = $246.20 :: (array(float64, 2d, C, nonconst), (int64 x 2), float64) -> none
# del $A1246.2
# del $246.5
# del $246.20
# jump 243
# label 303
# jump 134
A1[i, j] += A1[i, k] * A1[k, j]
# --- LINE 27 ---
# --- LINE 28 ---
# === Backward substitution === #
# --- LINE 29 ---
x[n-1] = 1
# --- LINE 30 ---
for k in range(n-2, -1, -1):
# --- LINE 31 ---
for i in range(k+1, n):
# --- LINE 32 ---
x[k] += x[i] * A1[i, k]
# --- LINE 33 ---
# --- LINE 34 ---
# === Normalization === #
# --- LINE 35 ---
norm = np.sum(x)
# --- LINE 36 ---
for k in range(n):
# --- LINE 37 ---
x[k] /= norm
# --- LINE 38 ---
# --- LINE 39 ---
return x
# Loop at line 30
# Has 1 overloads
# File: <ipython-input-12-f7e6a02b961c>
# --- LINE 1 ---
@jit
# --- LINE 2 ---
def gth_solve_jit3(A):
# --- LINE 3 ---
A1 = np.array(A, dtype=np.float64)
# --- LINE 4 ---
# --- LINE 5 ---
if len(A1.shape) != 2 or A1.shape[0] != A1.shape[1]:
# --- LINE 6 ---
# raise ValueError('matrix must be square') # Not supported
# --- LINE 7 ---
raise ValueError
# --- LINE 8 ---
# --- LINE 9 ---
n = A1.shape[0]
# --- LINE 10 ---
# --- LINE 11 ---
x = np.zeros(n, dtype=np.float64)
# --- LINE 12 ---
# --- LINE 13 ---
# === Reduction === #
# --- LINE 14 ---
for k in range(n-1):
# --- LINE 15 ---
scale = np.sum(A1[k, k+1:n])
# --- LINE 16 ---
#if scale <= 0:
# --- LINE 17 ---
# # There is one (and only one) recurrent class contained in
# --- LINE 18 ---
# # {0, ..., k};
# --- LINE 19 ---
# # compute the solution associated with that recurrent class.
# --- LINE 20 ---
# n = k+1
# --- LINE 21 ---
# break
# --- LINE 22 ---
for i in range(k+1, n):
# --- LINE 23 ---
A1[i, k] /= scale
# --- LINE 24 ---
# --- LINE 25 ---
for j in range(k+1, n):
# --- LINE 26 ---
A1[i, j] += A1[i, k] * A1[k, j]
# --- LINE 27 ---
# --- LINE 28 ---
# === Backward substitution === #
# --- LINE 29 ---
x[n-1] = 1
# --- LINE 30 ---
# label 321
# 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
# $321.1 = global(range: <built-in function range>) :: range
# $const321.3 = const(int, 2) :: int32
# $321.4 = n.1 - $const321.3 :: int64
# del $const321.3
# $const321.5 = const(int, -1) :: int32
# $const321.6 = const(int, -1) :: int32
# $321.7 = call $321.1($321.4, $const321.5, $const321.6, ) :: (int64, int64, int64) -> range_state64
# del $const321.6
# del $const321.5
# del $321.4
# del $321.1
# $321.8 = getiter(value=$321.7) :: range_iter64
# del $321.7
# $phi344.1 = $321.8 :: range_iter64
# del $321.8
# jump 344
# label 344
# $344.2 = iternext(value=$phi344.1) :: pair<int64, bool>
# $344.3 = pair_first(value=$344.2) :: int64
# $344.4 = pair_second(value=$344.2) :: bool
# del $344.2
# $phi347.1 = $344.3 :: int64
# del $344.3
# branch $344.4, 347, 417
# label 347
# k = $phi347.1 :: int64
# del $phi347.1
for k in range(n-2, -1, -1):
# --- LINE 31 ---
# label 370
# $370.2 = iternext(value=$phi370.1) :: pair<int64, bool>
# $370.3 = pair_first(value=$370.2) :: int64
# $370.4 = pair_second(value=$370.2) :: bool
# del $370.2
# $phi373.1 = $370.3 :: int64
# del $370.3
# branch $370.4, 373, 413
# label 373
# i = $phi373.1 :: int64
# del $phi373.1
# jump 350
# label 350
# $350.1 = global(range: <built-in function range>) :: range
# $const350.3 = const(int, 1) :: int32
# $350.4 = k + $const350.3 :: int64
# del $const350.3
# $350.6 = call $350.1($350.4, n.1, ) :: (int64, int64) -> range_state64
# del $350.4
# del $350.1
# $350.7 = getiter(value=$350.6) :: range_iter64
# del $350.6
# $phi370.1 = $350.7 :: range_iter64
# del $350.7
# jump 370
for i in range(k+1, n):
# --- LINE 32 ---
# label 418
# $const418.1 = const(NoneType, None) :: none
# $418.2 = cast(value=$const418.1) :: none
# del $const418.1
# return $418.2
# label 417
# del x.1
# del n.1
# del A1.1
# del $phi347.1
# del $phi344.1
# del $344.4
# jump 418
# $x373.2 = x.1 :: array(float64, 1d, C, nonconst)
# $k373.3 = k :: int64
# $373.6 = getitem(index=k, value=x.1) :: float64
# $373.9 = getitem(index=i, value=x.1) :: float64
# $373.13 = build_tuple(items=[Var(i, <ipython-input-12-f7e6a02b961c> (31)), Var(k, <ipython-input-12-f7e6a02b961c> (30))]) :: (int64 x 2)
# del i
# $373.14 = getitem(index=$373.13, value=A1.1) :: float64
# del $373.13
# $373.15 = $373.9 * $373.14 :: float64
# del $373.9
# del $373.14
# $373.16 = inplace_binop(rhs=$373.15, lhs=$373.6, fn=+) :: float64
# del $373.6
# del $373.15
# $x373.2[$k373.3] = $373.16 :: (array(float64, 1d, C, nonconst), int64, float64) -> none
# del $x373.2
# del $k373.3
# del $373.16
# jump 370
# label 414
# jump 344
# label 413
# del k
# del $phi373.1
# del $phi370.1
# del $370.4
# jump 414
x[k] += x[i] * A1[i, k]
# --- LINE 33 ---
# --- LINE 34 ---
# === Normalization === #
# --- LINE 35 ---
norm = np.sum(x)
# --- LINE 36 ---
for k in range(n):
# --- LINE 37 ---
x[k] /= norm
# --- LINE 38 ---
# --- LINE 39 ---
return x
# Loop at line 36
# Has 1 overloads
# File: <ipython-input-12-f7e6a02b961c>
# --- LINE 1 ---
@jit
# --- LINE 2 ---
def gth_solve_jit3(A):
# --- LINE 3 ---
A1 = np.array(A, dtype=np.float64)
# --- LINE 4 ---
# --- LINE 5 ---
if len(A1.shape) != 2 or A1.shape[0] != A1.shape[1]:
# --- LINE 6 ---
# raise ValueError('matrix must be square') # Not supported
# --- LINE 7 ---
raise ValueError
# --- LINE 8 ---
# --- LINE 9 ---
n = A1.shape[0]
# --- LINE 10 ---
# --- LINE 11 ---
x = np.zeros(n, dtype=np.float64)
# --- LINE 12 ---
# --- LINE 13 ---
# === Reduction === #
# --- LINE 14 ---
for k in range(n-1):
# --- LINE 15 ---
scale = np.sum(A1[k, k+1:n])
# --- LINE 16 ---
#if scale <= 0:
# --- LINE 17 ---
# # There is one (and only one) recurrent class contained in
# --- LINE 18 ---
# # {0, ..., k};
# --- LINE 19 ---
# # compute the solution associated with that recurrent class.
# --- LINE 20 ---
# n = k+1
# --- LINE 21 ---
# break
# --- LINE 22 ---
for i in range(k+1, n):
# --- LINE 23 ---
A1[i, k] /= scale
# --- LINE 24 ---
# --- LINE 25 ---
for j in range(k+1, n):
# --- LINE 26 ---
A1[i, j] += A1[i, k] * A1[k, j]
# --- LINE 27 ---
# --- LINE 28 ---
# === Backward substitution === #
# --- LINE 29 ---
x[n-1] = 1
# --- LINE 30 ---
for k in range(n-2, -1, -1):
# --- LINE 31 ---
for i in range(k+1, n):
# --- LINE 32 ---
x[k] += x[i] * A1[i, k]
# --- LINE 33 ---
# --- LINE 34 ---
# === Normalization === #
# --- LINE 35 ---
norm = np.sum(x)
# --- LINE 36 ---
# label 433
# x.1 = x :: array(float64, 1d, C, nonconst)
# del x
# norm.1 = norm :: float64
# del norm
# n.1 = n :: int64
# del n
# $433.1 = global(range: <built-in function range>) :: range
# $433.3 = call $433.1(n.1, ) :: (int64,) -> range_state64
# del n.1
# del $433.1
# $433.4 = getiter(value=$433.3) :: range_iter64
# del $433.3
# $phi446.1 = $433.4 :: range_iter64
# del $433.4
# jump 446
# label 446
# $446.2 = iternext(value=$phi446.1) :: pair<int64, bool>
# $446.3 = pair_first(value=$446.2) :: int64
# $446.4 = pair_second(value=$446.2) :: bool
# del $446.2
# $phi449.1 = $446.3 :: int64
# del $446.3
# branch $446.4, 449, 471
# label 449
# k = $phi449.1 :: int64
# del $phi449.1
for k in range(n):
# --- LINE 37 ---
# label 472
# $const472.1 = const(NoneType, None) :: none
# $472.2 = cast(value=$const472.1) :: none
# del $const472.1
# return $472.2
# label 471
# del x.1
# del norm.1
# del $phi449.1
# del $phi446.1
# del $446.4
# jump 472
# $x449.2 = x.1 :: array(float64, 1d, C, nonconst)
# $k449.3 = k :: int64
# $449.6 = getitem(index=k, value=x.1) :: float64
# del k
# $449.8 = inplace_binop(rhs=norm.1, lhs=$449.6, fn=/?) :: float64
# del $449.6
# $x449.2[$k449.3] = $449.8 :: (array(float64, 1d, C, nonconst), int64, float64) -> none
# del $x449.2
# del $k449.3
# del $449.8
# jump 446
x[k] /= norm
# --- LINE 38 ---
# --- LINE 39 ---
return x
================================================================================
In [14]:
Content source: oyamad/gth_solve_numba
Similar notebooks: