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

Contents of /branches/diaplayground/cusplibrary/testing/ainv.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: 7396 byte(s)
added pristine copy of cusplibrary (apache license) to be used by ripley.

1 #include <unittest/unittest.h>
2
3 #include <cusp/precond/ainv.h>
4
5 #include <cusp/gallery/poisson.h>
6 #include <cusp/krylov/cg.h>
7
8 inline
9 __host__ __device__
10 unsigned int hash32(unsigned int a)
11 {
12 a = (a + 0x7ed55d16) + (a << 12);
13 a = (a ^ 0xc761c23c) ^ (a >> 19);
14 a = (a + 0x165667b1) + (a << 5);
15 a = (a + 0xd3a2646c) ^ (a << 9);
16 a = (a + 0xfd7046c5) + (a << 3);
17 a = (a ^ 0xb55a4f09) ^ (a >> 16);
18 return a;
19 }
20
21 struct hash_01
22 {
23 __host__ __device__
24 float operator()(const unsigned int& index) const
25 {
26 return (float)(hash32(index)) / ((float)0xffffffff);
27 }
28 };
29
30 void TestAINVHeap(void)
31 {
32 hash_01 rng;
33 int seed = 100;
34
35 for (int trial = 0; trial < 10000; trial ++) {
36 cusp::precond::detail::ainv_matrix_row<int, float> row;
37
38 int size =100;
39 int i;
40 for (i=0; i < size; i++)
41 row.insert(i, rng(seed++));
42
43 for (i=0; i < size; i++)
44 row.add_to_value(i, rng(seed++) - .5);
45
46 for (i=0; i < size/2; i++)
47 row.remove_min();
48
49 for (i=0; i < size; i++)
50 row.insert(i+size, rng(seed++));
51
52 ASSERT_EQUAL(row.validate_heap(), true);
53 ASSERT_EQUAL(row.validate_backpointers(), true);
54
55 }
56 }
57 DECLARE_UNITTEST(TestAINVHeap);
58
59
60 void TestAINVFactorization(void)
61 {
62 typedef int IndexType;
63 typedef float ValueType;
64 typedef cusp::device_memory MemorySpace;
65
66 // Create 2D Poisson problem
67 cusp::csr_matrix<IndexType,ValueType,MemorySpace> A;
68 cusp::gallery::poisson5pt(A, 10, 10);
69 A.values[0] = 10;
70 int N = A.num_rows;
71
72 // factor exactly
73 cusp::precond::scaled_bridson_ainv<ValueType,MemorySpace> M(A, 0, -1);
74
75
76
77 cusp::array1d<ValueType,MemorySpace> x(N);
78 cusp::array1d<ValueType,MemorySpace> b(N, 0);
79
80 thrust::transform(thrust::counting_iterator<unsigned int>(0),
81 thrust::counting_iterator<unsigned int>(N),
82 x.begin(),
83 hash_01());
84
85 cusp::array1d<ValueType,MemorySpace> x_solve = x;
86
87 // cg should converge in 1 iteration
88 // because we're in single precision, this isn't exact, but 1e-5 tolerance should give enough leeway.
89 cusp::default_monitor<ValueType> monitor(b, 1, 0, 1e-5);
90 cusp::krylov::cg(A, x_solve, b, monitor, M);
91
92 ASSERT_EQUAL(monitor.converged(), true);
93 }
94 DECLARE_UNITTEST(TestAINVFactorization);
95
96 void TestAINVSymmetry(void)
97 {
98 typedef int IndexType;
99 typedef float ValueType;
100 typedef cusp::device_memory MemorySpace;
101
102 // Create 2D Poisson problem
103 cusp::csr_matrix<IndexType,ValueType,MemorySpace> A;
104 cusp::gallery::poisson5pt(A, 100, 100);
105 A.values[0] = 10;
106 int N = A.num_rows;
107
108 cusp::array1d<ValueType,MemorySpace> x(N);
109 cusp::array1d<ValueType,MemorySpace> b(N, 0);
110
111 thrust::transform(thrust::counting_iterator<unsigned int>(0),
112 thrust::counting_iterator<unsigned int>(N),
113 x.begin(),
114 hash_01());
115
116 ValueType nrm1, nrm2;
117
118 // test symmetric version
119 {
120 cusp::array1d<ValueType,MemorySpace> x_solve = x;
121 cusp::precond::bridson_ainv<ValueType,MemorySpace> M(A, .1);
122
123 cusp::default_monitor<ValueType> monitor(b, 125, 0, 1e-5);
124 cusp::krylov::cg(A, x_solve, b, monitor, M);
125
126 nrm1 = monitor.residual_norm();
127 ASSERT_EQUAL(monitor.converged(), true);
128 }
129
130 // test non-symmetric version
131 {
132 cusp::array1d<ValueType,MemorySpace> x_solve = x;
133 cusp::precond::nonsym_bridson_ainv<ValueType,MemorySpace> M(A, .1);
134
135 cusp::default_monitor<ValueType> monitor(b, 125, 0, 1e-5);
136 cusp::krylov::cg(A, x_solve, b, monitor, M);
137
138 nrm2 = monitor.residual_norm();
139 ASSERT_EQUAL(monitor.converged(), true);
140 }
141
142 ASSERT_EQUAL(nrm1, nrm2);
143
144 // assert they returned identical results
145 }
146 DECLARE_UNITTEST(TestAINVSymmetry);
147
148 void TestAINVConvergence(void)
149 {
150 typedef int IndexType;
151 typedef float ValueType;
152 typedef cusp::device_memory MemorySpace;
153
154 // Create 2D Poisson problem
155 cusp::csr_matrix<IndexType,ValueType,MemorySpace> A;
156 cusp::gallery::poisson5pt(A, 100, 100);
157 A.values[0] = 10;
158 int N = A.num_rows;
159
160 cusp::array1d<ValueType,MemorySpace> x(N);
161 cusp::array1d<ValueType,MemorySpace> b(N, 0);
162
163 thrust::transform(thrust::counting_iterator<unsigned int>(0),
164 thrust::counting_iterator<unsigned int>(N),
165 x.begin(),
166 hash_01());
167
168
169
170 // test drop tolerance strategy
171 {
172 cusp::array1d<ValueType,MemorySpace> x_solve = x;
173 cusp::precond::scaled_bridson_ainv<ValueType,MemorySpace> M(A, .1);
174
175 cusp::default_monitor<ValueType> monitor(b, 125, 0, 1e-5);
176 cusp::krylov::cg(A, x_solve, b, monitor, M);
177
178 ASSERT_EQUAL(monitor.converged(), true);
179 }
180
181 // test sparsity strategy
182 {
183 cusp::array1d<ValueType,MemorySpace> x_solve = x;
184 cusp::precond::scaled_bridson_ainv<ValueType,MemorySpace> M(A, 0, 10);
185
186 cusp::default_monitor<ValueType> monitor(b, 70, 0, 1e-5);
187 cusp::krylov::cg(A, x_solve, b, monitor, M);
188
189 ASSERT_EQUAL(monitor.converged(), true);
190 }
191
192 // test both
193 {
194 cusp::array1d<ValueType,MemorySpace> x_solve = x;
195 cusp::precond::scaled_bridson_ainv<ValueType,MemorySpace> M(A, .01, 4);
196
197 cusp::default_monitor<ValueType> monitor(b, 120, 0, 1e-5);
198 cusp::krylov::cg(A, x_solve, b, monitor, M);
199
200 ASSERT_EQUAL(monitor.converged(), true);
201 }
202 // test lin dropping
203 {
204 cusp::array1d<ValueType,MemorySpace> x_solve = x;
205 cusp::precond::scaled_bridson_ainv<ValueType,MemorySpace> M(A, 0, -1, true, 4);
206
207 cusp::default_monitor<ValueType> monitor(b, 120, 0, 1e-5);
208 cusp::krylov::cg(A, x_solve, b, monitor, M);
209
210 ASSERT_EQUAL(monitor.converged(), true);
211 }
212
213
214 // test drop tolerance strategy
215 {
216 cusp::array1d<ValueType,MemorySpace> x_solve = x;
217 cusp::precond::bridson_ainv<ValueType,MemorySpace> M(A, .1);
218
219 cusp::default_monitor<ValueType> monitor(b, 125, 0, 1e-5);
220 cusp::krylov::cg(A, x_solve, b, monitor, M);
221
222 ASSERT_EQUAL(monitor.converged(), true);
223 }
224
225 // test sparsity strategy
226 {
227 cusp::array1d<ValueType,MemorySpace> x_solve = x;
228 cusp::precond::bridson_ainv<ValueType,MemorySpace> M(A, 0, 10);
229
230 cusp::default_monitor<ValueType> monitor(b, 70, 0, 1e-5);
231 cusp::krylov::cg(A, x_solve, b, monitor, M);
232
233 ASSERT_EQUAL(monitor.converged(), true);
234 }
235
236 // test both
237 {
238 cusp::array1d<ValueType,MemorySpace> x_solve = x;
239 cusp::precond::bridson_ainv<ValueType,MemorySpace> M(A, .01, 4);
240
241 cusp::default_monitor<ValueType> monitor(b, 120, 0, 1e-5);
242 cusp::krylov::cg(A, x_solve, b, monitor, M);
243
244 ASSERT_EQUAL(monitor.converged(), true);
245 }
246
247 // test lin dropping
248 {
249 cusp::array1d<ValueType,MemorySpace> x_solve = x;
250 cusp::precond::bridson_ainv<ValueType,MemorySpace> M(A, 0, -1, true, 4);
251
252 cusp::default_monitor<ValueType> monitor(b, 120, 0, 1e-5);
253 cusp::krylov::cg(A, x_solve, b, monitor, M);
254
255 ASSERT_EQUAL(monitor.converged(), true);
256 }
257
258 }
259 DECLARE_UNITTEST(TestAINVConvergence);
260

  ViewVC Help
Powered by ViewVC 1.1.26