ddemidov / vexcl
1
#define BOOST_TEST_MODULE SparseMatrixVectorProduct
2
#include <boost/test/unit_test.hpp>
3
#include <vexcl/vector.hpp>
4
#include <vexcl/multivector.hpp>
5
#include <vexcl/spmat.hpp>
6
#include <vexcl/function.hpp>
7
#include "context_setup.hpp"
8
#include "random_matrix.hpp"
9

10 6
BOOST_AUTO_TEST_CASE(vector_product)
11
{
12
    const size_t n = 1024;
13

14 6
    std::vector<size_t> row;
15 6
    std::vector<size_t> col;
16 6
    std::vector<double> val;
17

18 6
    random_matrix(n, n, 16, row, col, val);
19

20 6
    std::vector<double> x = random_vector<double>(n);
21

22 6
    vex::SpMat <double> A(ctx, n, n, row.data(), col.data(), val.data());
23 6
    vex::vector<double> X(ctx, x);
24 6
    vex::vector<double> Y(ctx, n);
25

26 6
    Y = A * X;
27

28 6
    check_sample(Y, [&](size_t idx, double a) {
29 6
            double sum = 0;
30 6
            for(size_t j = row[idx]; j < row[idx + 1]; j++)
31 6
                sum += val[j] * x[col[j]];
32

33 6
            BOOST_CHECK_CLOSE(a, sum, 1e-8);
34 6
            });
35

36 6
    Y -= A * X;
37

38 6
    check_sample(Y, [&](size_t, double a) { BOOST_CHECK_SMALL(a, 1e-8); });
39

40 6
    Y += 42 * (A * X);
41

42 6
    check_sample(Y, [&](size_t idx, double a) {
43
            double sum = 0;
44 6
            for(size_t j = row[idx]; j < row[idx + 1]; j++)
45 6
                sum += val[j] * x[col[j]];
46

47 6
            BOOST_CHECK_CLOSE(a, 42 * sum, 1e-8);
48 6
            });
49

50 6
    Y = X + A * X;
51

52 6
    check_sample(Y, [&](size_t idx, double a) {
53
            double sum = 0;
54 6
            for(size_t j = row[idx]; j < row[idx + 1]; j++)
55 6
                sum += val[j] * x[col[j]];
56

57 6
            BOOST_CHECK_CLOSE(a, x[idx] + sum, 1e-8);
58 6
            });
59
}
60

61 6
BOOST_AUTO_TEST_CASE(non_square_matrix)
62
{
63
    const size_t n = 1024;
64
    const size_t m = 2 * n;
65

66 6
    std::vector<size_t> row;
67 6
    std::vector<size_t> col;
68 6
    std::vector<double> val;
69

70 6
    random_matrix(n, m, 16, row, col, val);
71

72 6
    std::vector<double> x = random_vector<double>(m);
73

74 6
    vex::SpMat <double> A(ctx, n, m, row.data(), col.data(), val.data());
75 6
    vex::vector<double> X(ctx, x);
76 6
    vex::vector<double> Y(ctx, n);
77

78 6
    Y = A * X;
79

80 6
    check_sample(Y, [&](size_t idx, double a) {
81 6
            double sum = 0;
82 6
            for(size_t j = row[idx]; j < row[idx + 1]; j++)
83 6
                sum += val[j] * x[col[j]];
84

85 6
            BOOST_CHECK_CLOSE(a, sum, 1e-8);
86 6
            });
87
}
88

89 6
BOOST_AUTO_TEST_CASE(non_default_types)
90
{
91
    const size_t n = 1024;
92

93 6
    std::vector<unsigned> row;
94 6
    std::vector<int>      col;
95 6
    std::vector<double>   val;
96

97 6
    random_matrix(n, n, 16, row, col, val);
98

99 6
    std::vector<double> x = random_vector<double>(n);
100

101 6
    vex::SpMat <double, int, unsigned> A(ctx, n, n, row.data(), col.data(), val.data());
102 6
    vex::vector<double> X(ctx, x);
103 6
    vex::vector<double> Y(ctx, n);
104

105 6
    Y = A * X;
106

107 6
    check_sample(Y, [&](size_t idx, double a) {
108 6
            double sum = 0;
109 6
            for(size_t j = row[idx]; j < row[idx + 1]; j++)
110 6
                sum += val[j] * x[col[j]];
111

112 6
            BOOST_CHECK_CLOSE(a, sum, 1e-8);
113 6
            });
114
}
115

