Skip to content
项目
群组
代码片段
帮助
当前项目
正在载入...
登录 / 注册
切换导航面板
P
pytensor
项目
项目
详情
活动
周期分析
仓库
仓库
文件
提交
分支
标签
贡献者
图表
比较
统计图
议题
0
议题
0
列表
看板
标记
里程碑
合并请求
0
合并请求
0
CI / CD
CI / CD
流水线
作业
日程
统计图
Wiki
Wiki
代码片段
代码片段
成员
成员
折叠边栏
关闭边栏
活动
图像
聊天
创建新问题
作业
提交
问题看板
Open sidebar
testgroup
pytensor
Commits
f28d5cb9
提交
f28d5cb9
authored
11月 20, 2008
作者:
James Bergstra
浏览文件
操作
浏览文件
下载
电子邮件补丁
差异文件
added dot and gemm-related optimizations to FAST_RUN from autoencoder/aa.py
上级
a489cb3e
隐藏空白字符变更
内嵌
并排
正在显示
3 个修改的文件
包含
242 行增加
和
235 行删除
+242
-235
aa.py
benchmark/autoencoder/aa.py
+0
-233
elemwise.py
theano/tensor/elemwise.py
+0
-2
opt.py
theano/tensor/opt.py
+242
-0
没有找到文件。
benchmark/autoencoder/aa.py
浏览文件 @
f28d5cb9
...
...
@@ -18,239 +18,6 @@ from theano.tensor import blas, opt
# c : aa.cc
class
_Dot22
(
Op
):
"""Compute matrix-matrix product.
This is a specialization of the more general Dot()
"""
def
make_node
(
self
,
x
,
y
):
assert
x
.
type
in
T
.
float_matrix_types
assert
y
.
type
==
x
.
type
bz
=
[
x
.
type
.
broadcastable
[
0
],
y
.
type
.
broadcastable
[
1
]]
outputs
=
[
T
.
tensor
(
x
.
type
.
dtype
,
bz
)]
return
Apply
(
self
,
[
x
,
y
],
outputs
)
def
perform
(
self
,
node
,
(
x
,
y
),
(
z
,
)):
try
:
z
[
0
]
=
numpy
.
asarray
(
numpy
.
dot
(
x
,
y
))
except
ValueError
,
e
:
# The error raised by numpy has no shape information, we mean to add that
e
.
args
=
e
.
args
+
(
x
.
shape
,
y
.
shape
)
raise
def
__str__
(
self
):
return
"_dot22"
def
c_support_code
(
self
):
#return blas.cblas_header_text()
mod_str
=
"""
#ifndef MOD
#define MOD
%
#endif
"""
return
blas
.
blas_proto
()
+
mod_str
def
c_headers
(
self
):
return
[
'<iostream>'
]
def
c_libraries
(
self
):
return
blas
.
ldflags
()
def
c_code
(
self
,
node
,
name
,
(
_x
,
_y
),
(
_z
,
),
sub
):
return
"""
int unit = 0;
int type_num =
%(_x)
s->descr->type_num;
int type_size =
%(_x)
s->descr->elsize; // in bytes
npy_intp* Nx =
%(_x)
s->dimensions;
npy_intp* Ny =
%(_y)
s->dimensions;
npy_intp* Nz = 0; //
%(_z)
s->dimensions;
npy_intp* Sx =
%(_x)
s->strides;
npy_intp* Sy =
%(_y)
s->strides;
npy_intp* Sz = 0;//
%(_z)
s->strides;
//strides for x, y, z in dimensions 0, 1
int sx_0, sx_1, sy_0, sy_1, sz_0, sz_1;
if ((NULL ==
%(_z)
s)
|| (
%(_z)
s->dimensions[0] !=
%(_x)
s->dimensions[0])
|| (
%(_z)
s->dimensions[1] !=
%(_y)
s->dimensions[1]))
{
if (NULL !=
%(_z)
s) Py_XDECREF(
%(_z)
s);
npy_intp dims[2];
dims[0] =
%(_x)
s->dimensions[0];
dims[1] =
%(_y)
s->dimensions[1];
%(_z)
s = (PyArrayObject*)PyArray_SimpleNew(2, dims, type_num_
%(_x)
s);
if(!
%(_z)
s) {
PyErr_SetString(PyExc_MemoryError, "failed to alloc dot22 output");
%(fail)
s
}
}
Nz =
%(_z)
s->dimensions;
Sz =
%(_z)
s->strides;
if (
%(_x)
s->nd != 2) {PyErr_SetString(PyExc_NotImplementedError, "rank(x) != 2");
%(fail)
s;}
if (
%(_y)
s->nd != 2) {PyErr_SetString(PyExc_NotImplementedError, "rank(y) != 2");
%(fail)
s;}
if (
%(_z)
s->nd != 2) {PyErr_SetString(PyExc_NotImplementedError, "rank(z) != 2");
%(fail)
s;}
if ((
%(_x)
s->descr->type_num != PyArray_DOUBLE)
&& (
%(_x)
s->descr->type_num != PyArray_FLOAT))
{PyErr_SetString(PyExc_NotImplementedError, "type(x) is not double or float");
%(fail)
s;}
if ((
%(_y)
s->descr->type_num != PyArray_DOUBLE)
&& (
%(_y)
s->descr->type_num != PyArray_FLOAT))
{PyErr_SetString(PyExc_NotImplementedError, "type(y) is not double or float");
%(fail)
s;}
if ((
%(_z)
s->descr->type_num != PyArray_DOUBLE)
&& (
%(_z)
s->descr->type_num != PyArray_FLOAT))
{PyErr_SetString(PyExc_NotImplementedError, "type(z) is not double or float");
%(fail)
s;}
if ((
%(_x)
s->descr->type_num !=
%(_y)
s->descr->type_num)
||(
%(_x)
s->descr->type_num !=
%(_z)
s->descr->type_num))
{ PyErr_SetString(PyExc_NotImplementedError, "type(z), type(y), type(z) are not all the same");
%(fail)
s; }
if ((Nx[0] != Nz[0]) || (Nx[1] != Ny[0]) || (Ny[1] != Nz[1]))
{
PyErr_SetString(PyExc_ValueError, "Input dimensions do not agree");
%(fail)
s;
}
if ((Sx[0] < 1) || (Sx[1] < 1) || (Sx[0] MOD type_size) || (Sx[1] MOD type_size)
|| (Sy[0] < 1) || (Sy[1] < 1) || (Sy[0] MOD type_size) || (Sy[1] MOD type_size)
|| (Sz[0] < 1) || (Sz[1] < 1) || (Sz[0] MOD type_size) || (Sz[1] MOD type_size))
{
PyErr_SetString(PyExc_ValueError, "stride is not multiple of element size");
%(fail)
s;
}
/*
encode the stride structure of _x,_y,_z into a single integer
*/
unit |= ((Sx[1] == type_size) ? 0x0 : (Sx[0] == type_size) ? 0x1 : 0x2) << 8;
unit |= ((Sy[1] == type_size) ? 0x0 : (Sy[0] == type_size) ? 0x1 : 0x2) << 4;
unit |= ((Sz[1] == type_size) ? 0x0 : (Sz[0] == type_size) ? 0x1 : 0x2) << 0;
/* create appropriate strides for malformed matrices that are row or column
* vectors
*/
sx_0 = (Nx[0] > 1) ? Sx[0]/type_size : Nx[1];
sx_1 = (Nx[1] > 1) ? Sx[1]/type_size : Nx[0];
sy_0 = (Ny[0] > 1) ? Sy[0]/type_size : Ny[1];
sy_1 = (Ny[1] > 1) ? Sy[1]/type_size : Ny[0];
sz_0 = (Nz[0] > 1) ? Sz[0]/type_size : Nz[1];
sz_1 = (Nz[1] > 1) ? Sz[1]/type_size : Nz[0];
switch (type_num)
{
case PyArray_FLOAT:
{
float a = 1.0;
float b = 0.0;
float* x = (float*)PyArray_DATA(
%(_x)
s);
float* y = (float*)PyArray_DATA(
%(_y)
s);
float* z = (float*)PyArray_DATA(
%(_z)
s);
char N = 'N';
char T = 'T';
int Nz0 = Nz[0], Nz1 = Nz[1], Nx1 = Nx[1];
//std::cerr << (unit/256) MOD 16 << (unit / 16) MOD 16 << unit MOD 16<< '
\\
n';
switch(unit)
{
case 0x000: sgemm_(&N, &N, &Nz1, &Nz0, &Nx1, &a, y, &sy_0, x, &sx_0, &b, z, &sz_0); break;
case 0x100: sgemm_(&N, &T, &Nz1, &Nz0, &Nx1, &a, y, &sy_0, x, &sx_1, &b, z, &sz_0); break;
case 0x010: sgemm_(&T, &N, &Nz1, &Nz0, &Nx1, &a, y, &sy_1, x, &sx_0, &b, z, &sz_0); break;
case 0x110: sgemm_(&T, &T, &Nz1, &Nz0, &Nx1, &a, y, &sy_1, x, &sx_1, &b, z, &sz_0); break;
case 0x001: sgemm_(&T, &T, &Nz0, &Nz1, &Nx1, &a, x, &sx_0, y, &sy_0, &b, z, &sz_1); break;
case 0x101: sgemm_(&N, &T, &Nz0, &Nz1, &Nx1, &a, x, &sx_1, y, &sy_0, &b, z, &sz_1); break;
case 0x011: sgemm_(&T, &N, &Nz0, &Nz1, &Nx1, &a, x, &sx_0, y, &sy_1, &b, z, &sz_1); break;
case 0x111: sgemm_(&N, &N, &Nz0, &Nz1, &Nx1, &a, x, &sx_1, y, &sy_1, &b, z, &sz_1); break;
default: PyErr_SetString(PyExc_ValueError, "some matrix has no unit stride");
%(fail)
s;
};
#undef REAL
}
break;
case PyArray_DOUBLE:
{
double a = 1.0;
double b = 0.0;
double* x = (double*)PyArray_DATA(
%(_x)
s);
double* y = (double*)PyArray_DATA(
%(_y)
s);
double* z = (double*)PyArray_DATA(
%(_z)
s);
char N = 'N';
char T = 'T';
int Nz0 = Nz[0], Nz1 = Nz[1], Nx1 = Nx[1];
//std::cerr << (unit/256) MOD 16 << (unit / 16) MOD 16 << unit MOD 16<< '
\\
n';
switch(unit)
{
case 0x000: dgemm_(&N, &N, &Nz1, &Nz0, &Nx1, &a, y, &sy_0, x, &sx_0, &b, z, &sz_0); break;
case 0x100: dgemm_(&N, &T, &Nz1, &Nz0, &Nx1, &a, y, &sy_0, x, &sx_1, &b, z, &sz_0); break;
case 0x010: dgemm_(&T, &N, &Nz1, &Nz0, &Nx1, &a, y, &sy_1, x, &sx_0, &b, z, &sz_0); break;
case 0x110: dgemm_(&T, &T, &Nz1, &Nz0, &Nx1, &a, y, &sy_1, x, &sx_1, &b, z, &sz_0); break;
case 0x001: dgemm_(&T, &T, &Nz0, &Nz1, &Nx1, &a, x, &sx_0, y, &sy_0, &b, z, &sz_1); break;
case 0x101: dgemm_(&N, &T, &Nz0, &Nz1, &Nx1, &a, x, &sx_1, y, &sy_0, &b, z, &sz_1); break;
case 0x011: dgemm_(&T, &N, &Nz0, &Nz1, &Nx1, &a, x, &sx_0, y, &sy_1, &b, z, &sz_1); break;
case 0x111: dgemm_(&N, &N, &Nz0, &Nz1, &Nx1, &a, x, &sx_1, y, &sy_1, &b, z, &sz_1); break;
default: PyErr_SetString(PyExc_ValueError, "some matrix has no unit stride");
%(fail)
s;
};
#undef REAL
}
break;
}
"""
%
dict
(
locals
(),
**
sub
)
_dot22
=
_Dot22
()
@gof.local_optimizer
([
T
.
dot
])
def
local_dot_to_dot22
(
node
):
if
node
.
op
==
T
.
dot
:
return
[
_dot22
(
*
node
.
inputs
)]
else
:
return
False
T
.
opt
.
register_specialize
(
local_dot_to_dot22
)
@gof.local_optimizer
([
T
.
sub
])
def
local_sub_to_gemm
(
node
):
if
node
.
op
==
T
.
sub
:
subleft
,
subright
=
node
.
inputs
if
subright
.
owner
and
(
subright
.
owner
.
op
==
_dot22
):
dotleft
,
dotright
=
subright
.
owner
.
inputs
return
[
T
.
gemm
(
subleft
,
-
1.0
,
dotleft
,
dotright
,
1.0
)]
if
subright
.
owner
and
(
subright
.
owner
.
op
==
T
.
mul
):
mulleft
,
mulright
=
subright
.
owner
.
inputs
#TODO: we actually want to get any scalar here, not necessrily a constant
mulleft_const
=
opt
.
local_mul_canonizer
.
get_constant
(
mulleft
)
if
mulleft_const
is
not
None
:
assert
mulleft_const
.
size
()
==
1
mulleft_const
=
mulleft_const
.
flatten
()[
0
]
#subleft - (mulleft_const * ?)
if
mulright
.
owner
and
(
mulright
.
owner
.
op
==
T
.
add
):
#subleft - (mulleft_const * (? + ?))
addleft
,
addright
=
mulright
.
owner
.
inputs
if
addright
.
owner
and
addright
.
owner
.
op
==
T
.
DimShuffle
([
False
,
False
],
[
1
,
0
]):
#subleft - (mulleft_const * (? + ?.T))
raise
NotImplementedError
()
if
addright
.
owner
and
addright
.
owner
.
op
==
T
.
DimShuffle
([
False
,
False
],
[
1
,
0
],
inplace
=
True
):
#subleft - (mulleft_const * (? + ?.T))
transposed
=
addright
.
owner
.
inputs
[
0
]
if
transposed
.
owner
and
transposed
.
owner
.
op
==
_dot22
:
x
,
y
=
transposed
.
owner
.
inputs
#subleft - (mulleft_const * (addleft + dot(x, y).T))
if
addleft
.
owner
and
addleft
.
owner
.
op
==
_dot22
:
u
,
v
=
addleft
.
owner
.
inputs
#subleft - (mulleft_const * (dot(u,v) + dot(x, y).T))
return
[
T
.
gemm
(
T
.
gemm
(
subleft
,
-
mulleft_const
,
y
.
T
,
x
.
T
,
1.0
),
-
mulleft_const
,
u
,
v
,
1.0
)]
if
mulright
.
owner
and
(
mulright
.
owner
.
op
==
_dot22
):
dotleft
,
dotright
=
mulright
.
owner
.
inputs
#TODO: we actually want to get any scalar here, not necessrily a constant
mulleft_const
=
opt
.
local_mul_canonizer
.
get_constant
(
mulleft
)
if
mulleft_const
:
return
[
T
.
gemm
(
subleft
,
-
mulleft_const
,
dotleft
,
dotright
,
1.0
)]
if
mulleft
.
owner
and
(
mulleft
.
owner
.
op
==
_dot22
):
dotleft
,
dotright
=
mulleft
.
owner
.
inputs
#TODO: we actually want to get any scalar here, not necessrily a constant
mulright_const
=
opt
.
local_mul_canonizer
.
get_constant
(
mulright
)
if
mulright_const
:
return
[
T
.
gemm
(
subleft
,
-
mulright_const
,
dotleft
,
dotright
,
1.0
)]
return
False
T
.
opt
.
register_specialize
(
local_sub_to_gemm
)
if
0
:
class
Opt
(
object
):
merge
=
theano
.
gof
.
MergeOptimizer
()
...
...
theano/tensor/elemwise.py
浏览文件 @
f28d5cb9
...
...
@@ -205,10 +205,8 @@ class DimShuffle(Op):
if
self
.
inplace
:
print
"INPLACE"
get_base
=
[
'{ PyArrayObject * base =
%(input)
s'
,
'Py_INCREF((PyObject*)base)'
]
else
:
print
"NOT INPLACE"
get_base
=
[(
'{ PyArrayObject * base = (PyArrayObject*)PyArray_FromAny((PyObject*)
%(input)
s, NULL,'
'0, 0, NPY_ALIGNED|NPY_ENSURECOPY, NULL)'
)]
...
...
theano/tensor/opt.py
浏览文件 @
f28d5cb9
...
...
@@ -813,6 +813,248 @@ def constant_folding(node):
register_canonicalize
(
constant_folding
)
#################
# BLAS-related
#################
import
blas
class
_Dot22
(
gof
.
Op
):
"""Compute a matrix-matrix product.
This is a specialization of the more general Dot()
"""
def
make_node
(
self
,
x
,
y
):
assert
x
.
type
in
T
.
float_matrix_types
#makes sure x is a matrix
assert
y
.
type
==
x
.
type
#makes sure y is a matrix
bz
=
[
x
.
type
.
broadcastable
[
0
],
y
.
type
.
broadcastable
[
1
]]
outputs
=
[
T
.
tensor
(
x
.
type
.
dtype
,
bz
)]
return
gof
.
Apply
(
self
,
[
x
,
y
],
outputs
)
def
perform
(
self
,
node
,
(
x
,
y
),
(
z
,
)):
try
:
z
[
0
]
=
numpy
.
asarray
(
numpy
.
dot
(
x
,
y
))
except
ValueError
,
e
:
# The error raised by numpy has no shape information, we mean to add that
e
.
args
=
e
.
args
+
(
x
.
shape
,
y
.
shape
)
raise
def
__str__
(
self
):
return
"_dot22"
def
c_support_code
(
self
):
#return blas.cblas_header_text()
mod_str
=
"""
#ifndef MOD
#define MOD
%
#endif
"""
return
blas
.
blas_proto
()
+
mod_str
def
c_headers
(
self
):
return
[
'<iostream>'
]
def
c_libraries
(
self
):
return
blas
.
ldflags
()
def
c_code
(
self
,
node
,
name
,
(
_x
,
_y
),
(
_z
,
),
sub
):
return
"""
int unit = 0;
int type_num =
%(_x)
s->descr->type_num;
int type_size =
%(_x)
s->descr->elsize; // in bytes
npy_intp* Nx =
%(_x)
s->dimensions;
npy_intp* Ny =
%(_y)
s->dimensions;
npy_intp* Nz = 0; //
%(_z)
s->dimensions;
npy_intp* Sx =
%(_x)
s->strides;
npy_intp* Sy =
%(_y)
s->strides;
npy_intp* Sz = 0;//
%(_z)
s->strides;
//strides for x, y, z in dimensions 0, 1
int sx_0, sx_1, sy_0, sy_1, sz_0, sz_1;
if ((NULL ==
%(_z)
s)
|| (
%(_z)
s->dimensions[0] !=
%(_x)
s->dimensions[0])
|| (
%(_z)
s->dimensions[1] !=
%(_y)
s->dimensions[1]))
{
if (NULL !=
%(_z)
s) Py_XDECREF(
%(_z)
s);
npy_intp dims[2];
dims[0] =
%(_x)
s->dimensions[0];
dims[1] =
%(_y)
s->dimensions[1];
%(_z)
s = (PyArrayObject*)PyArray_SimpleNew(2, dims, type_num_
%(_x)
s);
if(!
%(_z)
s) {
PyErr_SetString(PyExc_MemoryError, "failed to alloc dot22 output");
%(fail)
s
}
}
Nz =
%(_z)
s->dimensions;
Sz =
%(_z)
s->strides;
if (
%(_x)
s->nd != 2) {PyErr_SetString(PyExc_NotImplementedError, "rank(x) != 2");
%(fail)
s;}
if (
%(_y)
s->nd != 2) {PyErr_SetString(PyExc_NotImplementedError, "rank(y) != 2");
%(fail)
s;}
if (
%(_z)
s->nd != 2) {PyErr_SetString(PyExc_NotImplementedError, "rank(z) != 2");
%(fail)
s;}
if ((
%(_x)
s->descr->type_num != PyArray_DOUBLE)
&& (
%(_x)
s->descr->type_num != PyArray_FLOAT))
{PyErr_SetString(PyExc_NotImplementedError, "type(x) is not double or float");
%(fail)
s;}
if ((
%(_y)
s->descr->type_num != PyArray_DOUBLE)
&& (
%(_y)
s->descr->type_num != PyArray_FLOAT))
{PyErr_SetString(PyExc_NotImplementedError, "type(y) is not double or float");
%(fail)
s;}
if ((
%(_z)
s->descr->type_num != PyArray_DOUBLE)
&& (
%(_z)
s->descr->type_num != PyArray_FLOAT))
{PyErr_SetString(PyExc_NotImplementedError, "type(z) is not double or float");
%(fail)
s;}
if ((
%(_x)
s->descr->type_num !=
%(_y)
s->descr->type_num)
||(
%(_x)
s->descr->type_num !=
%(_z)
s->descr->type_num))
{ PyErr_SetString(PyExc_NotImplementedError, "type(z), type(y), type(z) are not all the same");
%(fail)
s; }
if ((Nx[0] != Nz[0]) || (Nx[1] != Ny[0]) || (Ny[1] != Nz[1]))
{
PyErr_SetString(PyExc_ValueError, "Input dimensions do not agree");
%(fail)
s;
}
if ((Sx[0] < 1) || (Sx[1] < 1) || (Sx[0] MOD type_size) || (Sx[1] MOD type_size)
|| (Sy[0] < 1) || (Sy[1] < 1) || (Sy[0] MOD type_size) || (Sy[1] MOD type_size)
|| (Sz[0] < 1) || (Sz[1] < 1) || (Sz[0] MOD type_size) || (Sz[1] MOD type_size))
{
PyErr_SetString(PyExc_ValueError, "stride is not multiple of element size");
%(fail)
s;
}
/*
encode the stride structure of _x,_y,_z into a single integer
*/
unit |= ((Sx[1] == type_size) ? 0x0 : (Sx[0] == type_size) ? 0x1 : 0x2) << 8;
unit |= ((Sy[1] == type_size) ? 0x0 : (Sy[0] == type_size) ? 0x1 : 0x2) << 4;
unit |= ((Sz[1] == type_size) ? 0x0 : (Sz[0] == type_size) ? 0x1 : 0x2) << 0;
/* create appropriate strides for malformed matrices that are row or column
* vectors
*/
sx_0 = (Nx[0] > 1) ? Sx[0]/type_size : Nx[1];
sx_1 = (Nx[1] > 1) ? Sx[1]/type_size : Nx[0];
sy_0 = (Ny[0] > 1) ? Sy[0]/type_size : Ny[1];
sy_1 = (Ny[1] > 1) ? Sy[1]/type_size : Ny[0];
sz_0 = (Nz[0] > 1) ? Sz[0]/type_size : Nz[1];
sz_1 = (Nz[1] > 1) ? Sz[1]/type_size : Nz[0];
switch (type_num)
{
case PyArray_FLOAT:
{
float a = 1.0;
float b = 0.0;
float* x = (float*)PyArray_DATA(
%(_x)
s);
float* y = (float*)PyArray_DATA(
%(_y)
s);
float* z = (float*)PyArray_DATA(
%(_z)
s);
char N = 'N';
char T = 'T';
int Nz0 = Nz[0], Nz1 = Nz[1], Nx1 = Nx[1];
//std::cerr << (unit/256) MOD 16 << (unit / 16) MOD 16 << unit MOD 16<< '
\\
n';
switch(unit)
{
case 0x000: sgemm_(&N, &N, &Nz1, &Nz0, &Nx1, &a, y, &sy_0, x, &sx_0, &b, z, &sz_0); break;
case 0x100: sgemm_(&N, &T, &Nz1, &Nz0, &Nx1, &a, y, &sy_0, x, &sx_1, &b, z, &sz_0); break;
case 0x010: sgemm_(&T, &N, &Nz1, &Nz0, &Nx1, &a, y, &sy_1, x, &sx_0, &b, z, &sz_0); break;
case 0x110: sgemm_(&T, &T, &Nz1, &Nz0, &Nx1, &a, y, &sy_1, x, &sx_1, &b, z, &sz_0); break;
case 0x001: sgemm_(&T, &T, &Nz0, &Nz1, &Nx1, &a, x, &sx_0, y, &sy_0, &b, z, &sz_1); break;
case 0x101: sgemm_(&N, &T, &Nz0, &Nz1, &Nx1, &a, x, &sx_1, y, &sy_0, &b, z, &sz_1); break;
case 0x011: sgemm_(&T, &N, &Nz0, &Nz1, &Nx1, &a, x, &sx_0, y, &sy_1, &b, z, &sz_1); break;
case 0x111: sgemm_(&N, &N, &Nz0, &Nz1, &Nx1, &a, x, &sx_1, y, &sy_1, &b, z, &sz_1); break;
default: PyErr_SetString(PyExc_ValueError, "some matrix has no unit stride");
%(fail)
s;
};
#undef REAL
}
break;
case PyArray_DOUBLE:
{
double a = 1.0;
double b = 0.0;
double* x = (double*)PyArray_DATA(
%(_x)
s);
double* y = (double*)PyArray_DATA(
%(_y)
s);
double* z = (double*)PyArray_DATA(
%(_z)
s);
char N = 'N';
char T = 'T';
int Nz0 = Nz[0], Nz1 = Nz[1], Nx1 = Nx[1];
//std::cerr << (unit/256) MOD 16 << (unit / 16) MOD 16 << unit MOD 16<< '
\\
n';
switch(unit)
{
case 0x000: dgemm_(&N, &N, &Nz1, &Nz0, &Nx1, &a, y, &sy_0, x, &sx_0, &b, z, &sz_0); break;
case 0x100: dgemm_(&N, &T, &Nz1, &Nz0, &Nx1, &a, y, &sy_0, x, &sx_1, &b, z, &sz_0); break;
case 0x010: dgemm_(&T, &N, &Nz1, &Nz0, &Nx1, &a, y, &sy_1, x, &sx_0, &b, z, &sz_0); break;
case 0x110: dgemm_(&T, &T, &Nz1, &Nz0, &Nx1, &a, y, &sy_1, x, &sx_1, &b, z, &sz_0); break;
case 0x001: dgemm_(&T, &T, &Nz0, &Nz1, &Nx1, &a, x, &sx_0, y, &sy_0, &b, z, &sz_1); break;
case 0x101: dgemm_(&N, &T, &Nz0, &Nz1, &Nx1, &a, x, &sx_1, y, &sy_0, &b, z, &sz_1); break;
case 0x011: dgemm_(&T, &N, &Nz0, &Nz1, &Nx1, &a, x, &sx_0, y, &sy_1, &b, z, &sz_1); break;
case 0x111: dgemm_(&N, &N, &Nz0, &Nz1, &Nx1, &a, x, &sx_1, y, &sy_1, &b, z, &sz_1); break;
default: PyErr_SetString(PyExc_ValueError, "some matrix has no unit stride");
%(fail)
s;
};
#undef REAL
}
break;
}
"""
%
dict
(
locals
(),
**
sub
)
_dot22
=
_Dot22
()
@gof.local_optimizer
([
T
.
dot
])
def
local_dot_to_dot22
(
node
):
if
node
.
op
==
T
.
dot
:
return
[
_dot22
(
*
node
.
inputs
)]
else
:
return
False
register_specialize
(
local_dot_to_dot22
)
@gof.local_optimizer
([
T
.
sub
])
def
local_sub_to_gemm
(
node
):
"""This is a massive beast for recognizing all the ways that a subtraction could be
replaced by a GEMM
"""
if
node
.
op
==
T
.
sub
:
subleft
,
subright
=
node
.
inputs
#EXPRESSION: subleft - subright
if
subright
.
owner
and
(
subright
.
owner
.
op
==
_dot22
):
dotleft
,
dotright
=
subright
.
owner
.
inputs
return
[
T
.
gemm
(
subleft
,
-
1.0
,
dotleft
,
dotright
,
1.0
)]
if
subright
.
owner
and
(
subright
.
owner
.
op
==
T
.
mul
):
mulleft
,
mulright
=
subright
.
owner
.
inputs
#EXPRESSION: subleft - (mulleft * mulright)
#TODO: we actually want to get any scalar here, not necessrily a constant
mulleft_const
=
local_mul_canonizer
.
get_constant
(
mulleft
)
if
mulleft_const
is
not
None
and
mulleft_const
.
size
==
1
:
mulleft_const
=
mulleft_const
.
flatten
()[
0
]
#EXPRESSION: subleft - (mulleft_const * ?)
if
mulright
.
owner
and
(
mulright
.
owner
.
op
==
T
.
add
):
#EXPRESSION: subleft - (mulleft_const * (? + ?))
addleft
,
addright
=
mulright
.
owner
.
inputs
if
addright
.
owner
and
addright
.
owner
.
op
==
T
.
DimShuffle
([
False
,
False
],
[
1
,
0
]):
#EXPRESSION: subleft - (mulleft_const * (? + ?.T))
raise
NotImplementedError
()
if
addright
.
owner
and
addright
.
owner
.
op
==
T
.
DimShuffle
([
False
,
False
],
[
1
,
0
],
inplace
=
True
):
#EXPRESSION: subleft - (mulleft_const * (? + ?.T))
transposed
=
addright
.
owner
.
inputs
[
0
]
if
transposed
.
owner
and
transposed
.
owner
.
op
==
_dot22
:
x
,
y
=
transposed
.
owner
.
inputs
#EXPRESSION: subleft - (mulleft_const * (addleft + dot(x, y).T))
if
addleft
.
owner
and
addleft
.
owner
.
op
==
_dot22
:
u
,
v
=
addleft
.
owner
.
inputs
#EXPRESSION: subleft - (mulleft_const * (dot(u,v) + dot(x, y).T))
return
[
T
.
gemm
(
T
.
gemm
(
subleft
,
-
mulleft_const
,
y
.
T
,
x
.
T
,
1.0
),
-
mulleft_const
,
u
,
v
,
1.0
)]
if
mulright
.
owner
and
(
mulright
.
owner
.
op
==
_dot22
):
dotleft
,
dotright
=
mulright
.
owner
.
inputs
#EXPRESSION: subleft - (mulleft_const * dot(dotleft, dotright))
return
[
T
.
gemm
(
subleft
,
-
mulleft_const
,
dotleft
,
dotright
,
1.0
)]
mulright_const
=
local_mul_canonizer
.
get_constant
(
mulright
)
if
mulright_const
is
not
None
and
mulright_const
.
size
==
1
:
mulright_const
=
mulright_const
.
flatten
()[
0
]
#EXPRESSION: subleft - (? * mulright_const)
if
mulleft
.
owner
and
(
mulleft
.
owner
.
op
==
_dot22
):
dotleft
,
dotright
=
mulleft
.
owner
.
inputs
#EXPRESSION: subleft - (dot(dotleft, dotright) * mulright_const)
return
[
T
.
gemm
(
subleft
,
-
mulright_const
,
dotleft
,
dotright
,
1.0
)]
return
False
register_specialize
(
local_sub_to_gemm
)
...
...
编写
预览
Markdown
格式
0%
重试
或
添加新文件
添加附件
取消
您添加了
0
人
到此讨论。请谨慎行事。
请先完成此评论的编辑!
取消
请
注册
或者
登录
后发表评论