ddemidov / vexcl
1
#define BOOST_TEST_MODULE SparseMatrices
2
#include <boost/test/unit_test.hpp>
3
#include <vexcl/vector.hpp>
4
#include <vexcl/sparse/csr.hpp>
5
#include <vexcl/sparse/ell.hpp>
6
#include <vexcl/sparse/matrix.hpp>
7
#include <vexcl/sparse/distributed.hpp>
8

9
typedef std::array<std::array<double, 2>, 2> matrix_value;
10
typedef std::array<double, 2> vector_value;
11

12
matrix_value mconst(double c) {
13 6
    matrix_value a = {c, c, c, c};
14 6
    return a;
15
}
16

17
vector_value vconst(double c) {
18
    vector_value v = {c, c};
19 6
    return v;
20
}
21

22
namespace vex {
23

24
template <> struct is_cl_native< matrix_value > : std::true_type {};
25
template <> struct is_cl_native< vector_value > : std::true_type {};
26

27
template <> struct type_name_impl<matrix_value> {
28 6
    static std::string get() { return "double4"; }
29
};
30

31
template <> struct type_name_impl<vector_value> {
32 6
    static std::string get() { return "double2"; }
33
};
34

35
namespace sparse {
36

37
template <>
38
struct rhs_of< matrix_value > {
39
    typedef vector_value type;
40
};
41

42
template <>
43
struct spmv_ops_impl<matrix_value, vector_value> {
44
    static void decl_accum_var(backend::source_generator &src, const std::string &name)
45
    {
46 6
        src.new_line() << type_name<vector_value>() << " " << name << " = {0,0};";
47
    }
48

49
    static void append_product(backend::source_generator &src,
50
            const std::string &sum, const std::string &mat_val, const std::string &vec_val)
51
    {
52 6
        src.open("{");
53 6
        src.new_line() << type_name<vector_value>() << " v = " << vec_val << ";";
54 6
        src.new_line() << sum << ".x += " << mat_val << ".x * v.x + " << mat_val << ".y * v.y;";
55 6
        src.new_line() << sum << ".y += " << mat_val << ".z * v.x + " << mat_val << ".w * v.y;";
56 6
        src.close("}");
57
    }
58
};
59

60
} // namespace sparse
61
} // namespace vex
62

63
#include "context_setup.hpp"
64
#include "random_matrix.hpp"
65

66 6
BOOST_AUTO_TEST_CASE(csr)
67
{
68
    const size_t n = 1024;
69

70 6
    std::vector<vex::command_queue> q(1, ctx.queue(0));
71

72 6
    std::vector<int>    row;
73 6
    std::vector<int>    col;
74 6
    std::vector<double> val;
75

76 6
    random_matrix(n, n, 16, row, col, val);
77

78 6
    std::vector<double> x = random_vector<double>(n);
79

80 6
    vex::sparse::csr<double> A(q, n, n, row, col, val);
81 6
    vex::vector<double> X(q, x);
82 6
    vex::vector<double> Y(q, n);
83

84 6
    Y = A * X;
85

86 6
    check_sample(Y, [&](size_t idx, double a) {
87 6
            double sum = 0;
88 6
            for(int j = row[idx]; j < row[idx + 1]; j++)
89 6
                sum += val[j] * x[col[j]];
90

91 6
            BOOST_CHECK_CLOSE(a, sum, 1e-8);
92 6
            });
93
}
94

95 6
BOOST_AUTO_TEST_CASE(ell)
96
{
97
    const size_t n = 1024;
98

99 6
    std::vector<vex::command_queue> q(1, ctx.queue(0));
100

101 6
    std::vector<int>    row;
102 6
    std::vector<int>    col;
103 6
    std::vector<double> val;
104

105 6
    random_matrix(n, n, 16, row, col, val);
106

107 6
    std::vector<double> x = random_vector<double>(n);
108

109 6
    vex::sparse::ell<double> A(q, n, n, row, col, val);
110 6
    vex::vector<double> X(q, x);
111 6
    vex::vector<double> Y(q, n);
112

113 6
    Y = A * X;
114

115 6
    check_sample(Y, [&](size_t idx, double a) {
116 6
            double sum = 0;
117 6
            for(int j = row[idx]; j < row[idx + 1]; j++)
118 6
                sum += val[j] * x[col[j]];
119

120 6
            BOOST_CHECK_CLOSE(a, sum, 1e-8);
121 6
            });
122
}
123

124 6
BOOST_AUTO_TEST_CASE(matrix)
125
{
126
    const size_t n = 1024;
127

128 6
    std::vector<vex::command_queue> q(1, ctx.queue(0));
129

130 6
    std::vector<int>    row;
131 6
    std::vector<int>    col;
132 6
    std::vector<double> val;
133

134 6
    random_matrix(n, n, 16, row, col, val);
135

136 6
    std::vector<double> x = random_vector<double>(n);
137

138 6
    vex::sparse::matrix<double> A(q, n, n, row, col, val);
139 6
    vex::vector<double> X(q, x);
140 6
    vex::vector<double> Y(q, n);
141

142 6
    Y = A * X;
143

144 6
    check_sample(Y, [&](size_t idx, double a) {
145 6
            double sum = 0;
146 6
            for(int j = row[idx]; j < row[idx + 1]; j++)
147 6
                sum += val[j] * x[col[j]];
148

149 6
            BOOST_CHECK_CLOSE(a, sum, 1e-8);
150 6
            });
151
}
152