116 6
BOOST_AUTO_TEST_CASE(empty_rows)
117
{
118
    const size_t n = 1024;
119
    const size_t non_empty_part = 256;
120

121 6
    std::vector<size_t> row;
122 6
    std::vector<size_t> col;
123 6
    std::vector<double> val;
124

125 6
    row.reserve(n + 1);
126

127 6
    random_matrix(non_empty_part, n, 16, row, col, val);
128

129 6
    while(row.size() < n + 1) row.push_back(col.size());
130

131 6
    std::vector<double> x = random_vector<double>(n);
132

133 6
    vex::SpMat <double> A(ctx, n, n, row.data(), col.data(), val.data());
134 6
    vex::vector<double> X(ctx, x);
135 6
    vex::vector<double> Y(ctx, n);
136

137 6
    Y = A * X;
138

139 6
    check_sample(Y, [&](size_t idx, double a) {
140 6
            double sum = 0;
141 6
            for(size_t j = row[idx]; j < row[idx + 1]; j++)
142 6
                sum += val[j] * x[col[j]];
143

144 6
            BOOST_CHECK_CLOSE(a, sum, 1e-8);
145 6
            });
146
}
147

148 6
BOOST_AUTO_TEST_CASE(ccsr_vector_product)
149
{
150
    const size_t n = 32;
151
    const double h2i = (n - 1) * (n - 1);
152

153 6
    std::vector<size_t> idx;
154 6
    std::vector<size_t> row(3);
155 6
    std::vector<int>    col(8);
156 6
    std::vector<double> val(8);
157

158 6
    idx.reserve(n * n * n);
159

160 6
    row[0] = 0;
161 6
    row[1] = 1;
162 6
    row[2] = 8;
163

164 6
    col[0] = 0;
165 6
    val[0] = 1;
166

167 6
    col[1] = -static_cast<int>(n * n);
168 6
    col[2] = -static_cast<int>(n);
169 6
    col[3] =    -1;
170 6
    col[4] =     0;
171 6
    col[5] =     1;
172 6
    col[6] =     n;
173 6
    col[7] =  (n * n);
174

175 6
    val[1] = -h2i;
176 6
    val[2] = -h2i;
177 6
    val[3] = -h2i;
178 6
    val[4] =  h2i * 6;
179 6
    val[5] = -h2i;
180 6
    val[6] = -h2i;
181 6
    val[7] = -h2i;
182

183 6
    for(size_t k = 0; k < n; k++) {
184 6
        for(size_t j = 0; j < n; j++) {
185 6
            for(size_t i = 0; i < n; i++) {
186
                if (
187 6
                        i == 0 || i == (n - 1) ||
188 6
                        j == 0 || j == (n - 1) ||
189 6
                        k == 0 || k == (n - 1)
190
                   )
191
                {
192 6
                    idx.push_back(0);
193 6
                } else {
194 6
                    idx.push_back(1);
195
                }
196
            }
197
        }
198
    }
199

200 6
    std::vector<double> x = random_vector<double>(n * n * n);
201

202 6
    std::vector<vex::command_queue> queue(1, ctx.queue(0));
203

204 6
    vex::SpMatCCSR<double,int> A(queue[0], x.size(), row.size() - 1,
205 6
            idx.data(), row.data(), col.data(), val.data());
206

207 6
    vex::vector<double> X(queue, x);
208 6
    vex::vector<double> Y(queue, x.size());
209

210 6
    Y = A * X;
211

212 6
    check_sample(Y, [&](size_t ii, double a) {
213 6
            double sum = 0;
214 6
            size_t i = idx[ii];
215 6
            for(size_t j = row[i]; j < row[i + 1]; j++)
216 6
                sum += val[j] * x[ii + col[j]];
217

218 6
            BOOST_CHECK_CLOSE(a, sum, 1e-8);
219 6
            });
220

221 6
    Y = X + A * X;
222

223 6
    check_sample(Y, [&](size_t ii, double a) {
224
            double sum = 0;
225 6
            size_t i = idx[ii];
226 6
            for(size_t j = row[i]; j < row[i + 1]; j++)
227 6
                sum += val[j] * x[ii + col[j]];
228

229 6
            BOOST_CHECK_CLOSE(a, x[ii] + sum, 1e-8);
230 6
            });
231
}
232

