/[escript]/branches/diaplayground/cusplibrary/testing/cg_m.cu
ViewVC logotype

Contents of /branches/diaplayground/cusplibrary/testing/cg_m.cu

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4955 - (show annotations)
Tue May 20 04:33:15 2014 UTC (5 years, 3 months ago) by caltinay
File size: 2434 byte(s)
added pristine copy of cusplibrary (apache license) to be used by ripley.

1 #include <unittest/unittest.h>
2
3 #include <cusp/gallery/poisson.h>
4 #include <cusp/csr_matrix.h>
5 #include <cusp/krylov/cg_m.h>
6
7 template <class LinearOperator, class VectorType1, class VectorType2, class VectorType3>
8 void check_residuals(LinearOperator& A, VectorType1& xs, VectorType2& b, VectorType3& sigma)
9 {
10 typedef typename LinearOperator::value_type ValueType;
11 typedef typename LinearOperator::memory_space MemorySpace;
12
13 size_t N = A.num_rows;
14
15 for (size_t i = 0; i < sigma.size(); i++)
16 {
17 // compute residual = b - (A + \sigma * I) x
18 ValueType s = sigma[i];
19
20 cusp::array1d<ValueType, MemorySpace> residual(A.num_rows, 0.0f);
21
22 // TODO replace this with a array1d view of a array2d
23 cusp::array1d<ValueType, MemorySpace> x(xs.begin() + i * N, xs.begin() + (i + 1) * N);
24 cusp::multiply(A, x, residual);
25 cusp::blas::axpby(residual, x, residual, 1.0f, s);
26 cusp::blas::axpby(residual, b, residual, -1.0f, 1.0f);
27
28 ASSERT_EQUAL(cusp::blas::nrm2(residual) < 1e-4 * cusp::blas::nrm2(b), true);
29
30 //std::cout << "Residual for sigma = " << s << " is " << cusp::blas::nrm2(residual) << std::endl;
31 }
32 } // end check_residuals
33
34
35 template <class MemorySpace>
36 void TestConjugateGradientM(void)
37 {
38 // which floating point type to use
39 typedef float ValueType;
40
41 // create an empty sparse matrix structure (HYB format)
42 cusp::csr_matrix<int, ValueType, MemorySpace> A;
43
44 // create a 2d Poisson problem on a 10x10 mesh
45 cusp::gallery::poisson5pt(A, 10, 10);
46
47 // allocate storage for solution (x) and right hand side (b)
48 size_t N_s = 4;
49 cusp::array1d<ValueType, MemorySpace> x(A.num_rows*N_s, ValueType(0));
50 cusp::array1d<ValueType, MemorySpace> b(A.num_rows, ValueType(1));
51
52 // set sigma values
53 cusp::array1d<ValueType, MemorySpace> sigma(N_s);
54 sigma[0] = ValueType(0.1);
55 sigma[1] = ValueType(0.5);
56 sigma[2] = ValueType(1.0);
57 sigma[3] = ValueType(5.0);
58
59 // set stopping criteria:
60 // iteration_limit = 100
61 // relative_tolerance = 1e-6
62 cusp::default_monitor<ValueType> monitor(b, 100, 1e-6);
63
64 // solve the linear systems (A + \sigma_i * I) * x = b for each
65 // sigma_i with the Conjugate Gradient method
66 cusp::krylov::cg_m(A, x, b, sigma, monitor);
67
68 check_residuals(A, x, b, sigma);
69 }
70 DECLARE_HOST_DEVICE_UNITTEST(TestConjugateGradientM);
71

  ViewVC Help
Powered by ViewVC 1.1.26