1 /* =========================================================================
2 Copyright (c) 2010-2016, Institute for Microelectronics,
3 Institute for Analysis and Scientific Computing,
4 TU Wien.
5 Portions of this software are copyright by UChicago Argonne, LLC.
6
7 -----------------
8 ViennaCL - The Vienna Computing Library
9 -----------------
10
11 Project Head: Karl Rupp rupp@iue.tuwien.ac.at
12
13 (A list of authors and contributors can be found in the PDF manual)
14
15 License: MIT (X11), see file LICENSE in the base directory
16 ============================================================================= */
17
18
19
20 /** \file tests/src/svd.cpp Tests the singular value decomposition.
21 * \test Tests the singular value decomposition.
22 **/
23
24 #include <stdexcept>
25 #include <iostream>
26 #include <string>
27 #include <vector>
28 #include <cmath>
29
30 #include "viennacl/matrix.hpp"
31 #include "viennacl/linalg/prod.hpp"
32
33 #include "viennacl/linalg/svd.hpp"
34
35 #include "viennacl/tools/timer.hpp"
36
37
read_matrix_size(std::fstream & f,std::size_t & sz1,std::size_t & sz2)38 inline void read_matrix_size(std::fstream& f, std::size_t & sz1, std::size_t & sz2)
39 {
40 if (!f.is_open())
41 throw std::invalid_argument("File is not opened");
42
43 f >> sz1 >> sz2;
44 }
45
46
47 template<typename ScalarType>
read_matrix_body(std::fstream & f,viennacl::matrix<ScalarType> & A)48 void read_matrix_body(std::fstream& f, viennacl::matrix<ScalarType>& A)
49 {
50 if (!f.is_open())
51 throw std::invalid_argument("File is not opened");
52
53 boost::numeric::ublas::matrix<ScalarType> h_A(A.size1(), A.size2());
54
55 for (std::size_t i = 0; i < h_A.size1(); i++)
56 {
57 for (std::size_t j = 0; j < h_A.size2(); j++)
58 {
59 ScalarType val = 0.0;
60 f >> val;
61 h_A(i, j) = val;
62 }
63 }
64
65 viennacl::copy(h_A, A);
66 }
67
68
69 template<typename ScalarType>
read_vector_body(std::fstream & f,std::vector<ScalarType> & v)70 void read_vector_body(std::fstream& f, std::vector<ScalarType>& v)
71 {
72 if (!f.is_open())
73 throw std::invalid_argument("File is not opened");
74
75 for (std::size_t i = 0; i < v.size(); i++)
76 {
77 ScalarType val = 0.0;
78 f >> val;
79 v[i] = val;
80 }
81 }
82
83
84 template<typename ScalarType>
random_fill(std::vector<ScalarType> & in)85 void random_fill(std::vector<ScalarType>& in)
86 {
87 for (std::size_t i = 0; i < in.size(); i++)
88 in[i] = static_cast<ScalarType>(rand()) / ScalarType(RAND_MAX);
89 }
90
91
92 template<typename ScalarType>
check_bidiag(viennacl::matrix<ScalarType> & A)93 bool check_bidiag(viennacl::matrix<ScalarType>& A)
94 {
95 const ScalarType EPS = 0.0001f;
96
97 std::vector<ScalarType> aA(A.size1() * A.size2());
98 viennacl::fast_copy(A, &aA[0]);
99
100 for (std::size_t i = 0; i < A.size1(); i++)
101 {
102 for (std::size_t j = 0; j < A.size2(); j++)
103 {
104 ScalarType val = aA[i * A.size2() + j];
105 if ((fabs(val) > EPS) && (i != j) && ((i + 1) != j))
106 {
107 std::cout << "Failed at " << i << " " << j << " " << val << std::endl;
108 return false;
109 }
110 }
111 }
112
113 return true;
114 }
115
116 template<typename ScalarType>
matrix_compare(viennacl::matrix<ScalarType> & res,viennacl::matrix<ScalarType> & ref)117 ScalarType matrix_compare(viennacl::matrix<ScalarType>& res,
118 viennacl::matrix<ScalarType>& ref)
119 {
120 std::vector<ScalarType> res_std(res.internal_size());
121 std::vector<ScalarType> ref_std(ref.internal_size());
122
123 viennacl::fast_copy(res, &res_std[0]);
124 viennacl::fast_copy(ref, &ref_std[0]);
125
126 ScalarType diff = 0.0;
127 ScalarType mx = 0.0;
128
129 for (std::size_t i = 0; i < res_std.size(); i++)
130 {
131 diff = std::max(diff, std::abs(res_std[i] - ref_std[i]));
132 mx = std::max(mx, res_std[i]);
133 }
134
135 return diff / mx;
136 }
137
138
139 template<typename ScalarType>
sigmas_compare(viennacl::matrix<ScalarType> & res,std::vector<ScalarType> & ref)140 ScalarType sigmas_compare(viennacl::matrix<ScalarType>& res,
141 std::vector<ScalarType>& ref)
142 {
143 std::vector<ScalarType> res_std(ref.size());
144
145 for (std::size_t i = 0; i < ref.size(); i++)
146 res_std[i] = res(i, i);
147
148 std::sort(ref.begin(), ref.end());
149 std::sort(res_std.begin(), res_std.end());
150
151 ScalarType diff = 0.0;
152 ScalarType mx = 0.0;
153 for (std::size_t i = 0; i < ref.size(); i++)
154 {
155 diff = std::max(diff, std::abs(res_std[i] - ref[i]));
156 mx = std::max(mx, res_std[i]);
157 }
158
159 return diff / mx;
160 }
161
162
163 template<typename ScalarType>
test_svd(const std::string & fn,ScalarType EPS)164 void test_svd(const std::string & fn, ScalarType EPS)
165 {
166 std::size_t sz1, sz2;
167
168 //read matrix
169
170 // sz1 = 2048, sz2 = 2048;
171 // std::vector<ScalarType> in(sz1 * sz2);
172 // random_fill(in);
173
174 // read file
175 std::fstream f(fn.c_str(), std::fstream::in);
176 //read size of input matrix
177 read_matrix_size(f, sz1, sz2);
178
179 std::size_t to = std::min(sz1, sz2);
180
181 viennacl::matrix<ScalarType> Ai(sz1, sz2), Aref(sz1, sz2), QL(sz1, sz1), QR(sz2, sz2);
182 read_matrix_body(f, Ai);
183
184 std::vector<ScalarType> sigma_ref(to);
185 read_vector_body(f, sigma_ref);
186
187 f.close();
188
189 // viennacl::fast_copy(&in[0], &in[0] + in.size(), Ai);
190
191 Aref = Ai;
192
193 viennacl::tools::timer timer;
194 timer.start();
195
196 viennacl::linalg::svd(Ai, QL, QR);
197
198 viennacl::backend::finish();
199
200 double time_spend = timer.get();
201
202 viennacl::matrix<ScalarType> result1(sz1, sz2), result2(sz1, sz2);
203 result1 = viennacl::linalg::prod(QL, Ai);
204 result2 = viennacl::linalg::prod(result1, trans(QR));
205
206 ScalarType sigma_diff = sigmas_compare(Ai, sigma_ref);
207 ScalarType prods_diff = matrix_compare(result2, Aref);
208
209 bool sigma_ok = (fabs(sigma_diff) < EPS)
210 && (fabs(prods_diff) < std::sqrt(EPS)); //note: computing the product is not accurate down to 10^{-16}, so we allow for a higher tolerance here
211
212 printf("%6s [%dx%d] %40s sigma_diff = %.6f; prod_diff = %.6f; time = %.6f\n", sigma_ok?"[[OK]]":"[FAIL]", (int)Aref.size1(), (int)Aref.size2(), fn.c_str(), sigma_diff, prods_diff, time_spend);
213 if (!sigma_ok)
214 exit(EXIT_FAILURE);
215 }
216
217
218 template<typename ScalarType>
time_svd(std::size_t sz1,std::size_t sz2)219 void time_svd(std::size_t sz1, std::size_t sz2)
220 {
221 viennacl::matrix<ScalarType> Ai(sz1, sz2), QL(sz1, sz1), QR(sz2, sz2);
222
223 std::vector<ScalarType> in(Ai.internal_size1() * Ai.internal_size2());
224 random_fill(in);
225
226 viennacl::fast_copy(&in[0], &in[0] + in.size(), Ai);
227
228
229 viennacl::tools::timer timer;
230 timer.start();
231
232 viennacl::linalg::svd(Ai, QL, QR);
233 viennacl::backend::finish();
234 double time_spend = timer.get();
235
236 printf("[%dx%d] time = %.6f\n", static_cast<int>(sz1), static_cast<int>(sz2), time_spend);
237 }
238
239
240 template<typename ScalarType>
test(ScalarType epsilon)241 int test(ScalarType epsilon)
242 {
243
244 test_svd<ScalarType>(std::string("../examples/testdata/svd/qr.example"), epsilon);
245 test_svd<ScalarType>(std::string("../examples/testdata/svd/wiki.example"), epsilon);
246 test_svd<ScalarType>(std::string("../examples/testdata/svd/wiki.qr.example"), epsilon);
247 test_svd<ScalarType>(std::string("../examples/testdata/svd/pysvd.example"), epsilon);
248 test_svd<ScalarType>(std::string("../examples/testdata/svd/random.example"), epsilon);
249
250 time_svd<ScalarType>(500, 500);
251 time_svd<ScalarType>(1024, 1024);
252 time_svd<ScalarType>(2048, 512);
253 //time_svd<ScalarType>(2048, 2048);
254 //time_svd(4096, 4096); //takes too long for a standard sanity test. Feel free to uncomment
255
256 return EXIT_SUCCESS;
257 }
258
259 //
260 // -------------------------------------------------------------
261 //
main()262 int main()
263 {
264 std::cout << std::endl;
265 std::cout << "----------------------------------------------" << std::endl;
266 std::cout << "----------------------------------------------" << std::endl;
267 std::cout << "## Test :: Singular Value Decomposition" << std::endl;
268 std::cout << "----------------------------------------------" << std::endl;
269 std::cout << "----------------------------------------------" << std::endl;
270 std::cout << std::endl;
271
272 int retval = EXIT_SUCCESS;
273
274 std::cout << std::endl;
275 std::cout << "----------------------------------------------" << std::endl;
276 std::cout << std::endl;
277 {
278 typedef float NumericT;
279 NumericT epsilon = NumericT(1.0E-4);
280 std::cout << "# Testing setup:" << std::endl;
281 std::cout << " eps: " << epsilon << std::endl;
282 std::cout << " numeric: float" << std::endl;
283 retval = test<NumericT>(epsilon);
284 if ( retval == EXIT_SUCCESS )
285 std::cout << "# Test passed" << std::endl;
286 else
287 return retval;
288 }
289 std::cout << std::endl;
290 std::cout << "----------------------------------------------" << std::endl;
291 std::cout << std::endl;
292 if ( viennacl::ocl::current_device().double_support() )
293 {
294 {
295 typedef double NumericT;
296 NumericT epsilon = 1.0E-6; //Note: higher accuracy not possible, because data only available with floating point precision
297 std::cout << "# Testing setup:" << std::endl;
298 std::cout << " eps: " << epsilon << std::endl;
299 std::cout << " numeric: double" << std::endl;
300 retval = test<NumericT>(epsilon);
301 if ( retval == EXIT_SUCCESS )
302 std::cout << "# Test passed" << std::endl;
303 else
304 return retval;
305 }
306 std::cout << std::endl;
307 std::cout << "----------------------------------------------" << std::endl;
308 std::cout << std::endl;
309 }
310
311 std::cout << std::endl;
312 std::cout << "------- Test completed --------" << std::endl;
313 std::cout << std::endl;
314
315
316 return retval;
317 }
318
319
320