/[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 5588 - (show annotations)
Tue Apr 21 05:07:16 2015 UTC (5 years, 2 months ago) by caltinay
File MIME type: text/plain
File size: 13494 byte(s)
Update to copyright headers in cusp library.

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

  ViewVC Help
Powered by ViewVC 1.1.26