11# Copyright 2019-2023 ETH Zurich and the DaCe authors. All rights reserved.
2+ import ctypes
23import dace
34import numpy as np
45import pytest
@@ -18,7 +19,7 @@ def csr_to_dense_python(A: CSR, B: dace.float32[M, N]):
1819 for i in dace .map [0 :M ]:
1920 for idx in dace .map [A .indptr [i ]:A .indptr [i + 1 ]]:
2021 B [i , A .indices [idx ]] = A .data [idx ]
21-
22+
2223 rng = np .random .default_rng (42 )
2324 A = sparse .random (20 , 20 , density = 0.1 , format = 'csr' , dtype = np .float32 , random_state = rng )
2425 B = np .zeros ((20 , 20 ), dtype = np .float32 )
@@ -41,7 +42,7 @@ def test_write_structure():
4142 M , N , nnz = (dace .symbol (s ) for s in ('M' , 'N' , 'nnz' ))
4243 CSR = dace .data .Structure (dict (indptr = dace .int32 [M + 1 ], indices = dace .int32 [nnz ], data = dace .float32 [nnz ]),
4344 name = 'CSRMatrix' )
44-
45+
4546 @dace .program
4647 def dense_to_csr_python (A : dace .float32 [M , N ], B : CSR ):
4748 idx = 0
@@ -53,7 +54,7 @@ def dense_to_csr_python(A: dace.float32[M, N], B: CSR):
5354 B .indices [idx ] = j
5455 idx += 1
5556 B .indptr [M ] = idx
56-
57+
5758 rng = np .random .default_rng (42 )
5859 tmp = sparse .random (20 , 20 , density = 0.1 , format = 'csr' , dtype = np .float32 , random_state = rng )
5960 A = tmp .toarray ()
@@ -75,7 +76,7 @@ def test_local_structure():
7576 M , N , nnz = (dace .symbol (s ) for s in ('M' , 'N' , 'nnz' ))
7677 CSR = dace .data .Structure (dict (indptr = dace .int32 [M + 1 ], indices = dace .int32 [nnz ], data = dace .float32 [nnz ]),
7778 name = 'CSRMatrix' )
78-
79+
7980 @dace .program
8081 def dense_to_csr_local_python (A : dace .float32 [M , N ], B : CSR ):
8182 tmp = dace .define_local_structure (CSR )
@@ -91,7 +92,7 @@ def dense_to_csr_local_python(A: dace.float32[M, N], B: CSR):
9192 B .indptr [:] = tmp .indptr [:]
9293 B .indices [:] = tmp .indices [:]
9394 B .data [:] = tmp .data [:]
94-
95+
9596 rng = np .random .default_rng (42 )
9697 tmp = sparse .random (20 , 20 , density = 0.1 , format = 'csr' , dtype = np .float32 , random_state = rng )
9798 A = tmp .toarray ()
@@ -118,12 +119,11 @@ def __init__(self, diag, upper, lower):
118119 self .lower = lower
119120
120121 n , nblocks = dace .symbol ('n' ), dace .symbol ('nblocks' )
121- BlockTriDiagonal = dace .data .Structure (
122- dict (diagonal = dace .complex128 [nblocks , n , n ],
123- upper = dace .complex128 [nblocks , n , n ],
124- lower = dace .complex128 [nblocks , n , n ]),
125- name = 'BlockTriDiagonalMatrix' )
126-
122+ BlockTriDiagonal = dace .data .Structure (dict (diagonal = dace .complex128 [nblocks , n , n ],
123+ upper = dace .complex128 [nblocks , n , n ],
124+ lower = dace .complex128 [nblocks , n , n ]),
125+ name = 'BlockTriDiagonalMatrix' )
126+
127127 @dace .program
128128 def rgf_leftToRight (A : BlockTriDiagonal , B : BlockTriDiagonal , n_ : dace .int32 , nblocks_ : dace .int32 ):
129129
@@ -139,42 +139,41 @@ def rgf_leftToRight(A: BlockTriDiagonal, B: BlockTriDiagonal, n_: dace.int32, nb
139139 # 2. Forward substitution
140140 # From left to right
141141 for i in range (1 , nblocks_ ):
142- tmp [i ] = np .linalg .inv (A .diagonal [i ] - A .lower [i - 1 ] @ tmp [i - 1 ] @ A .upper [i - 1 ])
142+ tmp [i ] = np .linalg .inv (A .diagonal [i ] - A .lower [i - 1 ] @ tmp [i - 1 ] @ A .upper [i - 1 ])
143143 # 3. Initialisation of last element of B
144144 B .diagonal [- 1 ] = tmp [- 1 ]
145145
146146 # 4. Backward substitution
147147 # From right to left
148148
149- for i in range (nblocks_ - 2 , - 1 , - 1 ):
150- B .diagonal [i ] = tmp [i ] @ (identity + A .upper [i ] @ B .diagonal [i + 1 ] @ A .lower [i ] @ tmp [i ])
151- B .upper [i ] = - tmp [i ] @ A .upper [i ] @ B .diagonal [i + 1 ]
152- B .lower [i ] = np .transpose (B .upper [i ])
153-
149+ for i in range (nblocks_ - 2 , - 1 , - 1 ):
150+ B .diagonal [i ] = tmp [i ] @ (identity + A .upper [i ] @ B .diagonal [i + 1 ] @ A .lower [i ] @ tmp [i ])
151+ B .upper [i ] = - tmp [i ] @ A .upper [i ] @ B .diagonal [i + 1 ]
152+ B .lower [i ] = np .transpose (B .upper [i ])
153+
154154 rng = np .random .default_rng (42 )
155155
156156 A_diag = rng .random ((10 , 20 , 20 )) + 1j * rng .random ((10 , 20 , 20 ))
157157 A_upper = rng .random ((10 , 20 , 20 )) + 1j * rng .random ((10 , 20 , 20 ))
158- A_lower = rng .random ((10 , 20 , 20 )) + 1j * rng .random ((10 , 20 , 20 ))
158+ A_lower = rng .random ((10 , 20 , 20 )) + 1j * rng .random ((10 , 20 , 20 ))
159159 inpBTD = BlockTriDiagonal .dtype ._typeclass .as_ctypes ()(diagonal = A_diag .__array_interface__ ['data' ][0 ],
160160 upper = A_upper .__array_interface__ ['data' ][0 ],
161161 lower = A_lower .__array_interface__ ['data' ][0 ])
162-
162+
163163 B_diag = np .zeros ((10 , 20 , 20 ), dtype = np .complex128 )
164164 B_upper = np .zeros ((10 , 20 , 20 ), dtype = np .complex128 )
165165 B_lower = np .zeros ((10 , 20 , 20 ), dtype = np .complex128 )
166166 outBTD = BlockTriDiagonal .dtype ._typeclass .as_ctypes ()(diagonal = B_diag .__array_interface__ ['data' ][0 ],
167167 upper = B_upper .__array_interface__ ['data' ][0 ],
168168 lower = B_lower .__array_interface__ ['data' ][0 ])
169-
169+
170170 func = rgf_leftToRight .compile ()
171171 func (A = inpBTD , B = outBTD , n_ = A_diag .shape [1 ], nblocks_ = A_diag .shape [0 ], n = A_diag .shape [1 ], nblocks = A_diag .shape [0 ])
172172
173173 A = BTD (A_diag , A_upper , A_lower )
174- B = BTD (np .zeros ((10 , 20 , 20 ), dtype = np .complex128 ),
175- np .zeros ((10 , 20 , 20 ), dtype = np .complex128 ),
174+ B = BTD (np .zeros ((10 , 20 , 20 ), dtype = np .complex128 ), np .zeros ((10 , 20 , 20 ), dtype = np .complex128 ),
176175 np .zeros ((10 , 20 , 20 ), dtype = np .complex128 ))
177-
176+
178177 rgf_leftToRight .f (A , B , A_diag .shape [1 ], A_diag .shape [0 ])
179178
180179 assert np .allclose (B .diagonal , B_diag )
@@ -195,15 +194,15 @@ def csr_to_dense_python(A: CSR, B: dace.float32[M, N]):
195194 for i in dace .map [0 :M ]:
196195 for idx in dace .map [A .indptr [i ]:A .indptr [i + 1 ]]:
197196 B [i , A .indices [idx ]] = A .data [idx ]
198-
197+
199198 rng = np .random .default_rng (42 )
200199 A = sparse .random (20 , 20 , density = 0.1 , format = 'csr' , dtype = np .float32 , random_state = rng )
201200 ref = A .toarray ()
202201
203202 inpA = CSR .dtype ._typeclass .as_ctypes ()(indptr = A .indptr .__array_interface__ ['data' ][0 ],
204203 indices = A .indices .__array_interface__ ['data' ][0 ],
205204 data = A .data .__array_interface__ ['data' ][0 ])
206-
205+
207206 # TODO: The following doesn't work because we need to create a Structure data descriptor from the ctypes class.
208207 # csr_to_dense_python(inpA, B)
209208 naive = csr_to_dense_python .to_sdfg (simplify = False )
@@ -224,9 +223,99 @@ def csr_to_dense_python(A: CSR, B: dace.float32[M, N]):
224223 assert np .allclose (B , ref )
225224
226225
226+ def test_write_structure_in_map ():
227+ M = dace .symbol ('M' )
228+ N = dace .symbol ('N' )
229+ Bundle = dace .data .Structure (members = {
230+ "data" : dace .data .Array (dace .float32 , (M , N )),
231+ "size" : dace .data .Scalar (dace .int64 )
232+ },
233+ name = "BundleType" )
234+
235+ @dace .program
236+ def init_prog (bundle : Bundle , fill_value : int ) -> None :
237+ for index in dace .map [0 :bundle .size ]:
238+ bundle .data [index , :] = fill_value
239+
240+ data = np .zeros ((10 , 5 ), dtype = np .float32 )
241+ fill_value = 42
242+ inp_struct = Bundle .dtype .base_type .as_ctypes ()(
243+ data = data .__array_interface__ ['data' ][0 ],
244+ size = 9 ,
245+ )
246+ ref = np .zeros ((10 , 5 ), dtype = np .float32 )
247+ ref [:9 , :] = fill_value
248+
249+ init_prog .compile ()(inp_struct , fill_value , M = 10 , N = 5 )
250+
251+ assert np .allclose (data , ref )
252+
253+
254+ def test_readwrite_structure_in_map ():
255+ M = dace .symbol ('M' )
256+ N = dace .symbol ('N' )
257+ Bundle = dace .data .Structure (members = {
258+ "data" : dace .data .Array (dace .float32 , (M , N )),
259+ "data2" : dace .data .Array (dace .float32 , (M , N )),
260+ "size" : dace .data .Scalar (dace .int64 )
261+ },
262+ name = "BundleTypeTwoArrays" )
263+
264+ @dace .program
265+ def copy_prog (bundle : Bundle ) -> None :
266+ for index in dace .map [0 :bundle .size ]:
267+ bundle .data [index , :] = bundle .data2 [index , :] + 5
268+
269+ data = np .zeros ((10 , 5 ), dtype = np .float32 )
270+ data2 = np .ones ((10 , 5 ), dtype = np .float32 )
271+ inp_struct = Bundle .dtype .base_type .as_ctypes ()(
272+ data = data .__array_interface__ ['data' ][0 ],
273+ data2 = data2 .__array_interface__ ['data' ][0 ],
274+ size = ctypes .c_int64 (6 ),
275+ )
276+ ref = np .zeros ((10 , 5 ), dtype = np .float32 )
277+ ref [:6 , :] = 6.0
278+
279+ csdfg = copy_prog .compile ()
280+ csdfg .fast_call ((ctypes .byref (inp_struct ), ctypes .c_int (5 )), (ctypes .c_int (5 ),))
281+
282+ assert np .allclose (data , ref )
283+
284+
285+ def test_write_structure_in_loop ():
286+ M = dace .symbol ('M' )
287+ N = dace .symbol ('N' )
288+ Bundle = dace .data .Structure (members = {
289+ "data" : dace .data .Array (dace .float32 , (M , N )),
290+ "size" : dace .data .Scalar (dace .int64 )
291+ },
292+ name = "BundleType" )
293+
294+ @dace .program
295+ def init_prog (bundle : Bundle , fill_value : int ) -> None :
296+ for index in range (bundle .size ):
297+ bundle .data [index , :] = fill_value
298+
299+ data = np .zeros ((10 , 5 ), dtype = np .float32 )
300+ fill_value = 42
301+ inp_struct = Bundle .dtype .base_type .as_ctypes ()(
302+ data = data .__array_interface__ ['data' ][0 ],
303+ size = 6 ,
304+ )
305+ ref = np .zeros ((10 , 5 ), dtype = np .float32 )
306+ ref [:6 , :] = fill_value
307+
308+ init_prog .compile ()(inp_struct , fill_value , M = 10 , N = 5 )
309+
310+ assert np .allclose (data , ref )
311+
312+
227313if __name__ == '__main__' :
228314 test_read_structure ()
229315 test_write_structure ()
230316 test_local_structure ()
231317 test_rgf ()
232318 # test_read_structure_gpu()
319+ test_write_structure_in_map ()
320+ test_readwrite_structure_in_map ()
321+ test_write_structure_in_loop ()
0 commit comments