/[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 5095 - (hide annotations)
Fri Jul 11 05:54:25 2014 UTC (6 years ago) by caltinay
Original Path: branches/diaplayground/cusplibrary/cusp/cds_matrix.h
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 caltinay 5095 /*
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