233 6
BOOST_AUTO_TEST_CASE(inline_spmv)
234
{
235
    const size_t n = 1024;
236

237 6
    std::vector<vex::command_queue> queue(1, ctx.queue(0));
238

239 6
    std::vector<size_t> row;
240 6
    std::vector<size_t> col;
241 6
    std::vector<double> val;
242

243 6
    random_matrix(n, n, 16, row, col, val);
244

245 6
    std::vector<double> x = random_vector<double>(n);
246

247 6
    vex::SpMat <double> A(queue, n, n, row.data(), col.data(), val.data());
248 6
    vex::vector<double> X(queue, x);
249 6
    vex::vector<double> Y(queue, n);
250

251 6
    Y = sin(vex::make_inline(A * X));
252

253 6
    check_sample(Y, [&](size_t idx, double a) {
254
            double sum = 0;
255 6
            for(size_t j = row[idx]; j < row[idx + 1]; j++)
256 6
                sum += val[j] * x[col[j]];
257

258 6
            BOOST_CHECK_CLOSE(a, sin(sum), 1e-8);
259 6
            });
260
}
261

262 6
BOOST_AUTO_TEST_CASE(multivector_product)
263
{
264
    const size_t n = 1024;
265
    const size_t m = 2;
266

267
    typedef std::array<double, m> elem_t;
268

269 6
    std::vector<size_t> row;
270 6
    std::vector<size_t> col;
271 6
    std::vector<double> val;
272

273 6
    random_matrix(n, n, 16, row, col, val);
274

275 6
    std::vector<double> x = random_vector<double>(n * m);
276

277 6
    vex::SpMat <double> A(ctx, n, n, row.data(), col.data(), val.data());
278

279 6
    vex::multivector<double,m> X(ctx, x);
280 6
    vex::multivector<double,m> Y(ctx, n);
281

282 6
    Y = A * X;
283

284 6
    check_sample(Y, [&](size_t idx, elem_t a) {
285 6
            double sum[] = {0, 0};
286 6
            for(size_t j = row[idx]; j < row[idx + 1]; j++) {
287 6
                sum[0] += val[j] * x[0 + col[j]];
288 6
                sum[1] += val[j] * x[n + col[j]];
289
            }
290

291 6
            BOOST_CHECK_CLOSE(a[0], sum[0], 1e-8);
292 6
            BOOST_CHECK_CLOSE(a[1], sum[1], 1e-8);
293 6
            });
294

295 6
    Y = X + A * X;
296

297 6
    check_sample(Y, [&](size_t idx, elem_t a) {
298
            double sum[] = {0, 0};
299 6
            for(size_t j = row[idx]; j < row[idx + 1]; j++) {
300 6
                sum[0] += val[j] * x[0 + col[j]];
301 6
                sum[1] += val[j] * x[n + col[j]];
302
            }
303

304 6
            BOOST_CHECK_CLOSE(a[0], x[0 + idx] + sum[0], 1e-8);
305 6
            BOOST_CHECK_CLOSE(a[1], x[n + idx] + sum[1], 1e-8);
306 6
            });
307
}
308

309 6
BOOST_AUTO_TEST_CASE(inline_multivector_product)
310
{
311
    const size_t n = 1024;
312
    const size_t m = 2;
313

314
    typedef std::array<double, m> elem_t;
315

316 6
    std::vector<vex::command_queue> queue(1, ctx.queue(0));
317

318 6
    std::vector<size_t> row;
319 6
    std::vector<size_t> col;
320 6
    std::vector<double> val;
321

322 6
    random_matrix(n, n, 16, row, col, val);
323

324 6
    std::vector<double> x = random_vector<double>(n * m);
325

326 6
    vex::SpMat <double> A(queue, n, n, row.data(), col.data(), val.data());
327

328 6
    vex::multivector<double,m> X(queue, x);
329 6
    vex::multivector<double,m> Y(queue, n);
330

331 6
    Y = cos(vex::make_inline(A * X));
332

333 6
    check_sample(Y, [&](size_t idx, elem_t a) {
334
            double sum[] = {0, 0};
335 6
            for(size_t j = row[idx]; j < row[idx + 1]; j++) {
336 6
                sum[0] += val[j] * x[0 + col[j]];
337 6
                sum[1] += val[j] * x[n + col[j]];
338
            }
339

340 6
            BOOST_CHECK_CLOSE(a[0], cos(sum[0]), 1e-8);
341 6
            BOOST_CHECK_CLOSE(a[1], cos(sum[1]), 1e-8);
342 6
            });
343
}
344

