/[escript]/trunk/cusplibrary/cusp/dia_matrix.h
ViewVC logotype

Contents of /trunk/cusplibrary/cusp/dia_matrix.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5148 - (show annotations)
Mon Sep 15 01:25:23 2014 UTC (5 years, 9 months ago) by caltinay
File MIME type: text/plain
File size: 11658 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 /*! \file dia_matrix.h
18 * \brief Diagonal matrix format.
19 */
20
21 #pragma once
22
23 #include <cusp/detail/config.h>
24
25 #include <cusp/array1d.h>
26 #include <cusp/format.h>
27 #include <cusp/detail/matrix_base.h>
28 #include <cusp/detail/utils.h>
29
30 namespace cusp
31 {
32
33 /*! \addtogroup sparse_matrices Sparse Matrices
34 */
35
36 /*! \addtogroup sparse_matrix_containers Sparse Matrix Containers
37 * \ingroup sparse_matrices
38 * \{
39 */
40
41 // Forward definitions
42 struct column_major;
43 template<typename ValueType, class MemorySpace, class Orientation> class array2d;
44 template<typename Array, class Orientation> class array2d_view;
45 template <typename Array1, typename Array2, typename IndexType, typename ValueType, typename MemorySpace> class dia_matrix_view;
46
47 /*! \p dia_matrix : Diagonal matrix container
48 *
49 * \tparam IndexType Type used for matrix indices (e.g. \c int).
50 * \tparam ValueType Type used for matrix values (e.g. \c float).
51 * \tparam MemorySpace A memory space (e.g. \c cusp::host_memory or cusp::device_memory)
52 *
53 * \note The diagonal offsets should not contain duplicate entries.
54 *
55 * The following code snippet demonstrates how to create a 4-by-3
56 * \p dia_matrix on the host with 3 diagonals (6 total nonzeros)
57 * and then copies the matrix to the device.
58 *
59 * \code
60 * #include <cusp/dia_matrix.h>
61 * ...
62 *
63 * // allocate storage for (4,3) matrix with 6 nonzeros in 3 diagonals
64 * cusp::dia_matrix<int,float,cusp::host_memory> A(4,3,6,3);
65 *
66 * // initialize diagonal offsets
67 * A.diagonal_offsets[0] = -2;
68 * A.diagonal_offsets[1] = 0;
69 * A.diagonal_offsets[2] = 1;
70 *
71 * // initialize diagonal values
72 *
73 * // first diagonal
74 * A.values(0,2) = 0; // outside matrix
75 * A.values(1,2) = 0; // outside matrix
76 * A.values(2,0) = 40;
77 * A.values(3,0) = 60;
78 *
79 * // second diagonal
80 * A.values(0,1) = 10;
81 * A.values(1,1) = 0;
82 * A.values(2,1) = 50;
83 * A.values(3,1) = 50; // outside matrix
84 *
85 * // third diagonal
86 * A.values(0,2) = 20;
87 * A.values(1,2) = 30;
88 * A.values(2,2) = 0; // outside matrix
89 * A.values(3,2) = 0; // outside matrix
90 *
91 * // A now represents the following matrix
92 * // [10 20 0]
93 * // [ 0 0 30]
94 * // [40 0 50]
95 * // [ 0 60 0]
96 *
97 * // copy to the device
98 * cusp::dia_matrix<int,float,cusp::device_memory> B = A;
99 * \endcode
100 *
101 */
102 template <typename IndexType, typename ValueType, class MemorySpace>
103 class dia_matrix : public detail::matrix_base<IndexType,ValueType,MemorySpace,cusp::dia_format>
104 {
105 typedef cusp::detail::matrix_base<IndexType,ValueType,MemorySpace,cusp::dia_format> Parent;
106 public:
107 // TODO statically assert is_signed<IndexType>
108
109 /*! rebind matrix to a different MemorySpace
110 */
111 template<typename MemorySpace2>
112 struct rebind { typedef cusp::dia_matrix<IndexType, ValueType, MemorySpace2> type; };
113
114 /*! type of diagonal offsets array
115 */
116 typedef typename cusp::array1d<IndexType, MemorySpace> diagonal_offsets_array_type;
117
118 /*! type of values array
119 */
120 typedef typename cusp::array2d<ValueType, MemorySpace, cusp::column_major> values_array_type;
121
122 /*! equivalent container type
123 */
124 typedef typename cusp::dia_matrix<IndexType, ValueType, MemorySpace> container;
125
126 /*! equivalent view type
127 */
128 typedef typename cusp::dia_matrix_view<typename diagonal_offsets_array_type::view,
129 typename values_array_type::view,
130 IndexType, ValueType, MemorySpace> view;
131
132 /*! equivalent const_view type
133 */
134 typedef typename cusp::dia_matrix_view<typename diagonal_offsets_array_type::const_view,
135 typename values_array_type::const_view,
136 IndexType, ValueType, MemorySpace> const_view;
137
138 /*! Storage for the diagonal offsets.
139 */
140 diagonal_offsets_array_type diagonal_offsets;
141
142 /*! Storage for the nonzero entries of the DIA data structure.
143 */
144 values_array_type values;
145
146 /*! Construct an empty \p dia_matrix.
147 */
148 dia_matrix() {}
149
150 /*! Construct a \p dia_matrix with a specific shape, number of nonzero entries,
151 * and number of occupied diagonals.
152 *
153 * \param num_rows Number of rows.
154 * \param num_cols Number of columns.
155 * \param num_entries Number of nonzero matrix entries.
156 * \param num_diagonals Number of occupied diagonals.
157 * \param alignment Amount of padding used to align the data structure (default 32).
158 */
159 dia_matrix(size_t num_rows, size_t num_cols, size_t num_entries,
160 size_t num_diagonals, size_t alignment = 32)
161 : Parent(num_rows, num_cols, num_entries),
162 diagonal_offsets(num_diagonals)
163 {
164 // TODO use array2d constructor when it can accept pitch
165 values.resize(num_rows, num_diagonals, detail::round_up(num_rows, alignment));
166 }
167
168 /*! Construct a \p dia_matrix from another matrix.
169 *
170 * \param matrix Another sparse or dense matrix.
171 */
172 template <typename MatrixType>
173 dia_matrix(const MatrixType& matrix);
174
175 /*! Resize matrix dimensions and underlying storage
176 */
177 void resize(size_t num_rows, size_t num_cols, size_t num_entries,
178 size_t num_diagonals)
179 {
180 Parent::resize(num_rows, num_cols, num_entries);
181 diagonal_offsets.resize(num_diagonals);
182 values.resize(num_rows, num_diagonals);
183 }
184
185 /*! Resize matrix dimensions and underlying storage
186 */
187 void resize(size_t num_rows, size_t num_cols, size_t num_entries,
188 size_t num_diagonals, size_t alignment)
189 {
190 Parent::resize(num_rows, num_cols, num_entries);
191 diagonal_offsets.resize(num_diagonals);
192 values.resize(num_rows, num_diagonals, detail::round_up(num_rows, alignment));
193 }
194
195 /*! Swap the contents of two \p dia_matrix objects.
196 *
197 * \param matrix Another \p dia_matrix with the same IndexType and ValueType.
198 */
199 void swap(dia_matrix& matrix)
200 {
201 Parent::swap(matrix);
202 diagonal_offsets.swap(matrix.diagonal_offsets);
203 values.swap(matrix.values);
204 }
205
206 /*! Assignment from another matrix.
207 *
208 * \param matrix Another sparse or dense matrix.
209 */
210 template <typename MatrixType>
211 dia_matrix& operator=(const MatrixType& matrix);
212 }; // class dia_matrix
213 /*! \}
214 */
215
216 /*! \addtogroup sparse_matrix_views Sparse Matrix Views
217 * \ingroup sparse_matrices
218 * \{
219 */
220
221 /*! \p dia_matrix_view : Diagonal matrix view
222 *
223 * \tparam Array1 Type of \c diagonal_offsets
224 * \tparam Array2 Type of \c values array view
225 * \tparam IndexType Type used for matrix indices (e.g. \c int).
226 * \tparam ValueType Type used for matrix values (e.g. \c float).
227 * \tparam MemorySpace A memory space (e.g. \c cusp::host_memory or cusp::device_memory)
228 *
229 */
230 template <typename Array1,
231 typename Array2,
232 typename IndexType = typename Array1::value_type,
233 typename ValueType = typename Array2::value_type,
234 typename MemorySpace = typename cusp::minimum_space<typename Array1::memory_space, typename Array2::memory_space>::type >
235 class dia_matrix_view : public detail::matrix_base<IndexType,ValueType,MemorySpace,cusp::dia_format>
236 {
237 typedef cusp::detail::matrix_base<IndexType,ValueType,MemorySpace,cusp::dia_format> Parent;
238 public:
239 /*! type of \c diagonal_offsets array
240 */
241 typedef Array1 diagonal_offsets_array_type;
242
243 /*! type of \c column_indices array
244 */
245 typedef Array2 values_array_type;
246
247 /*! equivalent container type
248 */
249 typedef typename cusp::dia_matrix<IndexType, ValueType, MemorySpace> container;
250
251 /*! equivalent view type
252 */
253 typedef typename cusp::dia_matrix_view<Array1, Array2, IndexType, ValueType, MemorySpace> view;
254
255 /*! Storage for the diagonal offsets.
256 */
257 diagonal_offsets_array_type diagonal_offsets;
258
259 /*! Storage for the nonzero entries of the DIA data structure.
260 */
261 values_array_type values;
262
263 /*! Construct an empty \p dia_matrix_view.
264 */
265 dia_matrix_view() {}
266
267 template <typename OtherArray1, typename OtherArray2>
268 dia_matrix_view(size_t num_rows, size_t num_cols, size_t num_entries,
269 OtherArray1& diagonal_offsets, OtherArray2& values)
270 : Parent(num_rows, num_cols, num_entries), diagonal_offsets(diagonal_offsets), values(values) {}
271
272 template <typename OtherArray1, typename OtherArray2>
273 dia_matrix_view(size_t num_rows, size_t num_cols, size_t num_entries,
274 const OtherArray1& diagonal_offsets, const OtherArray2& values)
275 : Parent(num_rows, num_cols, num_entries), diagonal_offsets(diagonal_offsets), values(values) {}
276
277 template <typename Matrix>
278 dia_matrix_view(Matrix& A)
279 : Parent(A), diagonal_offsets(A.diagonal_offsets), values(A.values) {}
280
281 template <typename Matrix>
282 dia_matrix_view(const Matrix& A)
283 : Parent(A), diagonal_offsets(A.diagonal_offsets), values(A.values) {}
284
285 /*! Resize matrix dimensions and underlying storage
286 */
287 void resize(size_t num_rows, size_t num_cols, size_t num_entries,
288 size_t num_diagonals)
289 {
290 Parent::resize(num_rows, num_cols, num_entries);
291 diagonal_offsets.resize(num_diagonals);
292 values.resize(num_rows, num_diagonals);
293 }
294
295 /*! Resize matrix dimensions and underlying storage
296 */
297 void resize(size_t num_rows, size_t num_cols, size_t num_entries,
298 size_t num_diagonals, size_t alignment)
299 {
300 Parent::resize(num_rows, num_cols, num_entries);
301 diagonal_offsets.resize(num_diagonals);
302 values.resize(num_rows, num_diagonals, detail::round_up(num_rows, alignment));
303 }
304 }; // class dia_matrix_view
305
306
307 template <typename Array1,
308 typename Array2>
309 dia_matrix_view<Array1,Array2>
310 make_dia_matrix_view(size_t num_rows,
311 size_t num_cols,
312 size_t num_entries,
313 Array1 diagonal_offsets,
314 Array2 values);
315
316 template <typename Array1,
317 typename Array2,
318 typename IndexType,
319 typename ValueType,
320 typename MemorySpace>
321 dia_matrix_view<Array1,Array2,IndexType,ValueType,MemorySpace>
322 make_dia_matrix_view(const dia_matrix_view<Array1,Array2,IndexType,ValueType,MemorySpace>& m);
323
324 template <typename IndexType, typename ValueType, class MemorySpace>
325 typename dia_matrix<IndexType,ValueType,MemorySpace>::view
326 make_dia_matrix_view(dia_matrix<IndexType,ValueType,MemorySpace>& m);
327
328 template <typename IndexType, typename ValueType, class MemorySpace>
329 typename dia_matrix<IndexType,ValueType,MemorySpace>::const_view
330 make_dia_matrix_view(const dia_matrix<IndexType,ValueType,MemorySpace>& m);
331 /*! \} // end Views
332 */
333
334 } // end namespace cusp
335
336 #include <cusp/array2d.h>
337 #include <cusp/detail/dia_matrix.inl>
338

  ViewVC Help
Powered by ViewVC 1.1.26