source: git/Singular/svd/blas.h @ aad4ca4

spielwiese
Last change on this file since aad4ca4 was aad4ca4, checked in by Hans Schoenemann <hannes@…>, 6 years ago
format: Warn -> WarnS, trailing spaces
  • Property mode set to 100644
File size: 18.1 KB
Line 
1/*************************************************************************
2Copyright (c) 2005-2007, Sergey Bochkanov (ALGLIB project).
3
4Redistribution and use in source and binary forms, with or without
5modification, are permitted provided that the following conditions are
6met:
7
8- Redistributions of source code must retain the above copyright
9  notice, this list of conditions and the following disclaimer.
10
11- Redistributions in binary form must reproduce the above copyright
12  notice, this list of conditions and the following disclaimer listed
13  in this license in the documentation and/or other materials
14  provided with the distribution.
15
16- Neither the name of the copyright holders nor the names of its
17  contributors may be used to endorse or promote products derived from
18  this software without specific prior written permission.
19
20THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
23A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
24OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
25SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
26LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
27DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
28THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
29(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
30OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31*************************************************************************/
32
33#ifndef _blas_h
34#define _blas_h
35
36#include "ap.h"
37#include "amp.h"
38namespace blas
39{
40    template<unsigned int Precision>
41    amp::ampf<Precision> vectornorm2(const ap::template_1d_array< amp::ampf<Precision> >& x,
42        int i1,
43        int i2);
44    template<unsigned int Precision>
45    int vectoridxabsmax(const ap::template_1d_array< amp::ampf<Precision> >& x,
46        int i1,
47        int i2);
48    template<unsigned int Precision>
49    int columnidxabsmax(const ap::template_2d_array< amp::ampf<Precision> >& x,
50        int i1,
51        int i2,
52        int j);
53    template<unsigned int Precision>
54    int rowidxabsmax(const ap::template_2d_array< amp::ampf<Precision> >& x,
55        int j1,
56        int j2,
57        int i);
58    template<unsigned int Precision>
59    amp::ampf<Precision> upperhessenberg1norm(const ap::template_2d_array< amp::ampf<Precision> >& a,
60        int i1,
61        int i2,
62        int j1,
63        int j2,
64        ap::template_1d_array< amp::ampf<Precision> >& work);
65    template<unsigned int Precision>
66    void copymatrix(const ap::template_2d_array< amp::ampf<Precision> >& a,
67        int is1,
68        int is2,
69        int js1,
70        int js2,
71        ap::template_2d_array< amp::ampf<Precision> >& b,
72        int id1,
73        int id2,
74        int jd1,
75        int jd2);
76    template<unsigned int Precision>
77    void inplacetranspose(ap::template_2d_array< amp::ampf<Precision> >& a,
78        int i1,
79        int i2,
80        int j1,
81        int j2,
82        ap::template_1d_array< amp::ampf<Precision> >& work);
83    template<unsigned int Precision>
84    void copyandtranspose(const ap::template_2d_array< amp::ampf<Precision> >& a,
85        int is1,
86        int is2,
87        int js1,
88        int js2,
89        ap::template_2d_array< amp::ampf<Precision> >& b,
90        int id1,
91        int id2,
92        int jd1,
93        int jd2);
94    template<unsigned int Precision>
95    void matrixvectormultiply(const ap::template_2d_array< amp::ampf<Precision> >& a,
96        int i1,
97        int i2,
98        int j1,
99        int j2,
100        bool trans,
101        const ap::template_1d_array< amp::ampf<Precision> >& x,
102        int ix1,
103        int ix2,
104        amp::ampf<Precision> alpha,
105        ap::template_1d_array< amp::ampf<Precision> >& y,
106        int iy1,
107        int iy2,
108        amp::ampf<Precision> beta);
109    template<unsigned int Precision>
110    amp::ampf<Precision> pythag2(amp::ampf<Precision> x,
111        amp::ampf<Precision> y);
112    template<unsigned int Precision>
113    void matrixmatrixmultiply(const ap::template_2d_array< amp::ampf<Precision> >& a,
114        int ai1,
115        int ai2,
116        int aj1,
117        int aj2,
118        bool transa,
119        const ap::template_2d_array< amp::ampf<Precision> >& b,
120        int bi1,
121        int bi2,
122        int bj1,
123        int bj2,
124        bool transb,
125        amp::ampf<Precision> alpha,
126        ap::template_2d_array< amp::ampf<Precision> >& c,
127        int ci1,
128        int ci2,
129        int cj1,
130        int cj2,
131        amp::ampf<Precision> beta,
132        ap::template_1d_array< amp::ampf<Precision> >& work);
133
134
135    template<unsigned int Precision>
136    amp::ampf<Precision> vectornorm2(const ap::template_1d_array< amp::ampf<Precision> >& x,
137        int i1,
138        int i2)
139    {
140        amp::ampf<Precision> result;
141        int n;
142        int ix;
143        amp::ampf<Precision> absxi;
144        amp::ampf<Precision> scl;
145        amp::ampf<Precision> ssq;
146
147
148        n = i2-i1+1;
149        if( n<1 )
150        {
151            result = 0;
152            return result;
153        }
154        if( n==1 )
155        {
156            result = amp::abs<Precision>(x(i1));
157            return result;
158        }
159        scl = 0;
160        ssq = 1;
161        for(ix=i1; ix<=i2; ix++)
162        {
163            if( x(ix)!=0 )
164            {
165                absxi = amp::abs<Precision>(x(ix));
166                if( scl<absxi )
167                {
168                    ssq = 1+ssq*amp::sqr<Precision>(scl/absxi);
169                    scl = absxi;
170                }
171                else
172                {
173                    ssq = ssq+amp::sqr<Precision>(absxi/scl);
174                }
175            }
176        }
177        result = scl*amp::sqrt<Precision>(ssq);
178        return result;
179    }
180
181
182    template<unsigned int Precision>
183    int vectoridxabsmax(const ap::template_1d_array< amp::ampf<Precision> >& x,
184        int i1,
185        int i2)
186    {
187        int result;
188        int i;
189        amp::ampf<Precision> a;
190
191
192        result = i1;
193        a = amp::abs<Precision>(x(result));
194        for(i=i1+1; i<=i2; i++)
195        {
196            if( amp::abs<Precision>(x(i))>amp::abs<Precision>(x(result)) )
197            {
198                result = i;
199            }
200        }
201        return result;
202    }
203
204
205    template<unsigned int Precision>
206    int columnidxabsmax(const ap::template_2d_array< amp::ampf<Precision> >& x,
207        int i1,
208        int i2,
209        int j)
210    {
211        int result;
212        int i;
213        amp::ampf<Precision> a;
214
215
216        result = i1;
217        a = amp::abs<Precision>(x(result,j));
218        for(i=i1+1; i<=i2; i++)
219        {
220            if( amp::abs<Precision>(x(i,j))>amp::abs<Precision>(x(result,j)) )
221            {
222                result = i;
223            }
224        }
225        return result;
226    }
227
228
229    template<unsigned int Precision>
230    int rowidxabsmax(const ap::template_2d_array< amp::ampf<Precision> >& x,
231        int j1,
232        int j2,
233        int i)
234    {
235        int result;
236        int j;
237        amp::ampf<Precision> a;
238
239
240        result = j1;
241        a = amp::abs<Precision>(x(i,result));
242        for(j=j1+1; j<=j2; j++)
243        {
244            if( amp::abs<Precision>(x(i,j))>amp::abs<Precision>(x(i,result)) )
245            {
246                result = j;
247            }
248        }
249        return result;
250    }
251
252
253    template<unsigned int Precision>
254    amp::ampf<Precision> upperhessenberg1norm(const ap::template_2d_array< amp::ampf<Precision> >& a,
255        int i1,
256        int i2,
257        int j1,
258        int j2,
259        ap::template_1d_array< amp::ampf<Precision> >& work)
260    {
261        amp::ampf<Precision> result;
262        int i;
263        int j;
264
265
266        ap::ap_error::make_assertion(i2-i1==j2-j1);
267        for(j=j1; j<=j2; j++)
268        {
269            work(j) = 0;
270        }
271        for(i=i1; i<=i2; i++)
272        {
273            for(j=ap::maxint(j1, j1+i-i1-1); j<=j2; j++)
274            {
275                work(j) = work(j)+amp::abs<Precision>(a(i,j));
276            }
277        }
278        result = 0;
279        for(j=j1; j<=j2; j++)
280        {
281            result = amp::maximum<Precision>(result, work(j));
282        }
283        return result;
284    }
285
286
287    template<unsigned int Precision>
288    void copymatrix(const ap::template_2d_array< amp::ampf<Precision> >& a,
289        int is1,
290        int is2,
291        int js1,
292        int js2,
293        ap::template_2d_array< amp::ampf<Precision> >& b,
294        int id1,
295        int id2,
296        int jd1,
297        int jd2)
298    {
299        int isrc;
300        int idst;
301
302
303        if( is1>is2 || js1>js2 )
304        {
305            return;
306        }
307        ap::ap_error::make_assertion(is2-is1==id2-id1);
308        ap::ap_error::make_assertion(js2-js1==jd2-jd1);
309        for(isrc=is1; isrc<=is2; isrc++)
310        {
311            idst = isrc-is1+id1;
312            ap::vmove(b.getrow(idst, jd1, jd2), a.getrow(isrc, js1, js2));
313        }
314    }
315
316
317    template<unsigned int Precision>
318    void inplacetranspose(ap::template_2d_array< amp::ampf<Precision> >& a,
319        int i1,
320        int i2,
321        int j1,
322        int j2,
323        ap::template_1d_array< amp::ampf<Precision> >& work)
324    {
325        int i;
326        int j;
327        int ips;
328        int jps;
329        int l;
330
331
332        if( i1>i2 || j1>j2 )
333        {
334            return;
335        }
336        ap::ap_error::make_assertion(i1-i2==j1-j2);
337        for(i=i1; i<=i2-1; i++)
338        {
339            j = j1+i-i1;
340            ips = i+1;
341            jps = j1+ips-i1;
342            l = i2-i;
343            ap::vmove(work.getvector(1, l), a.getcolumn(j, ips, i2));
344            ap::vmove(a.getcolumn(j, ips, i2), a.getrow(i, jps, j2));
345            ap::vmove(a.getrow(i, jps, j2), work.getvector(1, l));
346        }
347    }
348
349
350    template<unsigned int Precision>
351    void copyandtranspose(const ap::template_2d_array< amp::ampf<Precision> >& a,
352        int is1,
353        int is2,
354        int js1,
355        int js2,
356        ap::template_2d_array< amp::ampf<Precision> >& b,
357        int id1,
358        int id2,
359        int jd1,
360        int jd2)
361    {
362        int isrc;
363        int jdst;
364
365
366        if( is1>is2 || js1>js2 )
367        {
368            return;
369        }
370        ap::ap_error::make_assertion(is2-is1==jd2-jd1);
371        ap::ap_error::make_assertion(js2-js1==id2-id1);
372        for(isrc=is1; isrc<=is2; isrc++)
373        {
374            jdst = isrc-is1+jd1;
375            ap::vmove(b.getcolumn(jdst, id1, id2), a.getrow(isrc, js1, js2));
376        }
377    }
378
379
380    template<unsigned int Precision>
381    void matrixvectormultiply(const ap::template_2d_array< amp::ampf<Precision> >& a,
382        int i1,
383        int i2,
384        int j1,
385        int j2,
386        bool trans,
387        const ap::template_1d_array< amp::ampf<Precision> >& x,
388        int ix1,
389        int ix2,
390        amp::ampf<Precision> alpha,
391        ap::template_1d_array< amp::ampf<Precision> >& y,
392        int iy1,
393        int iy2,
394        amp::ampf<Precision> beta)
395    {
396        int i;
397        amp::ampf<Precision> v;
398
399
400        if( !trans )
401        {
402
403            //
404            // y := alpha*A*x + beta*y;
405            //
406            if( i1>i2 || j1>j2 )
407            {
408                return;
409            }
410            ap::ap_error::make_assertion(j2-j1==ix2-ix1);
411            ap::ap_error::make_assertion(i2-i1==iy2-iy1);
412
413            //
414            // beta*y
415            //
416            if( beta==0 )
417            {
418                for(i=iy1; i<=iy2; i++)
419                {
420                    y(i) = 0;
421                }
422            }
423            else
424            {
425                ap::vmul(y.getvector(iy1, iy2), beta);
426            }
427
428            //
429            // alpha*A*x
430            //
431            for(i=i1; i<=i2; i++)
432            {
433                v = ap::vdotproduct(a.getrow(i, j1, j2), x.getvector(ix1, ix2));
434                y(iy1+i-i1) = y(iy1+i-i1)+alpha*v;
435            }
436        }
437        else
438        {
439
440            //
441            // y := alpha*A'*x + beta*y;
442            //
443            if( i1>i2 || j1>j2 )
444            {
445                return;
446            }
447            ap::ap_error::make_assertion(i2-i1==ix2-ix1);
448            ap::ap_error::make_assertion(j2-j1==iy2-iy1);
449
450            //
451            // beta*y
452            //
453            if( beta==0 )
454            {
455                for(i=iy1; i<=iy2; i++)
456                {
457                    y(i) = 0;
458                }
459            }
460            else
461            {
462                ap::vmul(y.getvector(iy1, iy2), beta);
463            }
464
465            //
466            // alpha*A'*x
467            //
468            for(i=i1; i<=i2; i++)
469            {
470                v = alpha*x(ix1+i-i1);
471                ap::vadd(y.getvector(iy1, iy2), a.getrow(i, j1, j2), v);
472            }
473        }
474    }
475
476
477    template<unsigned int Precision>
478    amp::ampf<Precision> pythag2(amp::ampf<Precision> x,
479        amp::ampf<Precision> y)
480    {
481        amp::ampf<Precision> result;
482        amp::ampf<Precision> w;
483        amp::ampf<Precision> xabs;
484        amp::ampf<Precision> yabs;
485        amp::ampf<Precision> z;
486
487
488        xabs = amp::abs<Precision>(x);
489        yabs = amp::abs<Precision>(y);
490        w = amp::maximum<Precision>(xabs, yabs);
491        z = amp::minimum<Precision>(xabs, yabs);
492        if( z==0 )
493        {
494            result = w;
495        }
496        else
497        {
498            result = w*amp::sqrt<Precision>(1+amp::sqr<Precision>(z/w));
499        }
500        return result;
501    }
502
503
504    template<unsigned int Precision>
505    void matrixmatrixmultiply(const ap::template_2d_array< amp::ampf<Precision> >& a,
506        int ai1,
507        int ai2,
508        int aj1,
509        int aj2,
510        bool transa,
511        const ap::template_2d_array< amp::ampf<Precision> >& b,
512        int bi1,
513        int bi2,
514        int bj1,
515        int bj2,
516        bool transb,
517        amp::ampf<Precision> alpha,
518        ap::template_2d_array< amp::ampf<Precision> >& c,
519        int ci1,
520        int ci2,
521        int cj1,
522        int cj2,
523        amp::ampf<Precision> beta,
524        ap::template_1d_array< amp::ampf<Precision> >& work)
525    {
526        int arows;
527        int acols;
528        int brows;
529        int bcols;
530        int crows;
531        int ccols;
532        int i;
533        int j;
534        int k;
535        int l;
536        int r;
537        amp::ampf<Precision> v;
538
539
540
541        //
542        // Setup
543        //
544        if( !transa )
545        {
546            arows = ai2-ai1+1;
547            acols = aj2-aj1+1;
548        }
549        else
550        {
551            arows = aj2-aj1+1;
552            acols = ai2-ai1+1;
553        }
554        if( !transb )
555        {
556            brows = bi2-bi1+1;
557            bcols = bj2-bj1+1;
558        }
559        else
560        {
561            brows = bj2-bj1+1;
562            bcols = bi2-bi1+1;
563        }
564        ap::ap_error::make_assertion(acols==brows);
565        if( arows<=0 || acols<=0 || brows<=0 || bcols<=0 )
566        {
567            return;
568        }
569        crows = arows;
570        ccols = bcols;
571
572        //
573        // Test WORK
574        //
575        i = ap::maxint(arows, acols);
576        i = ap::maxint(brows, i);
577        i = ap::maxint(i, bcols);
578        work(1) = 0;
579        work(i) = 0;
580
581        //
582        // Prepare C
583        //
584        if( beta==0 )
585        {
586            for(i=ci1; i<=ci2; i++)
587            {
588                for(j=cj1; j<=cj2; j++)
589                {
590                    c(i,j) = 0;
591                }
592            }
593        }
594        else
595        {
596            for(i=ci1; i<=ci2; i++)
597            {
598                ap::vmul(c.getrow(i, cj1, cj2), beta);
599            }
600        }
601
602        //
603        // A*B
604        //
605        if( !transa && !transb )
606        {
607            for(l=ai1; l<=ai2; l++)
608            {
609                for(r=bi1; r<=bi2; r++)
610                {
611                    v = alpha*a(l,aj1+r-bi1);
612                    k = ci1+l-ai1;
613                    ap::vadd(c.getrow(k, cj1, cj2), b.getrow(r, bj1, bj2), v);
614                }
615            }
616            return;
617        }
618
619        //
620        // A*B'
621        //
622        if( !transa && transb )
623        {
624            if( arows*acols<brows*bcols )
625            {
626                for(r=bi1; r<=bi2; r++)
627                {
628                    for(l=ai1; l<=ai2; l++)
629                    {
630                        v = ap::vdotproduct(a.getrow(l, aj1, aj2), b.getrow(r, bj1, bj2));
631                        c(ci1+l-ai1,cj1+r-bi1) = c(ci1+l-ai1,cj1+r-bi1)+alpha*v;
632                    }
633                }
634                return;
635            }
636            else
637            {
638                for(l=ai1; l<=ai2; l++)
639                {
640                    for(r=bi1; r<=bi2; r++)
641                    {
642                        v = ap::vdotproduct(a.getrow(l, aj1, aj2), b.getrow(r, bj1, bj2));
643                        c(ci1+l-ai1,cj1+r-bi1) = c(ci1+l-ai1,cj1+r-bi1)+alpha*v;
644                    }
645                }
646                return;
647            }
648        }
649
650        //
651        // A'*B
652        //
653        if( transa && !transb )
654        {
655            for(l=aj1; l<=aj2; l++)
656            {
657                for(r=bi1; r<=bi2; r++)
658                {
659                    v = alpha*a(ai1+r-bi1,l);
660                    k = ci1+l-aj1;
661                    ap::vadd(c.getrow(k, cj1, cj2), b.getrow(r, bj1, bj2), v);
662                }
663            }
664            return;
665        }
666
667        //
668        // A'*B'
669        //
670        if( transa && transb )
671        {
672            if( arows*acols<brows*bcols )
673            {
674                for(r=bi1; r<=bi2; r++)
675                {
676                    for(i=1; i<=crows; i++)
677                    {
678                        work(i) = amp::ampf<Precision>("0.0");
679                    }
680                    for(l=ai1; l<=ai2; l++)
681                    {
682                        v = alpha*b(r,bj1+l-ai1);
683                        k = cj1+r-bi1;
684                        ap::vadd(work.getvector(1, crows), a.getrow(l, aj1, aj2), v);
685                    }
686                    ap::vadd(c.getcolumn(k, ci1, ci2), work.getvector(1, crows));
687                }
688                return;
689            }
690            else
691            {
692                for(l=aj1; l<=aj2; l++)
693                {
694                    k = ai2-ai1+1;
695                    ap::vmove(work.getvector(1, k), a.getcolumn(l, ai1, ai2));
696                    for(r=bi1; r<=bi2; r++)
697                    {
698                        v = ap::vdotproduct(work.getvector(1, k), b.getrow(r, bj1, bj2));
699                        c(ci1+l-aj1,cj1+r-bi1) = c(ci1+l-aj1,cj1+r-bi1)+alpha*v;
700                    }
701                }
702                return;
703            }
704        }
705    }
706} // namespace
707
708#endif
Note: See TracBrowser for help on using the repository browser.