/[escript]/trunk/cusplibrary/cusp/gallery/random.h
ViewVC logotype

Contents of /trunk/cusplibrary/cusp/gallery/random.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5148 - (show annotations)
Mon Sep 15 01:25:23 2014 UTC (5 years, 10 months ago) by caltinay
File MIME type: text/plain
File size: 4153 byte(s)
Merging ripley diagonal storage + CUDA support into trunk.
Options file version has been incremented due to new options
'cuda' and 'nvccflags'.

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 random.h
18 * \brief Random matrix generators
19 */
20
21 #pragma once
22
23 #include <cusp/detail/config.h>
24
25 #include <cusp/cds_matrix.h>
26 #include <cusp/coo_matrix.h>
27 #include <thrust/unique.h>
28 #include <thrust/sort.h>
29
30 #include <stdlib.h> // XXX remove when we switch RNGs
31
32 namespace cusp
33 {
34 namespace gallery
35 {
36 /*! \addtogroup gallery Matrix Gallery
37 * \ingroup gallery
38 * \{
39 */
40
41 // TODO use thrust RNGs, add seed parameter defaulting to num_rows ^ num_cols ^ num_samples
42 // TODO document
43 template <class MatrixType>
44 void random(size_t num_rows, size_t num_cols, size_t num_samples, MatrixType& output)
45 {
46 typedef typename MatrixType::index_type IndexType;
47 typedef typename MatrixType::value_type ValueType;
48
49 cusp::coo_matrix<IndexType,ValueType,cusp::host_memory> coo(num_rows, num_cols, num_samples);
50
51 srand(num_rows ^ num_cols ^ num_samples);
52
53 for(size_t n = 0; n < num_samples; n++)
54 {
55 coo.row_indices[n] = rand() % num_rows;
56 coo.column_indices[n] = rand() % num_cols;
57 coo.values[n] = ValueType(1);
58 }
59
60 // sort indices by (row,column)
61 coo.sort_by_row_and_column();
62
63 size_t num_entries = thrust::unique(thrust::make_zip_iterator(thrust::make_tuple(coo.row_indices.begin(), coo.column_indices.begin())),
64 thrust::make_zip_iterator(thrust::make_tuple(coo.row_indices.end(), coo.column_indices.end())))
65 - thrust::make_zip_iterator(thrust::make_tuple(coo.row_indices.begin(), coo.column_indices.begin()));
66
67 coo.resize(num_rows, num_cols, num_entries);
68
69 output = coo;
70 }
71
72 template <class MatrixType>
73 void randomblock(size_t num_rows, size_t num_diagonals, size_t block_size, MatrixType& output)
74 {
75 typedef typename MatrixType::index_type IndexType;
76 typedef typename MatrixType::value_type ValueType;
77
78 if (num_rows % block_size != 0)
79 throw cusp::invalid_input_exception("number of rows must be a multiple of block size!");
80
81 cusp::cds_matrix<IndexType,ValueType,cusp::host_memory> cds(num_rows, 0, num_diagonals, block_size);
82
83 srand(num_rows ^ num_diagonals);
84
85 // instead of entirely random, let's try and even out the number of
86 // subdiagonals and superdiagonals
87 const size_t max_offset = num_rows/32 - 1;
88 for(size_t n = 0; n < num_diagonals/2; n++)
89 {
90 const int offset = 1 + rand() % (max_offset-1);
91 cds.diagonal_offsets[2*n] = -offset;
92 cds.diagonal_offsets[2*n+1] = offset;
93 }
94 // for odd number of diagonals add main diagonal
95 if (num_diagonals%2 == 1)
96 cds.diagonal_offsets[num_diagonals-1]=0;
97
98 std::sort(cds.diagonal_offsets.begin(),cds.diagonal_offsets.end());
99
100 size_t num_entries = 0;
101
102 for(size_t n = 0; n < num_diagonals; n++)
103 {
104 const int offset = cds.diagonal_offsets[n];
105 const size_t num_blocks = num_rows/block_size-std::abs(offset);
106 num_entries += num_blocks*block_size*block_size;
107 const size_t first = std::max(0, -offset*(int)block_size);
108 for(size_t block = 0; block < num_blocks; block++)
109 {
110 for(size_t i = 0; i < block_size; i++)
111 {
112 for(size_t j = 0; j < block_size; j++)
113 {
114 cds.values(first+block*block_size+i, n*block_size+j) = ValueType(1);
115 }
116 }
117 }
118 }
119
120 cds.resize(num_rows, num_entries, block_size, num_diagonals);
121
122 output = cds;
123 }
124 /*! \}
125 */
126
127 } // end namespace gallery
128 } // end namespace cusp
129

  ViewVC Help
Powered by ViewVC 1.1.26