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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5209 - (show annotations)
Tue Oct 21 06:50:20 2014 UTC (5 years, 7 months ago) by caltinay
File MIME type: text/plain
File size: 13440 byte(s)
Adding memory optimization for symmetric CDS matrices (in ripley).
Only main and upper diagonal blocks are now stored.
SpMV not updated on GPUs yet but we get a speed-up on Xeons.


1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2014 by University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17 /*! \file cds_matrix.h
18 * \brief Column Diagonal Block 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 cds_matrix_view;
46
47 /*! \p cds_matrix : Diagonal block matrix container (always square)
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 8-by-8
56 * \p cds_matrix on the host with 2 diagonal 2x2 blocks (16 total nonzeros)
57 * and then copies the matrix to the device.
58 *
59 * \code
60 * #include <cusp/cds_matrix.h>
61 * ...
62 *
63 * // allocate storage for (8,8) matrix with 16 nonzeros in 2 diagonal blocks
64 * // of size 2x2
65 * cusp::cds_matrix<int,float,cusp::host_memory> A(8,16,2,2);
66 *
67 * // initialize diagonal offsets
68 * A.diagonal_offsets[0] = -2;
69 * A.diagonal_offsets[1] = 2;
70 *
71 * // initialize diagonal values
72 *
73 * // first diagonal block
74 * A.values( 0,0) = 0; // outside matrix
75 * A.values( 1,0) = 0; // outside matrix
76 * A.values( 2,0) = 0; // outside matrix
77 * A.values( 3,0) = 0; // outside matrix
78 * A.values( 4,0) = 0; // outside matrix
79 * A.values( 5,0) = 0; // outside matrix
80 * A.values( 6,0) = 0; // outside matrix
81 * A.values( 7,0) = 0; // outside matrix
82 * A.values( 8,0) = -10;
83 * A.values( 9,0) = -20;
84 * A.values(10,0) = -30;
85 * A.values(11,0) = -40;
86 * A.values(12,0) = -50;
87 * A.values(13,0) = -60;
88 * A.values(14,0) = -70;
89 * A.values(15,0) = -80;
90 *
91 * // second diagonal block
92 * A.values( 0,1) = 10;
93 * A.values( 1,1) = 20;
94 * A.values( 2,1) = 30;
95 * A.values( 3,1) = 40;
96 * A.values( 4,1) = 50;
97 * A.values( 5,1) = 60;
98 * A.values( 6,1) = 70;
99 * A.values( 7,1) = 80;
100 * A.values( 8,1) = 0; // outside matrix
101 * A.values( 9,1) = 0; // outside matrix
102 * A.values(10,1) = 0; // outside matrix
103 * A.values(11,1) = 0; // outside matrix
104 * A.values(12,1) = 0; // outside matrix
105 * A.values(13,1) = 0; // outside matrix
106 * A.values(14,1) = 0; // outside matrix
107 * A.values(15,1) = 0; // outside matrix
108 *
109 * // A now represents the following matrix
110 * // [ 0 0 0 0 10 30 0 0]
111 * // [ 0 0 0 0 20 40 0 0]
112 * // [ 0 0 0 0 0 0 50 70]
113 * // [ 0 0 0 0 0 0 60 80]
114 * // [-10 -30 0 0 0 0 0 0]
115 * // [-20 -40 0 0 0 0 0 0]
116 * // [ 0 0 -50 -70 0 0 0 0]
117 * // [ 0 0 -60 -80 0 0 0 0]
118 *
119 * // copy to the device
120 * cusp::cds_matrix<int,float,cusp::device_memory> B = A;
121 * \endcode
122 *
123 */
124 template <typename IndexType, typename ValueType, class MemorySpace>
125 class cds_matrix : public detail::matrix_base<IndexType,ValueType,MemorySpace,cusp::cds_format>
126 {
127 typedef cusp::detail::matrix_base<IndexType,ValueType,MemorySpace,cusp::cds_format> Parent;
128 public:
129 // TODO statically assert is_signed<IndexType>
130
131 /*! rebind matrix to a different MemorySpace
132 */
133 template<typename MemorySpace2>
134 struct rebind { typedef cusp::cds_matrix<IndexType, ValueType, MemorySpace2> type; };
135
136 /*! type of diagonal offsets array
137 */
138 typedef typename cusp::array1d<IndexType, MemorySpace> diagonal_offsets_array_type;
139
140 /*! type of values array
141 */
142 typedef typename cusp::array2d<ValueType, MemorySpace, cusp::column_major> values_array_type;
143
144 /*! equivalent container type
145 */
146 typedef typename cusp::cds_matrix<IndexType, ValueType, MemorySpace> container;
147
148 /*! equivalent view type
149 */
150 typedef typename cusp::cds_matrix_view<typename diagonal_offsets_array_type::view,
151 typename values_array_type::view,
152 IndexType, ValueType, MemorySpace> view;
153
154 /*! equivalent const_view type
155 */
156 typedef typename cusp::cds_matrix_view<typename diagonal_offsets_array_type::const_view,
157 typename values_array_type::const_view,
158 IndexType, ValueType, MemorySpace> const_view;
159
160 /*! Storage for the diagonal block offsets.
161 */
162 diagonal_offsets_array_type diagonal_offsets;
163
164 /*! Storage for the nonzero entries of the CDS data structure.
165 */
166 values_array_type values;
167
168 /*! Diagonal block size.
169 */
170 size_t block_size;
171
172 /*! Indicator for matrix symmetry. If true, only non-negative diagonal
173 * offsets are stored.
174 */
175 bool symmetric;
176
177 /*! Construct an empty \p cds_matrix.
178 */
179 cds_matrix() : block_size(1), symmetric(false) {}
180
181 /*! Construct a \p cds_matrix with a specific shape, number of nonzero entries,
182 * and number of occupied diagonals.
183 *
184 * \param num_rows Number of rows & columns.
185 * \param num_entries Number of nonzero matrix entries.
186 * \param num_diagonals Number of occupied diagonals.
187 * \param alignment Amount of padding used to align the data structure (default 32).
188 */
189 cds_matrix(size_t num_rows, size_t num_entries,
190 size_t num_diagonals, size_t blocksize, bool is_symmetric,
191 size_t alignment = 32)
192 : Parent(num_rows, num_rows, num_entries),
193 diagonal_offsets(num_diagonals),
194 block_size(blocksize),
195 symmetric(is_symmetric)
196 {
197 // TODO use array2d constructor when it can accept pitch
198 values.resize(num_rows, num_diagonals*blocksize, detail::round_up(num_rows, alignment));
199 }
200
201 /*! Construct a \p cds_matrix from another matrix.
202 *
203 * \param matrix Another sparse or dense matrix.
204 */
205 template <typename MatrixType>
206 cds_matrix(const MatrixType& matrix);
207
208 /*! Resize matrix dimensions and underlying storage
209 */
210 void resize(size_t num_rows, size_t num_entries, size_t num_diagonals,
211 size_t blocksize)
212 {
213 Parent::resize(num_rows, num_rows, num_entries);
214 diagonal_offsets.resize(num_diagonals);
215 block_size = blocksize;
216 values.resize(num_rows, num_diagonals*blocksize);
217 }
218
219 /*! Resize matrix dimensions and underlying storage
220 */
221 void resize(size_t num_rows, size_t num_entries, size_t num_diagonals,
222 size_t blocksize, size_t alignment)
223 {
224 Parent::resize(num_rows, num_rows, num_entries);
225 diagonal_offsets.resize(num_diagonals);
226 block_size = blocksize;
227 values.resize(num_rows, num_diagonals*blocksize, detail::round_up(num_rows, alignment));
228 }
229
230 /*! Swap the contents of two \p cds_matrix objects.
231 *
232 * \param matrix Another \p cds_matrix with the same IndexType and ValueType.
233 */
234 void swap(cds_matrix& matrix)
235 {
236 Parent::swap(matrix);
237 diagonal_offsets.swap(matrix.diagonal_offsets);
238 thrust::swap(block_size, matrix.block_size);
239 thrust::swap(symmetric, matrix.symmetric);
240 values.swap(matrix.values);
241 }
242
243 /*! Assignment from another matrix.
244 *
245 * \param matrix Another sparse or dense matrix.
246 */
247 template <typename MatrixType>
248 cds_matrix& operator=(const MatrixType& matrix);
249 }; // class cds_matrix
250 /*! \}
251 */
252
253 /*! \addtogroup sparse_matrix_views Sparse Matrix Views
254 * \ingroup sparse_matrices
255 * \{
256 */
257
258 /*! \p cds_matrix_view : Diagonal block matrix view
259 *
260 * \tparam Array1 Type of \c diagonal_offsets
261 * \tparam Array2 Type of \c values array view
262 * \tparam IndexType Type used for matrix indices (e.g. \c int).
263 * \tparam ValueType Type used for matrix values (e.g. \c float).
264 * \tparam MemorySpace A memory space (e.g. \c cusp::host_memory or cusp::device_memory)
265 *
266 */
267 template <typename Array1,
268 typename Array2,
269 typename IndexType = typename Array1::value_type,
270 typename ValueType = typename Array2::value_type,
271 typename MemorySpace = typename cusp::minimum_space<typename Array1::memory_space, typename Array2::memory_space>::type >
272 class cds_matrix_view : public detail::matrix_base<IndexType,ValueType,MemorySpace,cusp::cds_format>
273 {
274 typedef cusp::detail::matrix_base<IndexType,ValueType,MemorySpace,cusp::cds_format> Parent;
275 public:
276 /*! type of \c diagonal_offsets array
277 */
278 typedef Array1 diagonal_offsets_array_type;
279
280 /*! type of \c column_indices array
281 */
282 typedef Array2 values_array_type;
283
284 /*! equivalent container type
285 */
286 typedef typename cusp::cds_matrix<IndexType, ValueType, MemorySpace> container;
287
288 /*! equivalent view type
289 */
290 typedef typename cusp::cds_matrix_view<Array1, Array2, IndexType, ValueType, MemorySpace> view;
291
292 /*! Storage for the diagonal offsets.
293 */
294 diagonal_offsets_array_type diagonal_offsets;
295
296 /*! Storage for the nonzero entries of the DIA data structure.
297 */
298 values_array_type values;
299
300 /*! Size of the diagonal blocks.
301 */
302 size_t block_size;
303
304 /*! Construct an empty \p cds_matrix_view.
305 */
306 cds_matrix_view() {}
307
308 template <typename OtherArray1, typename OtherArray2>
309 cds_matrix_view(size_t num_rows, size_t num_entries, size_t blocksize,
310 OtherArray1& diagonal_offsets, OtherArray2& values)
311 : Parent(num_rows, num_rows, num_entries), block_size(blocksize),
312 diagonal_offsets(diagonal_offsets), values(values) {}
313
314 template <typename OtherArray1, typename OtherArray2>
315 cds_matrix_view(size_t num_rows, size_t num_entries, size_t blocksize,
316 const OtherArray1& diagonal_offsets, const OtherArray2& values)
317 : Parent(num_rows, num_rows, num_entries), block_size(blocksize),
318 diagonal_offsets(diagonal_offsets), values(values) {}
319
320 template <typename Matrix>
321 cds_matrix_view(Matrix& A)
322 : Parent(A), block_size(A.block_size), diagonal_offsets(A.diagonal_offsets), values(A.values) {}
323
324 template <typename Matrix>
325 cds_matrix_view(const Matrix& A)
326 : Parent(A), block_size(A.block_size), diagonal_offsets(A.diagonal_offsets), values(A.values) {}
327
328 /*! Resize matrix dimensions and underlying storage
329 */
330 void resize(size_t num_rows, size_t num_entries, size_t num_diagonals,
331 size_t blocksize)
332 {
333 Parent::resize(num_rows, num_rows, num_entries);
334 diagonal_offsets.resize(num_diagonals);
335 block_size = blocksize;
336 values.resize(num_rows, num_diagonals*blocksize);
337 }
338
339 /*! Resize matrix dimensions and underlying storage
340 */
341 void resize(size_t num_rows, size_t num_entries, size_t num_diagonals,
342 size_t blocksize, size_t alignment)
343 {
344 Parent::resize(num_rows, num_rows, num_entries);
345 diagonal_offsets.resize(num_diagonals);
346 block_size = blocksize;
347 values.resize(num_rows, num_diagonals*blocksize, detail::round_up(num_rows, alignment));
348 }
349 }; // class cds_matrix_view
350
351
352 template <typename Array1,
353 typename Array2>
354 cds_matrix_view<Array1,Array2>
355 make_cds_matrix_view(size_t num_rows,
356 size_t num_entries,
357 size_t block_size,
358 Array1 diagonal_offsets,
359 Array2 values);
360
361 template <typename Array1,
362 typename Array2,
363 typename IndexType,
364 typename ValueType,
365 typename MemorySpace>
366 cds_matrix_view<Array1,Array2,IndexType,ValueType,MemorySpace>
367 make_cds_matrix_view(const cds_matrix_view<Array1,Array2,IndexType,ValueType,MemorySpace>& m);
368
369 template <typename IndexType, typename ValueType, class MemorySpace>
370 typename cds_matrix<IndexType,ValueType,MemorySpace>::view
371 make_cds_matrix_view(cds_matrix<IndexType,ValueType,MemorySpace>& m);
372
373 template <typename IndexType, typename ValueType, class MemorySpace>
374 typename cds_matrix<IndexType,ValueType,MemorySpace>::const_view
375 make_cds_matrix_view(const cds_matrix<IndexType,ValueType,MemorySpace>& m);
376 /*! \} // end Views
377 */
378
379 } // end namespace cusp
380
381 #include <cusp/array2d.h>
382 #include <cusp/detail/cds_matrix.inl>
383

  ViewVC Help
Powered by ViewVC 1.1.26