153 6
BOOST_AUTO_TEST_CASE(distributed)
154
{
155
    const int n = 1024;
156

157 6
    std::vector<int>    ptr;
158 6
    std::vector<int>    col;
159 6
    std::vector<double> val;
160

161 6
    ptr.push_back(0);
162 6
    for(int i = 0; i < n; ++i) {
163 6
        if (i > 0) {
164 6
            col.push_back(i-1);
165 6
            val.push_back(-1);
166
        }
167 6
        col.push_back(i);
168 6
        val.push_back(2);
169 6
        if (i + 1 < n) {
170 6
            col.push_back(i+1);
171 6
            val.push_back(-1);
172
        }
173

174 6
        ptr.push_back(static_cast<int>(col.size()));
175
    }
176

177 6
    vex::sparse::distributed<vex::sparse::ell<double>> A(ctx, n, n, ptr, col, val);
178

179 6
    std::vector<double> x = random_vector<double>(n);
180 6
    vex::vector<double> X(ctx, x);
181 6
    vex::vector<double> Y(ctx, n);
182

183 6
    Y = A * X;
184

185 6
    for(int i = 0; i < n; ++i) {
186 6
        double y = Y[i];
187 6
        double sum = 0;
188 6
        for(int j = ptr[i]; j < ptr[i + 1]; j++)
189 6
            sum += val[j] * x[col[j]];
190

191 6
        BOOST_CHECK_CLOSE(y, sum, 1e-8);
192
    }
193
}
194

195 6
BOOST_AUTO_TEST_CASE(distributed_single)
196
{
197 6
    std::vector<vex::command_queue> q(1, ctx.queue(0));
198

199
    const int n = 1024;
200

201 6
    std::vector<int>    ptr;
202 6
    std::vector<int>    col;
203 6
    std::vector<double> val;
204

205 6
    ptr.push_back(0);
206 6
    for(int i = 0; i < n; ++i) {
207 6
        if (i > 0) {
208 6
            col.push_back(i-1);
209 6
            val.push_back(-1);
210
        }
211 6
        col.push_back(i);
212 6
        val.push_back(2);
213 6
        if (i + 1 < n) {
214 6
            col.push_back(i+1);
215 6
            val.push_back(-1);
216
        }
217

218 6
        ptr.push_back(static_cast<int>(col.size()));
219
    }
220

221 6
    vex::sparse::distributed<vex::sparse::ell<double>> A(q, n, n, ptr, col, val);
222

223 6
    std::vector<double> x = random_vector<double>(n);
224 6
    vex::vector<double> X(q, x);
225 6
    vex::vector<double> Y(q, n);
226

227 6
    Y = A * X;
228

229 6
    for(int i = 0; i < n; ++i) {
230 6
        double y = Y[i];
231 6
        double sum = 0;
232 6
        for(int j = ptr[i]; j < ptr[i + 1]; j++)
233 6
            sum += val[j] * x[col[j]];
234

235 6
        BOOST_CHECK_CLOSE(y, sum, 1e-8);
236
    }
237
}
238

239 6
BOOST_AUTO_TEST_CASE(custom_values)
240
{
241
    const int n = 1024;
242 6
    std::vector<vex::command_queue> q(1, ctx.queue(0));
243

244 6
    std::vector<int> ptr, col;
245 6
    std::vector<matrix_value> val;
246

247 6
    ptr.push_back(0);
248 6
    for(int i = 0; i < n; ++i) {
249 6
        if (i > 0) {
250 6
            col.push_back(i-1);
251 6
            val.push_back(mconst(-1));
252
        }
253 6
        col.push_back(i);
254 6
        val.push_back(mconst(2));
255 6
        if (i + 1 < n) {
256 6
            col.push_back(i+1);
257 6
            val.push_back(mconst(-1));
258
        }
259

260 6
        ptr.push_back(static_cast<int>(col.size()));
261
    }
262

263 6
    vex::sparse::matrix<matrix_value> A(q, n, n, ptr, col, val);
264

265 6
    std::vector<vector_value> x(n, vconst(1));
266 6
    vex::vector<vector_value> X(q, x);
267 6
    vex::vector<vector_value> Y(q, n);
268

269 6
    Y = A * X;
270

271 6
    for(int i = 0; i < n; ++i) {
272 6
        vector_value y = Y[i];
273 6
        vector_value sum = {0,0};
274 6
        for(int j = ptr[i]; j < ptr[i + 1]; j++) {
275 6
            sum[0] += val[j][0][0] * x[col[j]][0] + val[j][0][1] * x[col[j]][1];
276 6
            sum[1] += val[j][1][0] * x[col[j]][0] + val[j][1][1] * x[col[j]][1];
277
        }
278

279 6
        BOOST_CHECK_CLOSE(y[0], sum[0], 1e-8);
280 6
        BOOST_CHECK_CLOSE(y[1], sum[1], 1e-8);
281
    }
282
}
283

284

285 6
BOOST_AUTO_TEST_SUITE_END()

Read our documentation on viewing source code .

Loading