345 6
BOOST_AUTO_TEST_CASE(ccsr_multivector_product)
346
{
347
    const size_t n = 32;
348
    const size_t N = n * n * n;
349
    const double h2i = (n - 1) * (n - 1);
350

351
    typedef std::array<double, 2> elem_t;
352

353 6
    std::vector<size_t> idx;
354 6
    std::vector<size_t> row(3);
355 6
    std::vector<int>    col(8);
356 6
    std::vector<double> val(8);
357

358 6
    idx.reserve(N);
359

360 6
    row[0] = 0;
361 6
    row[1] = 1;
362 6
    row[2] = 8;
363

364 6
    col[0] = 0;
365 6
    val[0] = 1;
366

367 6
    col[1] = -static_cast<int>(n * n);
368 6
    col[2] = -static_cast<int>(n);
369 6
    col[3] =    -1;
370 6
    col[4] =     0;
371 6
    col[5] =     1;
372 6
    col[6] =     n;
373 6
    col[7] =  (n * n);
374

375 6
    val[1] = -h2i;
376 6
    val[2] = -h2i;
377 6
    val[3] = -h2i;
378 6
    val[4] =  h2i * 6;
379 6
    val[5] = -h2i;
380 6
    val[6] = -h2i;
381 6
    val[7] = -h2i;
382

383 6
    for(size_t k = 0; k < n; k++) {
384 6
        for(size_t j = 0; j < n; j++) {
385 6
            for(size_t i = 0; i < n; i++) {
386
                if (
387 6
                        i == 0 || i + 1 == n ||
388 6
                        j == 0 || j + 1 == n ||
389 6
                        k == 0 || k + 1 == n
390
                   )
391
                {
392 6
                    idx.push_back(0);
393 6
                } else {
394 6
                    idx.push_back(1);
395
                }
396
            }
397
        }
398
    }
399

400 6
    std::vector<double> x = random_vector<double>(N * 2);
401

402 6
    std::vector<vex::command_queue> queue(1, ctx.queue(0));
403

404 6
    vex::SpMatCCSR<double,int> A(queue[0], N, row.size() - 1,
405 6
            idx.data(), row.data(), col.data(), val.data());
406

407 6
    vex::multivector<double,2> X(queue, x);
408 6
    vex::multivector<double,2> Y(queue, N);
409

410 6
    Y = A * X;
411

412 6
    check_sample(Y, [&](size_t ii, elem_t a) {
413 6
            double sum[] = {0, 0};
414 6
            size_t i = idx[ii];
415 6
            for(size_t j = row[i]; j < row[i + 1]; j++) {
416 6
                sum[0] += val[j] * x[0 + ii + col[j]];
417 6
                sum[1] += val[j] * x[N + ii + col[j]];
418
            }
419

420 6
            BOOST_CHECK_CLOSE(a[0], sum[0], 1e-8);
421 6
            BOOST_CHECK_CLOSE(a[1], sum[1], 1e-8);
422 6
            });
423

424 6
    Y = X + A * X;
425

426 6
    check_sample(Y, [&](size_t ii, elem_t a) {
427
            double sum[] = {0, 0};
428 6
            size_t i = idx[ii];
429 6
            for(size_t j = row[i]; j < row[i + 1]; j++) {
430 6
                sum[0] += val[j] * x[0 + ii + col[j]];
431 6
                sum[1] += val[j] * x[N + ii + col[j]];
432
            }
433

434 6
            BOOST_CHECK_CLOSE(a[0], x[0 + ii] + sum[0], 1e-8);
435 6
            BOOST_CHECK_CLOSE(a[1], x[N + ii] + sum[1], 1e-8);
436 6
            });
437
}
438

439
#if defined(VEXCL_BACKEND_OPENCL) || defined(VEXCL_BACKEND_COMPUTE)
440 4
BOOST_AUTO_TEST_CASE(vector_valued_matrix)
441
{
442
    const size_t n = 1024;
443

444 4
    std::vector<size_t>     row;
445 4
    std::vector<size_t>     col;
446 4
    std::vector<cl_double2> val;
447

448 4
    random_matrix(n, n, 16, row, col, val);
449

450 4
    std::vector<cl_double2> x = random_vector<cl_double2>(n);
451

452 4
    vex::SpMat <cl_double2> A(ctx, n, n, row.data(), col.data(), val.data());
453 4
    vex::vector<cl_double2> X(ctx, x);
454 4
    vex::vector<cl_double2> Y(ctx, n);
455

456 4
    Y = A * X;
457

458 4
    check_sample(Y, [&](size_t idx, cl_double2 a) {
459 4
            cl_double2 sum = {{0, 0}};
460 4
            for(size_t j = row[idx]; j < row[idx + 1]; j++)
461 4
                sum += val[j] * x[col[j]];
462

463 4
            BOOST_CHECK_CLOSE(a.s[0], sum.s[0], 1e-8);
464 4
            BOOST_CHECK_CLOSE(a.s[1], sum.s[1], 1e-8);
465 4
            });
466
}
467

468
cl_double2 make_d2(double v) {
469
    cl_double2 d2 = {{v, v}};
470 4
    return d2;
471
}
472

473 4
BOOST_AUTO_TEST_CASE(vector_valued_ccsr_matrix)
474
{
475
    const size_t n = 32;
476
    const double h2i = (n - 1) * (n - 1);
477

478 4
    std::vector<size_t> idx;
479 4
    std::vector<size_t> row(3);
480 4
    std::vector<int>    col(8);
481 4
    std::vector<cl_double2> val(8);
482

483 4
    idx.reserve(n * n * n);
484

485 4
    row[0] = 0;
486 4
    row[1] = 1;
487 4
    row[2] = 8;
488

489 4
    col[0] = 0;
490 4
    val[0] = make_d2(1);
491

492 4
    col[1] = -static_cast<int>(n * n);
493 4
    col[2] = -static_cast<int>(n);
494 4
    col[3] =    -1;
495 4
    col[4] =     0;
496 4
    col[5] =     1;
497 4
    col[6] =     n;
498 4
    col[7] =  (n * n);
499

500 4
    val[1] = make_d2(-h2i);
501 4
    val[2] = make_d2(-h2i);
502 4
    val[3] = make_d2(-h2i);
503 4
    val[4] = make_d2( h2i * 6);
504 4
    val[5] = make_d2(-h2i);
505 4
    val[6] = make_d2(-h2i);
506 4
    val[7] = make_d2(-h2i);
507

508 4
    for(size_t k = 0; k < n; k++) {
509 4
        for(size_t j = 0; j < n; j++) {
510 4
            for(size_t i = 0; i < n; i++) {
511
                if (
512 4
                        i == 0 || i == (n - 1) ||
513 4
                        j == 0 || j == (n - 1) ||
514 4
                        k == 0 || k == (n - 1)
515
                   )
516
                {
517 4
                    idx.push_back(0);
518 4
                } else {
519 4
                    idx.push_back(1);
520
                }
521
            }
522
        }
523
    }
524

525 4
    std::vector<cl_double2> x = random_vector<cl_double2>(n * n * n);
526

527 4
    std::vector<vex::command_queue> queue(1, ctx.queue(0));
528

529 4
    vex::SpMatCCSR<cl_double2,int> A(queue[0], x.size(), row.size() - 1,
530 4
            idx.data(), row.data(), col.data(), val.data());
531

532 4
    vex::vector<cl_double2> X(queue, x);
533 4
    vex::vector<cl_double2> Y(queue, x.size());
534

535 4
    Y = A * X;
536

537 4
    check_sample(Y, [&](size_t ii, cl_double2 a) {
538 4
            cl_double2 sum = {{0, 0}};
539 4
            size_t i = idx[ii];
540 4
            for(size_t j = row[i]; j < row[i + 1]; j++)
541 4
                sum += val[j] * x[ii + col[j]];
542

543 4
            BOOST_CHECK_CLOSE(a.s[0], sum.s[0], 1e-8);
544 4
            BOOST_CHECK_CLOSE(a.s[1], sum.s[1], 1e-8);
545 4
            });
546
}
547
#endif
548

549 6
BOOST_AUTO_TEST_SUITE_END()

Read our documentation on viewing source code .

Loading