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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5209 - (hide annotations)
Tue Oct 21 06:50:20 2014 UTC (5 years, 8 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 caltinay 5209
2     /*****************************************************************************
3 caltinay 5095 *
4 caltinay 5209 * Copyright (c) 2014 by University of Queensland
5     * http://www.uq.edu.au
6 caltinay 5095 *
7 caltinay 5209 * 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 caltinay 5095 *
11 caltinay 5209 * 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 caltinay 5095
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 caltinay 5209 /*! Indicator for matrix symmetry. If true, only non-negative diagonal
173     * offsets are stored.
174     */
175     bool symmetric;
176    
177 caltinay 5095 /*! Construct an empty \p cds_matrix.
178     */
179 caltinay 5209 cds_matrix() : block_size(1), symmetric(false) {}
180 caltinay 5095
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 caltinay 5209 size_t num_diagonals, size_t blocksize, bool is_symmetric,
191     size_t alignment = 32)
192 caltinay 5095 : Parent(num_rows, num_rows, num_entries),
193     diagonal_offsets(num_diagonals),
194 caltinay 5209 block_size(blocksize),
195     symmetric(is_symmetric)
196 caltinay 5095 {
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 caltinay 5209 thrust::swap(symmetric, matrix.symmetric);
240 caltinay 5095 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