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.inspect_types()
gth_solve_jit2 (array(float64, 2d, A, nonconst),)
--------------------------------------------------------------------------------
# File: <ipython-input-2-b14bcb5c0c5e>
# --- 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 ---
# label 134
# $134.2 = iternext(value=$phi134.1) :: pyobject
# $134.3 = pair_first(value=$134.2) :: pyobject
# $134.4 = pair_second(value=$134.2) :: pyobject
# del $134.2
# $phi137.1 = $134.3 :: pyobject
# del $134.3
# branch $134.4, 137, 332
# label 137
# k = $phi137.1 :: pyobject
# del $phi137.1
# jump 117
# label 117
# $117.1 = global(range: <built-in function range>) :: pyobject
# $const117.3 = const(int, 1) :: pyobject
# $117.4 = n - $const117.3 :: pyobject
# del $const117.3
# $117.5 = call $117.1($117.4, ) :: pyobject
# del $117.4
# del $117.1
# $117.6 = getiter(value=$117.5) :: pyobject
# del $117.5
# $phi134.1 = $117.6 :: pyobject
# 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'>) :: pyobject
# $137.3 = getattr(attr=sum, value=$137.2) :: pyobject
# del $137.2
# $k137.5 = k :: pyobject
# $const137.7 = const(int, 1) :: pyobject
# $137.8 = k + $const137.7 :: pyobject
# del $const137.7
# $137.10 = global(slice: <type 'slice'>) :: pyobject
# $137.11 = call $137.10($137.8, n, ) :: pyobject
# del $137.8
# del $137.10
# $137.12 = build_tuple(items=[Var($k137.5, <ipython-input-2-b14bcb5c0c5e> (15)), Var($137.11, <ipython-input-2-b14bcb5c0c5e> (15))]) :: pyobject
# del $k137.5
# del $137.11
# $137.13 = getitem(index=$137.12, value=A1) :: pyobject
# del $137.12
# $137.14 = call $137.3($137.13, ) :: pyobject
# del $137.3
# del $137.13
# scale = $137.14 :: pyobject
# del $137.14
scale = np.sum(A1[k, k+1:n])
# --- LINE 16 ---
# $const137.16 = const(int, 0) :: pyobject
# $137.17 = scale <= $const137.16 :: pyobject
# del $const137.16
# branch $137.17, 187, 201
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 ---
# label 187
# del scale
# del $phi134.1
# del $137.17
# del $134.4
# $const187.2 = const(int, 1) :: pyobject
# $187.3 = k + $const187.2 :: pyobject
# del $const187.2
# n = $187.3 :: pyobject
# del $187.3
n = k+1
# --- LINE 21 ---
# jump 333
break
# --- LINE 22 ---
# label 201
# $201.1 = global(range: <built-in function range>) :: pyobject
# $const201.3 = const(int, 1) :: pyobject
# $201.4 = k + $const201.3 :: pyobject
# del $const201.3
# $201.6 = call $201.1($201.4, n, ) :: pyobject
# del $201.4
# del $201.1
# $201.7 = getiter(value=$201.6) :: pyobject
# del $201.6
# $phi221.1 = $201.7 :: pyobject
# del $201.7
# jump 221
# label 221
# $221.2 = iternext(value=$phi221.1) :: pyobject
# $221.3 = pair_first(value=$221.2) :: pyobject
# $221.4 = pair_second(value=$221.2) :: pyobject
# del $221.2
# $phi224.1 = $221.3 :: pyobject
# del $221.3
# branch $221.4, 224, 328
# label 224
# i = $phi224.1 :: pyobject
# del $phi224.1
for i in range(k+1, n):
# --- LINE 23 ---
# $A1224.2 = A1 :: pyobject
# $224.5 = build_tuple(items=[Var(i, <ipython-input-2-b14bcb5c0c5e> (22)), Var(k, <ipython-input-2-b14bcb5c0c5e> (14))]) :: pyobject
# $224.8 = getitem(index=$224.5, value=A1) :: pyobject
# $224.10 = inplace_binop(rhs=scale, lhs=$224.8, fn=/?) :: pyobject
# del $224.8
# $A1224.2[$224.5] = $224.10 :: pyobject
# del $A1224.2
# del $224.5
# del $224.10
A1[i, k] /= scale
# --- LINE 24 ---
# --- LINE 25 ---
# label 269
# $269.2 = iternext(value=$phi269.1) :: pyobject
# $269.3 = pair_first(value=$269.2) :: pyobject
# $269.4 = pair_second(value=$269.2) :: pyobject
# del $269.2
# $phi272.1 = $269.3 :: pyobject
# del $269.3
# branch $269.4, 272, 324
# label 272
# j = $phi272.1 :: pyobject
# del $phi272.1
# jump 249
# label 249
# $249.1 = global(range: <built-in function range>) :: pyobject
# $const249.3 = const(int, 1) :: pyobject
# $249.4 = k + $const249.3 :: pyobject
# del $const249.3
# $249.6 = call $249.1($249.4, n, ) :: pyobject
# del $249.4
# del $249.1
# $249.7 = getiter(value=$249.6) :: pyobject
# del $249.6
# $phi269.1 = $249.7 :: pyobject
# del $249.7
# jump 269
for j in range(k+1, n):
# --- LINE 26 ---
# $A1272.2 = A1 :: pyobject
# $272.5 = build_tuple(items=[Var(i, <ipython-input-2-b14bcb5c0c5e> (22)), Var(j, <ipython-input-2-b14bcb5c0c5e> (25))]) :: pyobject
# $272.8 = getitem(index=$272.5, value=A1) :: pyobject
# $272.12 = build_tuple(items=[Var(i, <ipython-input-2-b14bcb5c0c5e> (22)), Var(k, <ipython-input-2-b14bcb5c0c5e> (14))]) :: pyobject
# $272.13 = getitem(index=$272.12, value=A1) :: pyobject
# del $272.12
# $272.17 = build_tuple(items=[Var(k, <ipython-input-2-b14bcb5c0c5e> (14)), Var(j, <ipython-input-2-b14bcb5c0c5e> (25))]) :: pyobject
# del j
# $272.18 = getitem(index=$272.17, value=A1) :: pyobject
# del $272.17
# $272.19 = $272.13 * $272.18 :: pyobject
# del $272.18
# del $272.13
# $272.20 = inplace_binop(rhs=$272.19, lhs=$272.8, fn=+) :: pyobject
# del $272.8
# del $272.19
# $A1272.2[$272.5] = $272.20 :: pyobject
# del $A1272.2
# del $272.5
# del $272.20
# jump 269
# label 329
# jump 134
# label 324
# del $phi272.1
# del $phi269.1
# del $269.4
# jump 325
# label 325
# jump 221
# label 328
# del scale
# del $phi224.1
# del $phi221.1
# del $221.4
# jump 329
# label 332
# del $phi137.1
# del $phi134.1
# del $137.17
# del $134.4
A1[i, j] += A1[i, k] * A1[k, j]
# --- LINE 27 ---
# --- LINE 28 ---
# === Backward substitution === #
# --- LINE 29 ---
# jump 333
# label 333
# $const333.1 = const(int, 1) :: pyobject
# $const333.4 = const(int, 1) :: pyobject
# $333.5 = n - $const333.4 :: pyobject
# del $const333.4
# x[$333.5] = $const333.1 :: pyobject
# del $const333.1
# del $333.5
x[n-1] = 1
# --- LINE 30 ---
# label 350.1
# $const350.1.1 = const(LiftedLoop, LiftedLoop(<function gth_solve_jit2 at 0x108e181b8>)) :: pyobject
# $350.1.7 = call $const350.1.1(A1, i, k, n, x, ) :: pyobject
# del A1
# del $const350.1.1
# $350.1.10 = exhaust_iter(count=2, value=$350.1.7) :: pyobject
# del $350.1.7
# $350.1.8 = static_getitem(index=0, value=$350.1.10) :: pyobject
# $350.1.9 = static_getitem(index=1, value=$350.1.10) :: pyobject
# del $350.1.10
# i = $350.1.8 :: pyobject
# del i
# del $350.1.8
# k = $350.1.9 :: pyobject
# del $350.1.9
# jump 444
# jump 350.1
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 444
# $444.1 = global(np: <module 'numpy' from '/usr/local/lib/python2.7/site-packages/numpy/__init__.pyc'>) :: pyobject
# $444.2 = getattr(attr=sum, value=$444.1) :: pyobject
# del $444.1
# $444.4 = call $444.2(x, ) :: pyobject
# del $444.2
# norm = $444.4 :: pyobject
# del $444.4
norm = np.sum(x)
# --- LINE 36 ---
# jump 462.1
# label 462.1
# $const462.1.1 = const(LiftedLoop, LiftedLoop(<function gth_solve_jit2 at 0x108e181b8>)) :: pyobject
# $462.1.6 = call $const462.1.1(x, k, norm, n, ) :: pyobject
# del norm
# del n
# del $const462.1.1
# $462.1.8 = exhaust_iter(count=1, value=$462.1.6) :: pyobject
# del $462.1.6
# $462.1.7 = static_getitem(index=0, value=$462.1.8) :: pyobject
# del $462.1.8
# k = $462.1.7 :: pyobject
# del k
# del $462.1.7
# jump 498
for k in range(n):
# --- LINE 37 ---
x[k] /= norm
# --- LINE 38 ---
# --- LINE 39 ---
# label 498
# $498.2 = cast(value=x) :: pyobject
# del x
# return $498.2
return x
# The function contains lifted loops
# Loop at line 30
# Has 0 overloads
# Loop at line 36
# Has 0 overloads
================================================================================
In [4]:
gth_solve_jit2([[0.4, 0.6], [0.2, 0.8]])
Out[4]:
array([ 0.25, 0.75])
In [5]:
gth_solve_jit2.inspect_types()
gth_solve_jit2 (array(float64, 2d, A, nonconst),)
--------------------------------------------------------------------------------
# File: <ipython-input-2-b14bcb5c0c5e>
# --- 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 ---
# label 134
# $134.2 = iternext(value=$phi134.1) :: pyobject
# $134.3 = pair_first(value=$134.2) :: pyobject
# $134.4 = pair_second(value=$134.2) :: pyobject
# del $134.2
# $phi137.1 = $134.3 :: pyobject
# del $134.3
# branch $134.4, 137, 332
# label 137
# k = $phi137.1 :: pyobject
# del $phi137.1
# jump 117
# label 117
# $117.1 = global(range: <built-in function range>) :: pyobject
# $const117.3 = const(int, 1) :: pyobject
# $117.4 = n - $const117.3 :: pyobject
# del $const117.3
# $117.5 = call $117.1($117.4, ) :: pyobject
# del $117.4
# del $117.1
# $117.6 = getiter(value=$117.5) :: pyobject
# del $117.5
# $phi134.1 = $117.6 :: pyobject
# 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'>) :: pyobject
# $137.3 = getattr(attr=sum, value=$137.2) :: pyobject
# del $137.2
# $k137.5 = k :: pyobject
# $const137.7 = const(int, 1) :: pyobject
# $137.8 = k + $const137.7 :: pyobject
# del $const137.7
# $137.10 = global(slice: <type 'slice'>) :: pyobject
# $137.11 = call $137.10($137.8, n, ) :: pyobject
# del $137.8
# del $137.10
# $137.12 = build_tuple(items=[Var($k137.5, <ipython-input-2-b14bcb5c0c5e> (15)), Var($137.11, <ipython-input-2-b14bcb5c0c5e> (15))]) :: pyobject
# del $k137.5
# del $137.11
# $137.13 = getitem(index=$137.12, value=A1) :: pyobject
# del $137.12
# $137.14 = call $137.3($137.13, ) :: pyobject
# del $137.3
# del $137.13
# scale = $137.14 :: pyobject
# del $137.14
scale = np.sum(A1[k, k+1:n])
# --- LINE 16 ---
# $const137.16 = const(int, 0) :: pyobject
# $137.17 = scale <= $const137.16 :: pyobject
# del $const137.16
# branch $137.17, 187, 201
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 ---
# label 187
# del scale
# del $phi134.1
# del $137.17
# del $134.4
# $const187.2 = const(int, 1) :: pyobject
# $187.3 = k + $const187.2 :: pyobject
# del $const187.2
# n = $187.3 :: pyobject
# del $187.3
n = k+1
# --- LINE 21 ---
# jump 333
break
# --- LINE 22 ---
# label 201
# $201.1 = global(range: <built-in function range>) :: pyobject
# $const201.3 = const(int, 1) :: pyobject
# $201.4 = k + $const201.3 :: pyobject
# del $const201.3
# $201.6 = call $201.1($201.4, n, ) :: pyobject
# del $201.4
# del $201.1
# $201.7 = getiter(value=$201.6) :: pyobject
# del $201.6
# $phi221.1 = $201.7 :: pyobject
# del $201.7
# jump 221
# label 221
# $221.2 = iternext(value=$phi221.1) :: pyobject
# $221.3 = pair_first(value=$221.2) :: pyobject
# $221.4 = pair_second(value=$221.2) :: pyobject
# del $221.2
# $phi224.1 = $221.3 :: pyobject
# del $221.3
# branch $221.4, 224, 328
# label 224
# i = $phi224.1 :: pyobject
# del $phi224.1
for i in range(k+1, n):
# --- LINE 23 ---
# $A1224.2 = A1 :: pyobject
# $224.5 = build_tuple(items=[Var(i, <ipython-input-2-b14bcb5c0c5e> (22)), Var(k, <ipython-input-2-b14bcb5c0c5e> (14))]) :: pyobject
# $224.8 = getitem(index=$224.5, value=A1) :: pyobject
# $224.10 = inplace_binop(rhs=scale, lhs=$224.8, fn=/?) :: pyobject
# del $224.8
# $A1224.2[$224.5] = $224.10 :: pyobject
# del $A1224.2
# del $224.5
# del $224.10
A1[i, k] /= scale
# --- LINE 24 ---
# --- LINE 25 ---
# label 269
# $269.2 = iternext(value=$phi269.1) :: pyobject
# $269.3 = pair_first(value=$269.2) :: pyobject
# $269.4 = pair_second(value=$269.2) :: pyobject
# del $269.2
# $phi272.1 = $269.3 :: pyobject
# del $269.3
# branch $269.4, 272, 324
# label 272
# j = $phi272.1 :: pyobject
# del $phi272.1
# jump 249
# label 249
# $249.1 = global(range: <built-in function range>) :: pyobject
# $const249.3 = const(int, 1) :: pyobject
# $249.4 = k + $const249.3 :: pyobject
# del $const249.3
# $249.6 = call $249.1($249.4, n, ) :: pyobject
# del $249.4
# del $249.1
# $249.7 = getiter(value=$249.6) :: pyobject
# del $249.6
# $phi269.1 = $249.7 :: pyobject
# del $249.7
# jump 269
for j in range(k+1, n):
# --- LINE 26 ---
# $A1272.2 = A1 :: pyobject
# $272.5 = build_tuple(items=[Var(i, <ipython-input-2-b14bcb5c0c5e> (22)), Var(j, <ipython-input-2-b14bcb5c0c5e> (25))]) :: pyobject
# $272.8 = getitem(index=$272.5, value=A1) :: pyobject
# $272.12 = build_tuple(items=[Var(i, <ipython-input-2-b14bcb5c0c5e> (22)), Var(k, <ipython-input-2-b14bcb5c0c5e> (14))]) :: pyobject
# $272.13 = getitem(index=$272.12, value=A1) :: pyobject
# del $272.12
# $272.17 = build_tuple(items=[Var(k, <ipython-input-2-b14bcb5c0c5e> (14)), Var(j, <ipython-input-2-b14bcb5c0c5e> (25))]) :: pyobject
# del j
# $272.18 = getitem(index=$272.17, value=A1) :: pyobject
# del $272.17
# $272.19 = $272.13 * $272.18 :: pyobject
# del $272.18
# del $272.13
# $272.20 = inplace_binop(rhs=$272.19, lhs=$272.8, fn=+) :: pyobject
# del $272.8
# del $272.19
# $A1272.2[$272.5] = $272.20 :: pyobject
# del $A1272.2
# del $272.5
# del $272.20
# jump 269
# label 329
# jump 134
# label 324
# del $phi272.1
# del $phi269.1
# del $269.4
# jump 325
# label 325
# jump 221
# label 328
# del scale
# del $phi224.1
# del $phi221.1
# del $221.4
# jump 329
# label 332
# del $phi137.1
# del $phi134.1
# del $137.17
# del $134.4
A1[i, j] += A1[i, k] * A1[k, j]
# --- LINE 27 ---
# --- LINE 28 ---
# === Backward substitution === #
# --- LINE 29 ---
# jump 333
# label 333
# $const333.1 = const(int, 1) :: pyobject
# $const333.4 = const(int, 1) :: pyobject
# $333.5 = n - $const333.4 :: pyobject
# del $const333.4
# x[$333.5] = $const333.1 :: pyobject
# del $const333.1
# del $333.5
x[n-1] = 1
# --- LINE 30 ---
# label 350.1
# $const350.1.1 = const(LiftedLoop, LiftedLoop(<function gth_solve_jit2 at 0x108e181b8>)) :: pyobject
# $350.1.7 = call $const350.1.1(A1, i, k, n, x, ) :: pyobject
# del A1
# del $const350.1.1
# $350.1.10 = exhaust_iter(count=2, value=$350.1.7) :: pyobject
# del $350.1.7
# $350.1.8 = static_getitem(index=0, value=$350.1.10) :: pyobject
# $350.1.9 = static_getitem(index=1, value=$350.1.10) :: pyobject
# del $350.1.10
# i = $350.1.8 :: pyobject
# del i
# del $350.1.8
# k = $350.1.9 :: pyobject
# del $350.1.9
# jump 444
# jump 350.1
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 444
# $444.1 = global(np: <module 'numpy' from '/usr/local/lib/python2.7/site-packages/numpy/__init__.pyc'>) :: pyobject
# $444.2 = getattr(attr=sum, value=$444.1) :: pyobject
# del $444.1
# $444.4 = call $444.2(x, ) :: pyobject
# del $444.2
# norm = $444.4 :: pyobject
# del $444.4
norm = np.sum(x)
# --- LINE 36 ---
# jump 462.1
# label 462.1
# $const462.1.1 = const(LiftedLoop, LiftedLoop(<function gth_solve_jit2 at 0x108e181b8>)) :: pyobject
# $462.1.6 = call $const462.1.1(x, k, norm, n, ) :: pyobject
# del norm
# del n
# del $const462.1.1
# $462.1.8 = exhaust_iter(count=1, value=$462.1.6) :: pyobject
# del $462.1.6
# $462.1.7 = static_getitem(index=0, value=$462.1.8) :: pyobject
# del $462.1.8
# k = $462.1.7 :: pyobject
# del k
# del $462.1.7
# jump 498
for k in range(n):
# --- LINE 37 ---
x[k] /= norm
# --- LINE 38 ---
# --- LINE 39 ---
# label 498
# $498.2 = cast(value=x) :: pyobject
# del x
# return $498.2
return x
# The function contains lifted loops
# Loop at line 30
# Has 1 overloads
# File: <ipython-input-2-b14bcb5c0c5e>
# --- 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 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, 443
# label 373
# k.1 = $phi373.1 :: int64
# del $phi373.1
# label 347
# A1.1 = A1 :: array(float64, 2d, C, nonconst)
# del A1
# i.1 = i :: int64
# del i
# k.1 = k :: int64
# del k
# n.1 = n :: int64
# del n
# x.1 = x :: array(float64, 1d, C, nonconst)
# del x
# $347.1 = global(range: <built-in function range>) :: range
# $const347.3 = const(int, 2) :: int32
# $347.4 = n.1 - $const347.3 :: int64
# del $const347.3
# $const347.5 = const(int, -1) :: int32
# $const347.6 = const(int, -1) :: int32
# $347.7 = call $347.1($347.4, $const347.5, $const347.6, ) :: (int64, int64, int64) -> range_state64
# del $const347.6
# del $const347.5
# del $347.4
# del $347.1
# $347.8 = getiter(value=$347.7) :: range_iter64
# del $347.7
# $phi370.1 = $347.8 :: range_iter64
# del $347.8
# jump 370
for k in range(n-2, -1, -1):
# --- LINE 31 ---
# label 396
# $396.2 = iternext(value=$phi396.1) :: pair<int64, bool>
# $396.3 = pair_first(value=$396.2) :: int64
# $396.4 = pair_second(value=$396.2) :: bool
# del $396.2
# $phi399.1 = $396.3 :: int64
# del $396.3
# branch $396.4, 399, 439
# label 399
# i.1 = $phi399.1 :: int64
# del $phi399.1
# jump 376
# label 376
# $376.1 = global(range: <built-in function range>) :: range
# $const376.3 = const(int, 1) :: int32
# $376.4 = k.1 + $const376.3 :: int64
# del $const376.3
# $376.6 = call $376.1($376.4, n.1, ) :: (int64, int64) -> range_state64
# del $376.4
# del $376.1
# $376.7 = getiter(value=$376.6) :: range_iter64
# del $376.6
# $phi396.1 = $376.7 :: range_iter64
# del $376.7
# jump 396
for i in range(k+1, n):
# --- LINE 32 ---
# label 443
# del x.1
# del n.1
# del A1.1
# del $phi373.1
# del $phi370.1
# del $370.4
# jump 444
# $x399.2 = x.1 :: array(float64, 1d, C, nonconst)
# $k399.3 = k.1 :: int64
# $399.6 = getitem(index=k.1, value=x.1) :: float64
# $399.9 = getitem(index=i.1, value=x.1) :: float64
# $399.13 = build_tuple(items=[Var(i.1, <ipython-input-2-b14bcb5c0c5e> (30)), Var(k.1, <ipython-input-2-b14bcb5c0c5e> (30))]) :: (int64 x 2)
# $399.14 = getitem(index=$399.13, value=A1.1) :: float64
# del $399.13
# $399.15 = $399.9 * $399.14 :: float64
# del $399.9
# del $399.14
# $399.16 = inplace_binop(rhs=$399.15, lhs=$399.6, fn=+) :: float64
# del $399.6
# del $399.15
# $x399.2[$k399.3] = $399.16 :: (array(float64, 1d, C, nonconst), int64, float64) -> none
# del $x399.2
# del $k399.3
# del $399.16
# jump 396
# label 440
# jump 370
# label 439
# del $phi399.1
# del $phi396.1
# del $396.4
# jump 440
# label 444
# $444.3 = build_tuple(items=[Var(i.1, <ipython-input-2-b14bcb5c0c5e> (30)), Var(k.1, <ipython-input-2-b14bcb5c0c5e> (30))]) :: (int64 x 2)
# del k.1
# del i.1
# $444.4 = cast(value=$444.3) :: (int64 x 2)
# del $444.3
# return $444.4
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-b14bcb5c0c5e>
# --- 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 472
# $472.2 = iternext(value=$phi472.1) :: pair<int64, bool>
# $472.3 = pair_first(value=$472.2) :: int64
# $472.4 = pair_second(value=$472.2) :: bool
# del $472.2
# $phi475.1 = $472.3 :: int64
# del $472.3
# branch $472.4, 475, 497
# label 459
# x.1 = x :: array(float64, 1d, C, nonconst)
# del x
# k.1 = k :: int64
# del k
# norm.1 = norm :: float64
# del norm
# n.1 = n :: int64
# del n
# $459.1 = global(range: <built-in function range>) :: range
# $459.3 = call $459.1(n.1, ) :: (int64,) -> range_state64
# del n.1
# del $459.1
# $459.4 = getiter(value=$459.3) :: range_iter64
# del $459.3
# $phi472.1 = $459.4 :: range_iter64
# del $459.4
# jump 472
# label 475
# k.1 = $phi475.1 :: int64
# del $phi475.1
for k in range(n):
# --- LINE 37 ---
# label 497
# del x.1
# del norm.1
# del $phi475.1
# del $phi472.1
# del $472.4
# jump 498
# label 498
# $498.2 = build_tuple(items=[Var(k.1, <ipython-input-2-b14bcb5c0c5e> (36))]) :: (int64 x 1)
# del k.1
# $498.3 = cast(value=$498.2) :: (int64 x 1)
# del $498.2
# return $498.3
# $x475.2 = x.1 :: array(float64, 1d, C, nonconst)
# $k475.3 = k.1 :: int64
# $475.6 = getitem(index=k.1, value=x.1) :: float64
# $475.8 = inplace_binop(rhs=norm.1, lhs=$475.6, fn=/?) :: float64
# del $475.6
# $x475.2[$k475.3] = $475.8 :: (array(float64, 1d, C, nonconst), int64, float64) -> none
# del $x475.2
# del $k475.3
# del $475.8
# jump 472
x[k] /= norm
# --- LINE 38 ---
# --- LINE 39 ---
return x
================================================================================
In [6]:
sizes = [10, 50, 100] # [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 [7]:
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: 182 µs per loop
100000 loops, best of 3: 5.2 µs per loop
1000 loops, best of 3: 281 µs per loop
rand_matrices[1] (50 x 50)
1000 loops, best of 3: 1.12 ms per loop
10000 loops, best of 3: 63.1 µs per loop
10 loops, best of 3: 21.1 ms per loop
rand_matrices[2] (100 x 100)
100 loops, best of 3: 3 ms per loop
1000 loops, best of 3: 429 µs per loop
10 loops, best of 3: 164 ms per loop
In [8]:
import platform
print platform.platform()
Darwin-13.4.0-x86_64-i386-64bit
In [9]:
import sys
print sys.version
2.7.8 (default, Jul 2 2014, 10:14:46)
[GCC 4.2.1 Compatible Apple LLVM 5.1 (clang-503.0.40)]
In [10]:
print np.__version__
1.9.0
In [11]:
import numba
print numba.__version__
0.15.1
In [12]:
import llvm
print llvm.__version__
0.12.7-5-gc0ae9c2
In [ ]:
Content source: oyamad/gth_solve_numba
Similar notebooks: