/[escript]/trunk/cusplibrary/cusp/gallery/stencil.inl
ViewVC logotype

Contents of /trunk/cusplibrary/cusp/gallery/stencil.inl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5148 - (show annotations)
Mon Sep 15 01:25:23 2014 UTC (5 years, 10 months ago) by caltinay
File size: 6892 byte(s)
Merging ripley diagonal storage + CUDA support into trunk.
Options file version has been incremented due to new options
'cuda' and 'nvccflags'.

1 /*
2 * Copyright 2008-2009 NVIDIA Corporation
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16
17 #include <cusp/copy.h>
18 #include <cusp/dia_matrix.h>
19
20 #include <thrust/tuple.h>
21 #include <thrust/reduce.h>
22 #include <thrust/scan.h>
23 #include <thrust/count.h>
24 #include <thrust/functional.h>
25
26 namespace cusp
27 {
28 namespace gallery
29 {
30 namespace detail
31 {
32
33 template <typename StencilPoint, typename GridDimension, typename IndexType, int i, int n>
34 struct inside_grid_helper
35 {
36 __host__ __device__
37 static bool inside_grid(StencilPoint point, GridDimension grid, IndexType index)
38 {
39 IndexType x = index % thrust::get<i>(grid) + thrust::get<i>(thrust::get<0>(point));
40
41 if (x < 0 || x >= thrust::get<i>(grid))
42 return false;
43 else
44 return inside_grid_helper<StencilPoint,GridDimension,IndexType,i + 1,n>::inside_grid(point, grid, index / thrust::get<i>(grid));
45 }
46 };
47
48 template <typename StencilPoint, typename GridDimension, typename IndexType, int n>
49 struct inside_grid_helper<StencilPoint,GridDimension,IndexType,n,n>
50 {
51 __host__ __device__
52 static bool inside_grid(StencilPoint point, GridDimension grid, IndexType index)
53 {
54 return true;
55 }
56 };
57
58 template <typename StencilPoint, typename GridDimension, typename IndexType>
59 __host__ __device__
60 bool inside_grid(StencilPoint point, GridDimension grid, IndexType index)
61 {
62 return inside_grid_helper<StencilPoint,GridDimension,IndexType, 0, thrust::tuple_size<GridDimension>::value >::inside_grid(point, grid, index);
63 }
64
65
66 template <typename Tuple, typename UnaryFunction, int i, int size>
67 struct tuple_for_each_helper
68 {
69 static UnaryFunction for_each(Tuple& t, UnaryFunction f)
70 {
71 f(thrust::get<i>(t));
72
73 return tuple_for_each_helper<Tuple,UnaryFunction,i + 1,size>::for_each(t, f);
74 }
75 };
76
77 template <typename Tuple, typename UnaryFunction, int size>
78 struct tuple_for_each_helper<Tuple,UnaryFunction,size,size>
79 {
80 static UnaryFunction for_each(Tuple& t, UnaryFunction f)
81 {
82 return f;
83 }
84 };
85
86 template <typename Tuple, typename UnaryFunction>
87 UnaryFunction tuple_for_each(Tuple& t, UnaryFunction f)
88 {
89 return tuple_for_each_helper<Tuple, UnaryFunction, 0, thrust::tuple_size<Tuple>::value>::for_each(t, f);
90 }
91
92
93 template <typename Tuple, typename OutputIterator>
94 struct unpack_tuple_functor
95 {
96 OutputIterator output;
97 unpack_tuple_functor(OutputIterator output) : output(output) {}
98
99 template <typename T>
100 void operator()(T v)
101 {
102 *output++ = v;
103 }
104 };
105
106
107 template <typename Tuple, typename OutputIterator>
108 OutputIterator unpack_tuple(const Tuple & t, OutputIterator output)
109 {
110 return tuple_for_each(t, unpack_tuple_functor<Tuple,OutputIterator>(output)).output;
111 }
112
113
114 template <typename IndexType,
115 typename ValueType,
116 typename StencilPoint,
117 typename GridDimension>
118 struct fill_diagonal_entries
119 {
120 StencilPoint point;
121 GridDimension grid;
122
123 fill_diagonal_entries(StencilPoint point, GridDimension grid)
124 : point(point), grid(grid) {}
125
126 __host__ __device__
127 ValueType operator()(IndexType index)
128 {
129 if (inside_grid(point, grid, index))
130 return thrust::get<1>(point);
131 else
132 return ValueType(0);
133 }
134 };
135
136 } // end namespace detail
137
138 template <typename IndexType,
139 typename ValueType,
140 typename MemorySpace,
141 typename StencilPoint,
142 typename GridDimension>
143 void generate_matrix_from_stencil( cusp::dia_matrix<IndexType,ValueType,MemorySpace>& matrix,
144 const cusp::array1d<StencilPoint,cusp::host_memory>& stencil,
145 const GridDimension& grid)
146 {
147 IndexType num_dimensions = thrust::tuple_size<GridDimension>::value;
148
149 cusp::array1d<IndexType,cusp::host_memory> grid_indices(num_dimensions);
150 detail::unpack_tuple(grid, grid_indices.begin());
151
152 IndexType num_rows = thrust::reduce(grid_indices.begin(), grid_indices.end(), IndexType(1), thrust::multiplies<IndexType>());
153
154 IndexType num_diagonals = stencil.size();
155
156 cusp::array1d<IndexType,cusp::host_memory> strides(grid_indices.size());
157 thrust::exclusive_scan(grid_indices.begin(), grid_indices.end(), strides.begin(), IndexType(1), thrust::multiplies<IndexType>());
158
159 cusp::array1d<IndexType,cusp::host_memory> offsets(stencil.size(), 0);
160 for(size_t i = 0; i < offsets.size(); i++)
161 {
162 cusp::array1d<IndexType,cusp::host_memory> stencil_indices(num_dimensions);
163 detail::unpack_tuple(thrust::get<0>(stencil[i]), stencil_indices.begin());
164
165 for(IndexType j = 0; j < num_dimensions; j++)
166 {
167 offsets[i] += strides[j] * stencil_indices[j];
168 }
169 }
170
171 // TODO compute num_entries directly from stencil
172 matrix.resize(num_rows, num_rows, 0, num_diagonals); // XXX we set NNZ to zero for now
173
174 cusp::copy(offsets, matrix.diagonal_offsets);
175
176 // ideally we'd have row views and column views here
177 for(IndexType i = 0; i < num_diagonals; i++)
178 {
179 thrust::transform(thrust::counting_iterator<IndexType>(0),
180 thrust::counting_iterator<IndexType>(num_rows),
181 matrix.values.values.begin() + matrix.values.pitch * i,
182 detail::fill_diagonal_entries<IndexType,ValueType,StencilPoint,GridDimension>(stencil[i], grid));
183 }
184
185 matrix.num_entries = matrix.values.values.size() - thrust::count(matrix.values.values.begin(), matrix.values.values.end(), ValueType(0));
186 }
187
188 // TODO add an entry point and make this the default path
189 template <typename MatrixType,
190 typename StencilPoint,
191 typename GridDimension>
192 void generate_matrix_from_stencil( MatrixType& matrix,
193 const cusp::array1d<StencilPoint,cusp::host_memory>& stencil,
194 const GridDimension& grid)
195 {
196 CUSP_PROFILE_SCOPED();
197
198 typedef typename MatrixType::index_type IndexType;
199 typedef typename MatrixType::value_type ValueType;
200 typedef typename MatrixType::memory_space MemorySpace;
201
202 cusp::dia_matrix<IndexType,ValueType,MemorySpace> dia;
203 generate_matrix_from_stencil(dia, stencil, grid);
204
205 cusp::convert(dia, matrix);
206 }
207
208 } // end namespace gallery
209 } // end namespace cusp
210

  ViewVC Help
Powered by ViewVC 1.1.26