My Project
Loading...
Searching...
No Matches
testsvdunit.h
Go to the documentation of this file.
1
2#ifndef _testsvdunit_h
3#define _testsvdunit_h
4
5#include <stdio.h>
6#include "ap.h"
7#include "amp.h"
8#include "reflections.h"
9#include "bidiagonal.h"
10#include "qr.h"
11#include "lq.h"
12#include "blas.h"
13#include "rotations.h"
14#include "bdsvd.h"
15#include "svd.h"
16namespace testsvdunit
17{
18 template<unsigned int Precision>
19 bool testsvd(bool silent);
20 template<unsigned int Precision>
21 void fillsparsea(ap::template_2d_array< amp::ampf<Precision> >& a,
22 int m,
23 int n,
24 amp::ampf<Precision> sparcity);
25 template<unsigned int Precision>
26 void getsvderror(const ap::template_2d_array< amp::ampf<Precision> >& a,
27 int m,
28 int n,
29 const ap::template_2d_array< amp::ampf<Precision> >& u,
30 const ap::template_1d_array< amp::ampf<Precision> >& w,
31 const ap::template_2d_array< amp::ampf<Precision> >& vt,
32 amp::ampf<Precision>& materr,
33 amp::ampf<Precision>& orterr,
34 bool& wsorted);
35 template<unsigned int Precision>
36 void testsvdproblem(const ap::template_2d_array< amp::ampf<Precision> >& a,
37 int m,
38 int n,
39 amp::ampf<Precision>& materr,
40 amp::ampf<Precision>& orterr,
41 amp::ampf<Precision>& othererr,
42 bool& wsorted,
43 bool& wfailed);
44 template<unsigned int Precision>
45 bool testsvdunit_test_silent();
46 template<unsigned int Precision>
47 bool testsvdunit_test();
48
49
50 STATIC_VAR int failcount;
51 STATIC_VAR int succcount;
52
53
54 /*************************************************************************
55 Testing SVD decomposition subroutine
56 *************************************************************************/
57 template<unsigned int Precision>
58 bool testsvd(bool silent)
59 {
60 bool result;
61 ap::template_2d_array< amp::ampf<Precision> > a;
62 int m;
63 int n;
64 int maxmn;
65 int i;
66 int j;
67 int gpass;
68 int pass;
69 bool waserrors;
70 bool wsorted;
71 bool wfailed;
72 amp::ampf<Precision> materr;
73 amp::ampf<Precision> orterr;
74 amp::ampf<Precision> othererr;
75 amp::ampf<Precision> threshold;
76 amp::ampf<Precision> failthreshold;
77 amp::ampf<Precision> failr;
78
79
80 failcount = 0;
81 succcount = 0;
82 materr = 0;
83 orterr = 0;
84 othererr = 0;
85 wsorted = true;
86 wfailed = false;
87 waserrors = false;
88 maxmn = 30;
89 threshold = 5*100*amp::ampf<Precision>::getAlgoPascalEpsilon();
90 failthreshold = amp::ampf<Precision>("5.0E-3");
91 a.setbounds(0, maxmn-1, 0, maxmn-1);
92
93 //
94 // TODO: div by zero fail, convergence fail
95 //
96 for(gpass=1; gpass<=1; gpass++)
97 {
98
99 //
100 // zero matrix, several cases
101 //
102 for(i=0; i<=maxmn-1; i++)
103 {
104 for(j=0; j<=maxmn-1; j++)
105 {
106 a(i,j) = 0;
107 }
108 }
109 for(i=1; i<=ap::minint(5, maxmn); i++)
110 {
111 for(j=1; j<=ap::minint(5, maxmn); j++)
112 {
113 testsvdproblem<Precision>(a, i, j, materr, orterr, othererr, wsorted, wfailed);
114 }
115 }
116
117 //
118 // Long dense matrix
119 //
120 for(i=0; i<=maxmn-1; i++)
121 {
122 for(j=0; j<=ap::minint(5, maxmn)-1; j++)
123 {
124 a(i,j) = 2*amp::ampf<Precision>::getRandom()-1;
125 }
126 }
127 for(i=1; i<=maxmn; i++)
128 {
129 for(j=1; j<=ap::minint(5, maxmn); j++)
130 {
131 testsvdproblem<Precision>(a, i, j, materr, orterr, othererr, wsorted, wfailed);
132 }
133 }
134 for(i=0; i<=ap::minint(5, maxmn)-1; i++)
135 {
136 for(j=0; j<=maxmn-1; j++)
137 {
138 a(i,j) = 2*amp::ampf<Precision>::getRandom()-1;
139 }
140 }
141 for(i=1; i<=ap::minint(5, maxmn); i++)
142 {
143 for(j=1; j<=maxmn; j++)
144 {
145 testsvdproblem<Precision>(a, i, j, materr, orterr, othererr, wsorted, wfailed);
146 }
147 }
148
149 //
150 // Dense matrices
151 //
152 for(m=1; m<=ap::minint(10, maxmn); m++)
153 {
154 for(n=1; n<=ap::minint(10, maxmn); n++)
155 {
156 for(i=0; i<=m-1; i++)
157 {
158 for(j=0; j<=n-1; j++)
159 {
160 a(i,j) = 2*amp::ampf<Precision>::getRandom()-1;
161 }
162 }
163 testsvdproblem<Precision>(a, m, n, materr, orterr, othererr, wsorted, wfailed);
164 }
165 }
166
167 //
168 // Sparse matrices, very sparse matrices, incredible sparse matrices
169 //
170 for(m=1; m<=10; m++)
171 {
172 for(n=1; n<=10; n++)
173 {
174 for(pass=1; pass<=2; pass++)
175 {
176 fillsparsea<Precision>(a, m, n, amp::ampf<Precision>("0.8"));
177 testsvdproblem<Precision>(a, m, n, materr, orterr, othererr, wsorted, wfailed);
178 fillsparsea<Precision>(a, m, n, amp::ampf<Precision>("0.9"));
179 testsvdproblem<Precision>(a, m, n, materr, orterr, othererr, wsorted, wfailed);
180 fillsparsea<Precision>(a, m, n, amp::ampf<Precision>("0.95"));
181 testsvdproblem<Precision>(a, m, n, materr, orterr, othererr, wsorted, wfailed);
182 }
183 }
184 }
185 }
186
187 //
188 // report
189 //
190 failr = amp::ampf<Precision>(failcount)/(amp::ampf<Precision>(succcount+failcount));
191 waserrors = materr>threshold || orterr>threshold || othererr>threshold || !wsorted || failr>failthreshold;
192 if( !silent )
193 {
194 printf("TESTING SVD DECOMPOSITION\n");
195 printf("SVD decomposition error: %5.3le\n",
196 double(amp::ampf<Precision>(materr).toDouble()));
197 printf("SVD orthogonality error: %5.3le\n",
198 double(amp::ampf<Precision>(orterr).toDouble()));
199 printf("SVD with different parameters error: %5.3le\n",
200 double(amp::ampf<Precision>(othererr).toDouble()));
201 printf("Singular values order: ");
202 if( wsorted )
203 {
204 printf("OK\n");
205 }
206 else
207 {
208 printf("FAILED\n");
209 }
210 printf("Always converged: ");
211 if( !wfailed )
212 {
213 printf("YES\n");
214 }
215 else
216 {
217 printf("NO\n");
218 printf("Fail ratio: %5.3lf\n",
219 double(amp::ampf<Precision>(failr).toDouble()));
220 }
221 printf("Threshold: %5.3le\n",
222 double(amp::ampf<Precision>(threshold).toDouble()));
223 if( waserrors )
224 {
225 printf("TEST FAILED\n");
226 }
227 else
228 {
229 printf("TEST PASSED\n");
230 }
231 printf("\n\n");
232 }
233 result = !waserrors;
234 return result;
235 }
236
237
238 template<unsigned int Precision>
239 void fillsparsea(ap::template_2d_array< amp::ampf<Precision> >& a,
240 int m,
241 int n,
242 amp::ampf<Precision> sparcity)
243 {
244 int i;
245 int j;
246
247
248 for(i=0; i<=m-1; i++)
249 {
250 for(j=0; j<=n-1; j++)
251 {
252 if( amp::ampf<Precision>::getRandom()>=sparcity )
253 {
254 a(i,j) = 2*amp::ampf<Precision>::getRandom()-1;
255 }
256 else
257 {
258 a(i,j) = 0;
259 }
260 }
261 }
262 }
263
264
265 template<unsigned int Precision>
266 void getsvderror(const ap::template_2d_array< amp::ampf<Precision> >& a,
267 int m,
268 int n,
269 const ap::template_2d_array< amp::ampf<Precision> >& u,
270 const ap::template_1d_array< amp::ampf<Precision> >& w,
271 const ap::template_2d_array< amp::ampf<Precision> >& vt,
272 amp::ampf<Precision>& materr,
273 amp::ampf<Precision>& orterr,
274 bool& wsorted)
275 {
276 int i;
277 int j;
278 int k;
279 int minmn;
280 amp::ampf<Precision> locerr;
281 amp::ampf<Precision> sm;
282
283
284 minmn = ap::minint(m, n);
285
286 //
287 // decomposition error
288 //
289 locerr = 0;
290 for(i=0; i<=m-1; i++)
291 {
292 for(j=0; j<=n-1; j++)
293 {
294 sm = 0;
295 for(k=0; k<=minmn-1; k++)
296 {
297 sm = sm+w(k)*u(i,k)*vt(k,j);
298 }
299 locerr = amp::maximum<Precision>(locerr, amp::abs<Precision>(a(i,j)-sm));
300 }
301 }
302 materr = amp::maximum<Precision>(materr, locerr);
303
304 //
305 // orthogonality error
306 //
307 locerr = 0;
308 for(i=0; i<=minmn-1; i++)
309 {
310 for(j=i; j<=minmn-1; j++)
311 {
312 sm = ap::vdotproduct(u.getcolumn(i, 0, m-1), u.getcolumn(j, 0, m-1));
313 if( i!=j )
314 {
315 locerr = amp::maximum<Precision>(locerr, amp::abs<Precision>(sm));
316 }
317 else
318 {
319 locerr = amp::maximum<Precision>(locerr, amp::abs<Precision>(sm-1));
320 }
321 sm = ap::vdotproduct(vt.getrow(i, 0, n-1), vt.getrow(j, 0, n-1));
322 if( i!=j )
323 {
324 locerr = amp::maximum<Precision>(locerr, amp::abs<Precision>(sm));
325 }
326 else
327 {
328 locerr = amp::maximum<Precision>(locerr, amp::abs<Precision>(sm-1));
329 }
330 }
331 }
332 orterr = amp::maximum<Precision>(orterr, locerr);
333
334 //
335 // values order error
336 //
337 for(i=1; i<=minmn-1; i++)
338 {
339 if( w(i)>w(i-1) )
340 {
341 wsorted = false;
342 }
343 }
344 }
345
346
347 template<unsigned int Precision>
348 void testsvdproblem(const ap::template_2d_array< amp::ampf<Precision> >& a,
349 int m,
350 int n,
351 amp::ampf<Precision>& materr,
352 amp::ampf<Precision>& orterr,
353 amp::ampf<Precision>& othererr,
354 bool& wsorted,
355 bool& wfailed)
356 {
357 ap::template_2d_array< amp::ampf<Precision> > u;
358 ap::template_2d_array< amp::ampf<Precision> > vt;
359 ap::template_2d_array< amp::ampf<Precision> > u2;
360 ap::template_2d_array< amp::ampf<Precision> > vt2;
361 ap::template_1d_array< amp::ampf<Precision> > w;
362 ap::template_1d_array< amp::ampf<Precision> > w2;
363 int i;
364 int j;
365 int k;
366 int ujob;
367 int vtjob;
368 int memjob;
369 int ucheck;
370 int vtcheck;
371 amp::ampf<Precision> v;
372 amp::ampf<Precision> mx;
373
374
375
376 //
377 // Main SVD test
378 //
379 if( !svd::rmatrixsvd<Precision>(a, m, n, 2, 2, 2, w, u, vt) )
380 {
381 failcount = failcount+1;
382 wfailed = true;
383 return;
384 }
385 getsvderror<Precision>(a, m, n, u, w, vt, materr, orterr, wsorted);
386
387 //
388 // Additional SVD tests
389 //
390 for(ujob=0; ujob<=2; ujob++)
391 {
392 for(vtjob=0; vtjob<=2; vtjob++)
393 {
394 for(memjob=0; memjob<=2; memjob++)
395 {
396 if( !svd::rmatrixsvd<Precision>(a, m, n, ujob, vtjob, memjob, w2, u2, vt2) )
397 {
398 failcount = failcount+1;
399 wfailed = true;
400 return;
401 }
402 ucheck = 0;
403 if( ujob==1 )
404 {
405 ucheck = ap::minint(m, n);
406 }
407 if( ujob==2 )
408 {
409 ucheck = m;
410 }
411 vtcheck = 0;
412 if( vtjob==1 )
413 {
414 vtcheck = ap::minint(m, n);
415 }
416 if( vtjob==2 )
417 {
418 vtcheck = n;
419 }
420 for(i=0; i<=m-1; i++)
421 {
422 for(j=0; j<=ucheck-1; j++)
423 {
424 othererr = amp::maximum<Precision>(othererr, amp::abs<Precision>(u(i,j)-u2(i,j)));
425 }
426 }
427 for(i=0; i<=vtcheck-1; i++)
428 {
429 for(j=0; j<=n-1; j++)
430 {
431 othererr = amp::maximum<Precision>(othererr, amp::abs<Precision>(vt(i,j)-vt2(i,j)));
432 }
433 }
434 for(i=0; i<=ap::minint(m, n)-1; i++)
435 {
436 othererr = amp::maximum<Precision>(othererr, amp::abs<Precision>(w(i)-w2(i)));
437 }
438 }
439 }
440 }
441
442 //
443 // update counter
444 //
445 succcount = succcount+1;
446 }
447
448
449 /*************************************************************************
450 Silent unit test
451 *************************************************************************/
452 template<unsigned int Precision>
453 bool testsvdunit_test_silent()
454 {
455 bool result;
456
457
458 result = testsvd<Precision>(true);
459 return result;
460 }
461
462
463 /*************************************************************************
464 Unit test
465 *************************************************************************/
466 template<unsigned int Precision>
467 bool testsvdunit_test()
468 {
469 bool result;
470
471
472 result = testsvd<Precision>(false);
473 return result;
474 }
475} // namespace
476
477#endif
m
int m
Definition cfEzgcd.cc:128
i
int i
Definition cfEzgcd.cc:132
k
int k
Definition cfEzgcd.cc:99
amp::ampf::getAlgoPascalEpsilon
static const ampf getAlgoPascalEpsilon()
Definition amp.h:565
amp::ampf::getRandom
static const ampf getRandom()
Definition amp.h:593
ap::template_1d_array
Definition ap.h:657
ap::template_2d_array
Definition ap.h:807
ap::template_2d_array::setbounds
void setbounds(int iLow1, int iHigh1, int iLow2, int iHigh2)
Definition ap.h:890
result
return result
Definition facAbsBiFact.cc:76
w
const CanonicalForm & w
Definition facAbsFact.cc:51
v
const Variable & v
< [in] a sqrfree bivariate poly
Definition facBivar.h:39
j
int j
Definition facHensel.cc:110
STATIC_VAR
#define STATIC_VAR
Definition globaldefs.h:7
amp::abs
const ampf< Precision > abs(const ampf< Precision > &x)
Definition amp.h:713
amp::maximum
const ampf< Precision > maximum(const ampf< Precision > &x, const ampf< Precision > &y)
Definition amp.h:722
ap::vdotproduct
T vdotproduct(const_raw_vector< T > v1, const_raw_vector< T > v2)
Definition ap.h:181
ap::minint
int minint(int m1, int m2)
Definition ap.cpp:167
svd::rmatrixsvd
bool rmatrixsvd(ap::template_2d_array< amp::ampf< Precision > > a, int m, int n, int uneeded, int vtneeded, int additionalmemory, ap::template_1d_array< amp::ampf< Precision > > &w, ap::template_2d_array< amp::ampf< Precision > > &u, ap::template_2d_array< amp::ampf< Precision > > &vt)
Definition svd.h:140
testsvdunit::fillsparsea
void fillsparsea(ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, amp::ampf< Precision > sparcity)
Definition testsvdunit.h:239
testsvdunit::testsvdunit_test
bool testsvdunit_test()
Definition testsvdunit.h:467
testsvdunit::testsvd
bool testsvd(bool silent)
Definition testsvdunit.h:58
testsvdunit::getsvderror
void getsvderror(const ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, const ap::template_2d_array< amp::ampf< Precision > > &u, const ap::template_1d_array< amp::ampf< Precision > > &w, const ap::template_2d_array< amp::ampf< Precision > > &vt, amp::ampf< Precision > &materr, amp::ampf< Precision > &orterr, bool &wsorted)
Definition testsvdunit.h:266
testsvdunit::testsvdunit_test_silent
bool testsvdunit_test_silent()
Definition testsvdunit.h:453
testsvdunit::failcount
STATIC_VAR int failcount
Definition testsvdunit.h:50
testsvdunit::testsvdproblem
void testsvdproblem(const ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, amp::ampf< Precision > &materr, amp::ampf< Precision > &orterr, amp::ampf< Precision > &othererr, bool &wsorted, bool &wfailed)
Definition testsvdunit.h:348
testsvdunit::succcount
STATIC_VAR int succcount
Definition testsvdunit.h:51
Generated on for My Project by   doxygen 1.17.0 for Singular