/[escript]/branches/diaplayground/cusplibrary/cusp/cds_matrix.h
ViewVC logotype

Contents of /branches/diaplayground/cusplibrary/cusp/cds_matrix.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5095 - (show annotations)
Fri Jul 11 05:54:25 2014 UTC (5 years, 10 months ago) by caltinay
File MIME type: text/plain
File size: 13177 byte(s)
Added implementation of column diagonal storage (CDS) matrix format
and relevant operations within cusp to allow solving for block sizes > 1.
Added exception if trying to solve with more than one rank.
Minor clean up.

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 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 /*! Construct an empty \p cds_matrix.
173 */
174 cds_matrix() {}
175
176 /*! Construct a \p cds_matrix with a specific shape, number of nonzero entries,
177 * and number of occupied diagonals.
178 *
179 * \param num_rows Number of rows & columns.
180 * \param num_entries Number of nonzero matrix entries.
181 * \param num_diagonals Number of occupied diagonals.
182 * \param alignment Amount of padding used to align the data structure (default 32).
183 */
184 cds_matrix(size_t num_rows, size_t num_entries,
185 size_t num_diagonals, size_t blocksize, size_t alignment = 32)
186 : Parent(num_rows, num_rows, num_entries),
187 diagonal_offsets(num_diagonals),
188 block_size(blocksize)
189 {
190 // TODO use array2d constructor when it can accept pitch
191 values.resize(num_rows, num_diagonals*blocksize, detail::round_up(num_rows, alignment));
192 }
193
194 /*! Construct a \p cds_matrix from another matrix.
195 *
196 * \param matrix Another sparse or dense matrix.
197 */
198 template <typename MatrixType>
199 cds_matrix(const MatrixType& matrix);
200
201 /*! Resize matrix dimensions and underlying storage
202 */
203 void resize(size_t num_rows, size_t num_entries, size_t num_diagonals,
204 size_t blocksize)
205 {
206 Parent::resize(num_rows, num_rows, num_entries);
207 diagonal_offsets.resize(num_diagonals);
208 block_size = blocksize;
209 values.resize(num_rows, num_diagonals*blocksize);
210 }
211
212 /*! Resize matrix dimensions and underlying storage
213 */
214 void resize(size_t num_rows, size_t num_entries, size_t num_diagonals,
215 size_t blocksize, size_t alignment)
216 {
217 Parent::resize(num_rows, num_rows, num_entries);
218 diagonal_offsets.resize(num_diagonals);
219 block_size = blocksize;
220 values.resize(num_rows, num_diagonals*blocksize, detail::round_up(num_rows, alignment));
221 }
222
223 /*! Swap the contents of two \p cds_matrix objects.
224 *
225 * \param matrix Another \p cds_matrix with the same IndexType and ValueType.
226 */
227 void swap(cds_matrix& matrix)
228 {
229 Parent::swap(matrix);
230 diagonal_offsets.swap(matrix.diagonal_offsets);
231 thrust::swap(block_size, matrix.block_size);
232 values.swap(matrix.values);
233 }
234
235 /*! Assignment from another matrix.
236 *
237 * \param matrix Another sparse or dense matrix.
238 */
239 template <typename MatrixType>
240 cds_matrix& operator=(const MatrixType& matrix);
241 }; // class cds_matrix
242 /*! \}
243 */
244
245 /*! \addtogroup sparse_matrix_views Sparse Matrix Views
246 * \ingroup sparse_matrices
247 * \{
248 */
249
250 /*! \p cds_matrix_view : Diagonal block matrix view
251 *
252 * \tparam Array1 Type of \c diagonal_offsets
253 * \tparam Array2 Type of \c values array view
254 * \tparam IndexType Type used for matrix indices (e.g. \c int).
255 * \tparam ValueType Type used for matrix values (e.g. \c float).
256 * \tparam MemorySpace A memory space (e.g. \c cusp::host_memory or cusp::device_memory)
257 *
258 */
259 template <typename Array1,
260 typename Array2,
261 typename IndexType = typename Array1::value_type,
262 typename ValueType = typename Array2::value_type,
263 typename MemorySpace = typename cusp::minimum_space<typename Array1::memory_space, typename Array2::memory_space>::type >
264 class cds_matrix_view : public detail::matrix_base<IndexType,ValueType,MemorySpace,cusp::cds_format>
265 {
266 typedef cusp::detail::matrix_base<IndexType,ValueType,MemorySpace,cusp::cds_format> Parent;
267 public:
268 /*! type of \c diagonal_offsets array
269 */
270 typedef Array1 diagonal_offsets_array_type;
271
272 /*! type of \c column_indices array
273 */
274 typedef Array2 values_array_type;
275
276 /*! equivalent container type
277 */
278 typedef typename cusp::cds_matrix<IndexType, ValueType, MemorySpace> container;
279
280 /*! equivalent view type
281 */
282 typedef typename cusp::cds_matrix_view<Array1, Array2, IndexType, ValueType, MemorySpace> view;
283
284 /*! Storage for the diagonal offsets.
285 */
286 diagonal_offsets_array_type diagonal_offsets;
287
288 /*! Storage for the nonzero entries of the DIA data structure.
289 */
290 values_array_type values;
291
292 /*! Size of the diagonal blocks.
293 */
294 size_t block_size;
295
296 /*! Construct an empty \p cds_matrix_view.
297 */
298 cds_matrix_view() {}
299
300 template <typename OtherArray1, typename OtherArray2>
301 cds_matrix_view(size_t num_rows, size_t num_entries, size_t blocksize,
302 OtherArray1& diagonal_offsets, OtherArray2& values)
303 : Parent(num_rows, num_rows, num_entries), block_size(blocksize),
304 diagonal_offsets(diagonal_offsets), values(values) {}
305
306 template <typename OtherArray1, typename OtherArray2>
307 cds_matrix_view(size_t num_rows, size_t num_entries, size_t blocksize,
308 const OtherArray1& diagonal_offsets, const OtherArray2& values)
309 : Parent(num_rows, num_rows, num_entries), block_size(blocksize),
310 diagonal_offsets(diagonal_offsets), values(values) {}
311
312 template <typename Matrix>
313 cds_matrix_view(Matrix& A)
314 : Parent(A), block_size(A.block_size), diagonal_offsets(A.diagonal_offsets), values(A.values) {}
315
316 template <typename Matrix>
317 cds_matrix_view(const Matrix& A)
318 : Parent(A), block_size(A.block_size), diagonal_offsets(A.diagonal_offsets), values(A.values) {}
319
320 /*! Resize matrix dimensions and underlying storage
321 */
322 void resize(size_t num_rows, size_t num_entries, size_t num_diagonals,
323 size_t blocksize)
324 {
325 Parent::resize(num_rows, num_rows, num_entries);
326 diagonal_offsets.resize(num_diagonals);
327 block_size = blocksize;
328 values.resize(num_rows, num_diagonals*blocksize);
329 }
330
331 /*! Resize matrix dimensions and underlying storage
332 */
333 void resize(size_t num_rows, size_t num_entries, size_t num_diagonals,
334 size_t blocksize, size_t alignment)
335 {
336 Parent::resize(num_rows, num_rows, num_entries);
337 diagonal_offsets.resize(num_diagonals);
338 block_size = blocksize;
339 values.resize(num_rows, num_diagonals*blocksize, detail::round_up(num_rows, alignment));
340 }
341 }; // class cds_matrix_view
342
343
344 template <typename Array1,
345 typename Array2>
346 cds_matrix_view<Array1,Array2>
347 make_cds_matrix_view(size_t num_rows,
348 size_t num_entries,
349 size_t block_size,
350 Array1 diagonal_offsets,
351 Array2 values);
352
353 template <typename Array1,
354 typename Array2,
355 typename IndexType,
356 typename ValueType,
357 typename MemorySpace>
358 cds_matrix_view<Array1,Array2,IndexType,ValueType,MemorySpace>
359 make_cds_matrix_view(const cds_matrix_view<Array1,Array2,IndexType,ValueType,MemorySpace>& m);
360
361 template <typename IndexType, typename ValueType, class MemorySpace>
362 typename cds_matrix<IndexType,ValueType,MemorySpace>::view
363 make_cds_matrix_view(cds_matrix<IndexType,ValueType,MemorySpace>& m);
364
365 template <typename IndexType, typename ValueType, class MemorySpace>
366 typename cds_matrix<IndexType,ValueType,MemorySpace>::const_view
367 make_cds_matrix_view(const cds_matrix<IndexType,ValueType,MemorySpace>& m);
368 /*! \} // end Views
369 */
370
371 } // end namespace cusp
372
373 #include <cusp/array2d.h>
374 #include <cusp/detail/cds_matrix.inl>
375

  ViewVC Help
Powered by ViewVC 1.1.26