Skip to content
This repository was archived by the owner on Jan 26, 2026. It is now read-only.

Commit 015829d

Browse files
fschlimbfschlimb
committed
fixing empty(), getting closer to working with PRK stencil
1 parent c60c534 commit 015829d

File tree

7 files changed

+244
-51
lines changed

7 files changed

+244
-51
lines changed

‎ddptensor/ddptensor.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -55,4 +55,4 @@ def __getitem__(self, key):
5555
def __setitem__(self, key, value):
5656
key = key if isinstance(key, tuple) else (key,)
5757
key = [x if isinstance(x, slice) else slice(x, x + 1, 1) for x in key]
58-
self._t.__setitem__(key, value._t) # if isinstance(value, dtensor) else value)
58+
self._t.__setitem__(key, value._t if isinstance(value, dtensor) else value)

‎src/Creator.cpp‎

Lines changed: 6 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -127,12 +127,15 @@ struct DeferredFull : public Deferred {
127127
const PyScalar &val, ::imex::ptensor::DType &dtyp) {
128128
dtyp = jit::PT_DTYPE<T>::value;
129129

130-
if constexpr (std::is_floating_point_v<T>)
130+
if (is_none(val)) {
131+
return {};
132+
} else if constexpr (std::is_floating_point_v<T>) {
131133
return ::imex::createFloat<sizeof(T) * 8>(loc, builder, val._float);
132-
else if constexpr (std::is_same_v<bool, T>)
134+
} else if constexpr (std::is_same_v<bool, T>) {
133135
return ::imex::createInt<1>(loc, builder, val._int);
134-
else if constexpr (std::is_integral_v<T>)
136+
} else if constexpr (std::is_integral_v<T>) {
135137
return ::imex::createInt<sizeof(T) * 8>(loc, builder, val._int);
138+
}
136139
assert("Unsupported dtype in dispatch");
137140
return {};
138141
};

‎src/SetGetItem.cpp‎

Lines changed: 8 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@
66
*/
77

88
#include "ddptensor/SetGetItem.hpp"
9+
#include "ddptensor/Creator.hpp"
910
#include "ddptensor/DDPTensorImpl.hpp"
1011
#include "ddptensor/Factory.hpp"
1112
#include "ddptensor/Mediator.hpp"
@@ -258,8 +259,13 @@ struct DeferredSetItem : public Deferred {
258259
};
259260

260261
ddptensor *SetItem::__setitem__(ddptensor &a, const std::vector<py::slice> &v,
261-
const ddptensor &b) {
262-
return new ddptensor(defer<DeferredSetItem>(a.get(), b.get(), v));
262+
const py::object &b) {
263+
264+
auto bb = Creator::mk_future(b);
265+
auto res = new ddptensor(defer<DeferredSetItem>(a.get(), bb.first->get(), v));
266+
if (bb.second)
267+
delete bb.first;
268+
return res;
263269
}
264270

