source: git/Singular/svd/tests/testsvdunit.h @ 7fb70f

spielwiese
Last change on this file since 7fb70f was 7fb70f, checked in by Hans Schoenemann <hannes@…>, 7 years ago
add (from 3-1): system("svd",m)
  • Property mode set to 100644
File size: 14.6 KB
Line 
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 int failcount;
51    static 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
Note: See TracBrowser for help on using the repository browser.