265271
struct DeferredGetItem : public Deferred {

‎src/include/ddptensor/CppTypes.hpp‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@
88

99
#include "p2c_ids.hpp"
1010

11+
#include <cmath>
1112
#include <cstring>
1213
#include <numeric>
1314
#include <vector>
@@ -30,6 +31,8 @@ union PyScalar {
3031
double _float;
3132
};
3233

34+
inline bool is_none(PyScalar s) { return std::isnan(s._float); }
35+
3336
enum _RANKS : rank_type {
3437
NOOWNER = std::numeric_limits<rank_type>::max(),
3538
REPLICATED = std::numeric_limits<rank_type>::max() - 1,

‎src/include/ddptensor/PyTypes.hpp‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -15,6 +15,10 @@ template <typename T> py::object get_impl_dtype() {
1515

1616
inline PyScalar mk_scalar(const py::object &b, DTypeId dtype) {
1717
PyScalar s;
18+
if (b.is_none()) {
19+
s._float = std::numeric_limits<double>::quiet_NaN();
20+
return s;
21+
}
1822
switch (dtype) {
1923
case FLOAT64:
2024
case FLOAT32:

‎src/include/ddptensor/SetGetItem.hpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -23,5 +23,5 @@ struct GetItem {
2323

2424
struct SetItem {
2525
static ddptensor *__setitem__(ddptensor &a, const std::vector<py::slice> &v,
26-
const ddptensor &b);
26+
const py::object &b);
2727
};

‎test/stencil-2d.py‎

Lines changed: 221 additions & 44 deletions
Original file line numberDiff line numberDiff line change
@@ -1,44 +1,221 @@
1-
import ddptensor as aa
2-
3-
# import numpy as aa
4-
import numpy as np
5-
6-
aa.init(False)
7-
8-
n = 16
9-
r = 2
10-
A = aa.ones((n, n), dtype=aa.int64)
11-
B = aa.zeros((n, n), dtype=aa.int64)
12-
W = aa.ones(((2 * r + 1), (2 * r + 1)), dtype=aa.int64)
13-
for i in range(3):
14-
# if i:
15-
# C = B[2:n-2,2:n-2]
16-
# B[2:n-2,2:n-2] = C
17-
# else:
18-
# D = B[2:n-2,2:n-2]
19-
B[2 : n - 2, 2 : n - 2] = (
20-
B[2 : n - 2, 2 : n - 2]
21-
+ W[2, 2] * A[2 : n - 2, 2 : n - 2]
22-
+ W[2, 0] * A[2 : n - 2, 0 : n - 4]
23-
+ W[2, 1] * A[2 : n - 2, 1 : n - 3]
24-
+ W[2, 3] * A[2 : n - 2, 3 : n - 1]
25-
+ W[2, 4] * A[2 : n - 2, 4 : n - 0]
26-
+ W[0, 2] * A[0 : n - 4, 2 : n - 2]
27-
+ W[1, 2] * A[1 : n - 3, 2 : n - 2]
28-
+ W[3, 2] * A[3 : n - 1, 2 : n - 2]
29-
+ W[4, 2] * A[4 : n - 0, 2 : n - 2]
30-
)
31-
# A[2:n-2,2:n-2] \
32-
# + A[2:n-2,2:n-2] \
33-
# + A[2:n-2,0:n-4] \
34-
# + A[2:n-2,1:n-3] \
35-
# + A[2:n-2,3:n-1] \
36-
# + A[2:n-2,4:n-0] \
37-
# + A[0:n-4,2:n-2] \
38-
# + A[1:n-3,2:n-2] \
39-
# + A[3:n-1,2:n-2] \
40-
# + A[4:n-0,2:n-2]
41-
A[0:n, 0:n] = A + 1
42-
print(B)
43-
44-
aa.fini()
1+
#!/usr/bin/env python3
2+
#
3+
# Copyright (c) 2015, Intel Corporation
4+
#
5+
# Redistribution and use in source and binary forms, with or without
6+
# modification, are permitted provided that the following conditions
7+
# are met:
8+
#
9+
# * Redistributions of source code must retain the above copyright
10+
# notice, this list of conditions and the following disclaimer.
11+
# * Redistributions in binary form must reproduce the above
12+
# copyright notice, this list of conditions and the following
13+
# disclaimer in the documentation and/or other materials provided
14+
# with the distribution.
15+
# * Neither the name of Intel Corporation nor the names of its
16+
# contributors may be used to endorse or promote products
17+
# derived from this software without specific prior written
18+
# permission.
19+
#
20+
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21+
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22+
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
23+
# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
24+
# COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
25+
# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
26+
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
27+
# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
28+
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
29+
# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
30+
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31+
# POSSIBILITY OF SUCH DAMAGE.
32+
#
33+
#
34+
# *******************************************************************
35+
#
36+
# NAME: Stencil
37+
#
38+
# PURPOSE: This program tests the efficiency with which a space-invariant,
39+
# linear, symmetric filter (stencil) can be applied to a square
40+
# grid or image.
41+
#
42+
# USAGE: The program takes as input the linear
43+
# dimension of the grid, and the number of iterations on the grid
44+
#
45+
# <progname> <iterations> <grid size>
46+
#
47+
# The output consists of diagnostics to make sure the
48+
# algorithm worked, and of timing statistics.
49+
#
50+
# HISTORY: - Written by Rob Van der Wijngnumpyrt, February 2009.
51+
# - RvdW: Removed unrolling pragmas for clarity;
52+
# added constant to array "in" at end of each iteration to force
53+
# refreshing of neighbor data in parallel versions; August 2013
54+
# - Converted to Python by Jeff Hammond, February 2016.
55+
#
56+
# *******************************************************************
57+
58+
import sys
59+
60+
print(
61+
"Python version = ", str(sys.version_info.major) + "." + str(sys.version_info.minor)
62+
)
63+
if sys.version_info >= (3, 3):
64+
from time import process_time as timer
65+
else:
66+
from timeit import default_timer as timer
67+
68+
import ddptensor as numpy
69+
70+
# print('Numpy version = ', numpy.version.version)
71+
72+
73+
def main():
74+
# ********************************************************************
75+
# read and test input parameters
76+
# ********************************************************************
77+
78+
print("Parallel Research Kernels")
79+
print("Python stencil execution on 2D grid")
80+
81+
if len(sys.argv) < 3:
82+
print("argument count = ", len(sys.argv))
83+
sys.exit(
84+
"Usage: ./stencil <# iterations> <array dimension> [<star/stencil> <radius>]"
85+
)
86+
87+
iterations = int(sys.argv[1])
88+
if iterations < 1:
89+
sys.exit("ERROR: iterations must be >= 1")
90+
91+
n = int(sys.argv[2])
92+
if n < 1:
93+
sys.exit("ERROR: array dimension must be >= 1")
94+
95+
if len(sys.argv) > 3:
96+
pattern = sys.argv[3]
97+
else:
98+
pattern = "star"
99+
100+
if len(sys.argv) > 4:
101+
r = int(sys.argv[4])
102+
if r < 1:
103+
sys.exit("ERROR: Stencil radius should be positive")
104+
if (2 * r + 1) > n:
105+
sys.exit("ERROR: Stencil radius exceeds grid size")
106+
else:
107+
r = 2
108+
109+
print("Number of iterations = ", iterations)
110+
print("Grid size = ", n)
111+
print("Radius of stencil = ", r)
112+
if pattern == "star":
113+
print("Type of stencil = ", "star")
114+
else:
115+
print("Type of stencil = ", "stencil")
116+
print("Data type = double precision")
117+
print("Compact representation of stencil loop body")
118+
119+
# there is certainly a more Pythonic way to initialize W,
120+
# but it will have no impact on performance.
121+
W = numpy.zeros(((2 * r + 1), (2 * r + 1)), dtype=numpy.float64)
122+
if pattern == "star":
123+
stencil_size = 4 * r + 1
124+
for i in range(1, r + 1):
125+
W[r, r + i] = +1.0 / (2 * i * r)
126+
W[r + i, r] = +1.0 / (2 * i * r)
127+
W[r, r - i] = -1.0 / (2 * i * r)
128+
W[r - i, r] = -1.0 / (2 * i * r)
129+
130+
else:
131+
stencil_size = (2 * r + 1) ** 2.0
132+
for j in range(1, r + 1):
133+
for i in range(-j + 1, j):
134+
W[r + i, r + j] = +1.0 / (4 * j * (2 * j - 1) * r)
135+
W[r + i, r - j] = -1.0 / (4 * j * (2 * j - 1) * r)
136+
W[r + j, r + i] = +1.0 / (4 * j * (2 * j - 1) * r)
137+
W[r - j, r + i] = -1.0 / (4 * j * (2 * j - 1) * r)
138+
139+
W[r + j, r + j] = +1.0 / (4 * j * r)
140+
W[r - j, r - j] = -1.0 / (4 * j * r)
141+
142+
# A = numpy.fromfunction(lambda i,j: i+j, (n,n), dtype=float)
143+
A = numpy.empty((n, n), dtype=numpy.float64)
144+
for i in range(n):
145+
for j in range(n):
146+
A[i, j] = float(i + j)
147+
print(A.dtype)
148+
B = numpy.zeros((n, n), dtype=numpy.float64)
149+
150+
for k in range(iterations + 1):
151+
# start timer after a warmup iteration
152+
if k < 1:
153+
t0 = timer()
154+
155+
if pattern == "star":
156+
if r == 2:
157+
B[2 : n - 2, 2 : n - 2] = (
158+
B[2 : n - 2, 2 : n - 2]
159+
+ W[2, 2] * A[2 : n - 2, 2 : n - 2]
160+
+ W[2, 0] * A[2 : n - 2, 0 : n - 4]
161+
+ W[2, 1] * A[2 : n - 2, 1 : n - 3]
162+
+ W[2, 3] * A[2 : n - 2, 3 : n - 1]
163+
+ W[2, 4] * A[2 : n - 2, 4 : n - 0]
164+
+ W[0, 2] * A[0 : n - 4, 2 : n - 2]
165+
+ W[1, 2] * A[1 : n - 3, 2 : n - 2]
166+
+ W[3, 2] * A[3 : n - 1, 2 : n - 2]
167+
+ W[4, 2] * A[4 : n - 0, 2 : n - 2]
168+
)
169+
else:
170+
b = n - r
171+
B[r:b, r:b] = B[r:b, r:b] + W[r, r] * A[r:b, r:b]
172+
for s in range(1, r + 1):
173+
B[r:b, r:b] = (
174+
B[r:b, r:b]
175+
+ W[r, r - s] * A[r:b, r - s : b - s]
176+
+ W[r, r + s] * A[r:b, r + s : b + s]
177+
+ W[r - s, r] * A[r - s : b - s, r:b]
178+
+ W[r + s, r] * A[r + s : b + s, r:b]
179+
)
180+
else: # stencil
181+
if r > 0:
182+
b = n - r
183+
for s in range(-r, r + 1):
184+
for t in range(-r, r + 1):
185+
B[r:b, r:b] = (
186+
B[r:b, r:b]
187+
+ W[r + t, r + s] * A[r + t : b + t, r + s : b + s]
188+
)
189+
190+
A = A + 1.0
191+
192+
t1 = timer()
193+
stencil_time = t1 - t0
194+
195+
# ******************************************************************************
196+
# * Analyze and output results.
197+
# ******************************************************************************
198+
199+
print(W, B)
200+
# norm = numpy.linalg.norm(numpy.reshape(B,n*n),ord=1)
201+
# active_points = (n-2*r)**2
202+
# norm /= active_points
203+
204+
# epsilon=1.e-8
205+
206+
# # verify correctness
207+
# reference_norm = 2*(iterations+1)
208+
# if abs(norm-reference_norm) < epsilon:
209+
# print('Solution validates')
210+
# flops = (2*stencil_size+1) * active_points
211+
# avgtime = stencil_time/iterations
212+
# print('Rate (MFlops/s): ',1.e-6*flops/avgtime, ' Avg time (s): ',avgtime)
213+
# else:
214+
# print('ERROR: L1 norm = ', norm,' Reference L1 norm = ', reference_norm)
215+
# sys.exit()
216+
217+
218+
if __name__ == "__main__":
219+
numpy.init(False)
220+
main()
221+
numpy.fini()

0 commit comments

Comments
 (0)