Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.LinearRegression/3.2/alglib/blas.cs @ 2211

Last change on this file since 2211 was 2154, checked in by gkronber, 15 years ago

Added linear regression plugin. #697
File size: 16.9 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
33using System;
34
35class blas
36{
37    public static double vectornorm2(ref double[] x,
38        int i1,
39        int i2)
40    {
41        double result = 0;
42        int n = 0;
43        int ix = 0;
44        double absxi = 0;
45        double scl = 0;
46        double ssq = 0;
47
48        n = i2-i1+1;
49        if( n<1 )
50        {
51            result = 0;
52            return result;
53        }
54        if( n==1 )
55        {
56            result = Math.Abs(x[i1]);
57            return result;
58        }
59        scl = 0;
60        ssq = 1;
61        for(ix=i1; ix<=i2; ix++)
62        {
63            if( x[ix]!=0 )
64            {
65                absxi = Math.Abs(x[ix]);
66                if( scl<absxi )
67                {
68                    ssq = 1+ssq*AP.Math.Sqr(scl/absxi);
69                    scl = absxi;
70                }
71                else
72                {
73                    ssq = ssq+AP.Math.Sqr(absxi/scl);
74                }
75            }
76        }
77        result = scl*Math.Sqrt(ssq);
78        return result;
79    }
80
81
82    public static int vectoridxabsmax(ref double[] x,
83        int i1,
84        int i2)
85    {
86        int result = 0;
87        int i = 0;
88        double a = 0;
89
90        result = i1;
91        a = Math.Abs(x[result]);
92        for(i=i1+1; i<=i2; i++)
93        {
94            if( Math.Abs(x[i])>Math.Abs(x[result]) )
95            {
96                result = i;
97            }
98        }
99        return result;
100    }
101
102
103    public static int columnidxabsmax(ref double[,] x,
104        int i1,
105        int i2,
106        int j)
107    {
108        int result = 0;
109        int i = 0;
110        double a = 0;
111
112        result = i1;
113        a = Math.Abs(x[result,j]);
114        for(i=i1+1; i<=i2; i++)
115        {
116            if( Math.Abs(x[i,j])>Math.Abs(x[result,j]) )
117            {
118                result = i;
119            }
120        }
121        return result;
122    }
123
124
125    public static int rowidxabsmax(ref double[,] x,
126        int j1,
127        int j2,
128        int i)
129    {
130        int result = 0;
131        int j = 0;
132        double a = 0;
133
134        result = j1;
135        a = Math.Abs(x[i,result]);
136        for(j=j1+1; j<=j2; j++)
137        {
138            if( Math.Abs(x[i,j])>Math.Abs(x[i,result]) )
139            {
140                result = j;
141            }
142        }
143        return result;
144    }
145
146
147    public static double upperhessenberg1norm(ref double[,] a,
148        int i1,
149        int i2,
150        int j1,
151        int j2,
152        ref double[] work)
153    {
154        double result = 0;
155        int i = 0;
156        int j = 0;
157
158        System.Diagnostics.Debug.Assert(i2-i1==j2-j1, "UpperHessenberg1Norm: I2-I1<>J2-J1!");
159        for(j=j1; j<=j2; j++)
160        {
161            work[j] = 0;
162        }
163        for(i=i1; i<=i2; i++)
164        {
165            for(j=Math.Max(j1, j1+i-i1-1); j<=j2; j++)
166            {
167                work[j] = work[j]+Math.Abs(a[i,j]);
168            }
169        }
170        result = 0;
171        for(j=j1; j<=j2; j++)
172        {
173            result = Math.Max(result, work[j]);
174        }
175        return result;
176    }
177
178
179    public static void copymatrix(ref double[,] a,
180        int is1,
181        int is2,
182        int js1,
183        int js2,
184        ref double[,] b,
185        int id1,
186        int id2,
187        int jd1,
188        int jd2)
189    {
190        int isrc = 0;
191        int idst = 0;
192        int i_ = 0;
193        int i1_ = 0;
194
195        if( is1>is2 | js1>js2 )
196        {
197            return;
198        }
199        System.Diagnostics.Debug.Assert(is2-is1==id2-id1, "CopyMatrix: different sizes!");
200        System.Diagnostics.Debug.Assert(js2-js1==jd2-jd1, "CopyMatrix: different sizes!");
201        for(isrc=is1; isrc<=is2; isrc++)
202        {
203            idst = isrc-is1+id1;
204            i1_ = (js1) - (jd1);
205            for(i_=jd1; i_<=jd2;i_++)
206            {
207                b[idst,i_] = a[isrc,i_+i1_];
208            }
209        }
210    }
211
212
213    public static void inplacetranspose(ref double[,] a,
214        int i1,
215        int i2,
216        int j1,
217        int j2,
218        ref double[] work)
219    {
220        int i = 0;
221        int j = 0;
222        int ips = 0;
223        int jps = 0;
224        int l = 0;
225        int i_ = 0;
226        int i1_ = 0;
227
228        if( i1>i2 | j1>j2 )
229        {
230            return;
231        }
232        System.Diagnostics.Debug.Assert(i1-i2==j1-j2, "InplaceTranspose error: incorrect array size!");
233        for(i=i1; i<=i2-1; i++)
234        {
235            j = j1+i-i1;
236            ips = i+1;
237            jps = j1+ips-i1;
238            l = i2-i;
239            i1_ = (ips) - (1);
240            for(i_=1; i_<=l;i_++)
241            {
242                work[i_] = a[i_+i1_,j];
243            }
244            i1_ = (jps) - (ips);
245            for(i_=ips; i_<=i2;i_++)
246            {
247                a[i_,j] = a[i,i_+i1_];
248            }
249            i1_ = (1) - (jps);
250            for(i_=jps; i_<=j2;i_++)
251            {
252                a[i,i_] = work[i_+i1_];
253            }
254        }
255    }
256
257
258    public static void copyandtranspose(ref double[,] a,
259        int is1,
260        int is2,
261        int js1,
262        int js2,
263        ref double[,] b,
264        int id1,
265        int id2,
266        int jd1,
267        int jd2)
268    {
269        int isrc = 0;
270        int jdst = 0;
271        int i_ = 0;
272        int i1_ = 0;
273
274        if( is1>is2 | js1>js2 )
275        {
276            return;
277        }
278        System.Diagnostics.Debug.Assert(is2-is1==jd2-jd1, "CopyAndTranspose: different sizes!");
279        System.Diagnostics.Debug.Assert(js2-js1==id2-id1, "CopyAndTranspose: different sizes!");
280        for(isrc=is1; isrc<=is2; isrc++)
281        {
282            jdst = isrc-is1+jd1;
283            i1_ = (js1) - (id1);
284            for(i_=id1; i_<=id2;i_++)
285            {
286                b[i_,jdst] = a[isrc,i_+i1_];
287            }
288        }
289    }
290
291
292    public static void matrixvectormultiply(ref double[,] a,
293        int i1,
294        int i2,
295        int j1,
296        int j2,
297        bool trans,
298        ref double[] x,
299        int ix1,
300        int ix2,
301        double alpha,
302        ref double[] y,
303        int iy1,
304        int iy2,
305        double beta)
306    {
307        int i = 0;
308        double v = 0;
309        int i_ = 0;
310        int i1_ = 0;
311
312        if( !trans )
313        {
314           
315            //
316            // y := alpha*A*x + beta*y;
317            //
318            if( i1>i2 | j1>j2 )
319            {
320                return;
321            }
322            System.Diagnostics.Debug.Assert(j2-j1==ix2-ix1, "MatrixVectorMultiply: A and X dont match!");
323            System.Diagnostics.Debug.Assert(i2-i1==iy2-iy1, "MatrixVectorMultiply: A and Y dont match!");
324           
325            //
326            // beta*y
327            //
328            if( beta==0 )
329            {
330                for(i=iy1; i<=iy2; i++)
331                {
332                    y[i] = 0;
333                }
334            }
335            else
336            {
337                for(i_=iy1; i_<=iy2;i_++)
338                {
339                    y[i_] = beta*y[i_];
340                }
341            }
342           
343            //
344            // alpha*A*x
345            //
346            for(i=i1; i<=i2; i++)
347            {
348                i1_ = (ix1)-(j1);
349                v = 0.0;
350                for(i_=j1; i_<=j2;i_++)
351                {
352                    v += a[i,i_]*x[i_+i1_];
353                }
354                y[iy1+i-i1] = y[iy1+i-i1]+alpha*v;
355            }
356        }
357        else
358        {
359           
360            //
361            // y := alpha*A'*x + beta*y;
362            //
363            if( i1>i2 | j1>j2 )
364            {
365                return;
366            }
367            System.Diagnostics.Debug.Assert(i2-i1==ix2-ix1, "MatrixVectorMultiply: A and X dont match!");
368            System.Diagnostics.Debug.Assert(j2-j1==iy2-iy1, "MatrixVectorMultiply: A and Y dont match!");
369           
370            //
371            // beta*y
372            //
373            if( beta==0 )
374            {
375                for(i=iy1; i<=iy2; i++)
376                {
377                    y[i] = 0;
378                }
379            }
380            else
381            {
382                for(i_=iy1; i_<=iy2;i_++)
383                {
384                    y[i_] = beta*y[i_];
385                }
386            }
387           
388            //
389            // alpha*A'*x
390            //
391            for(i=i1; i<=i2; i++)
392            {
393                v = alpha*x[ix1+i-i1];
394                i1_ = (j1) - (iy1);
395                for(i_=iy1; i_<=iy2;i_++)
396                {
397                    y[i_] = y[i_] + v*a[i,i_+i1_];
398                }
399            }
400        }
401    }
402
403
404    public static double pythag2(double x,
405        double y)
406    {
407        double result = 0;
408        double w = 0;
409        double xabs = 0;
410        double yabs = 0;
411        double z = 0;
412
413        xabs = Math.Abs(x);
414        yabs = Math.Abs(y);
415        w = Math.Max(xabs, yabs);
416        z = Math.Min(xabs, yabs);
417        if( z==0 )
418        {
419            result = w;
420        }
421        else
422        {
423            result = w*Math.Sqrt(1+AP.Math.Sqr(z/w));
424        }
425        return result;
426    }
427
428
429    public static void matrixmatrixmultiply(ref double[,] a,
430        int ai1,
431        int ai2,
432        int aj1,
433        int aj2,
434        bool transa,
435        ref double[,] b,
436        int bi1,
437        int bi2,
438        int bj1,
439        int bj2,
440        bool transb,
441        double alpha,
442        ref double[,] c,
443        int ci1,
444        int ci2,
445        int cj1,
446        int cj2,
447        double beta,
448        ref double[] work)
449    {
450        int arows = 0;
451        int acols = 0;
452        int brows = 0;
453        int bcols = 0;
454        int crows = 0;
455        int ccols = 0;
456        int i = 0;
457        int j = 0;
458        int k = 0;
459        int l = 0;
460        int r = 0;
461        double v = 0;
462        int i_ = 0;
463        int i1_ = 0;
464
465       
466        //
467        // Setup
468        //
469        if( !transa )
470        {
471            arows = ai2-ai1+1;
472            acols = aj2-aj1+1;
473        }
474        else
475        {
476            arows = aj2-aj1+1;
477            acols = ai2-ai1+1;
478        }
479        if( !transb )
480        {
481            brows = bi2-bi1+1;
482            bcols = bj2-bj1+1;
483        }
484        else
485        {
486            brows = bj2-bj1+1;
487            bcols = bi2-bi1+1;
488        }
489        System.Diagnostics.Debug.Assert(acols==brows, "MatrixMatrixMultiply: incorrect matrix sizes!");
490        if( arows<=0 | acols<=0 | brows<=0 | bcols<=0 )
491        {
492            return;
493        }
494        crows = arows;
495        ccols = bcols;
496       
497        //
498        // Test WORK
499        //
500        i = Math.Max(arows, acols);
501        i = Math.Max(brows, i);
502        i = Math.Max(i, bcols);
503        work[1] = 0;
504        work[i] = 0;
505       
506        //
507        // Prepare C
508        //
509        if( beta==0 )
510        {
511            for(i=ci1; i<=ci2; i++)
512            {
513                for(j=cj1; j<=cj2; j++)
514                {
515                    c[i,j] = 0;
516                }
517            }
518        }
519        else
520        {
521            for(i=ci1; i<=ci2; i++)
522            {
523                for(i_=cj1; i_<=cj2;i_++)
524                {
525                    c[i,i_] = beta*c[i,i_];
526                }
527            }
528        }
529       
530        //
531        // A*B
532        //
533        if( !transa & !transb )
534        {
535            for(l=ai1; l<=ai2; l++)
536            {
537                for(r=bi1; r<=bi2; r++)
538                {
539                    v = alpha*a[l,aj1+r-bi1];
540                    k = ci1+l-ai1;
541                    i1_ = (bj1) - (cj1);
542                    for(i_=cj1; i_<=cj2;i_++)
543                    {
544                        c[k,i_] = c[k,i_] + v*b[r,i_+i1_];
545                    }
546                }
547            }
548            return;
549        }
550       
551        //
552        // A*B'
553        //
554        if( !transa & transb )
555        {
556            if( arows*acols<brows*bcols )
557            {
558                for(r=bi1; r<=bi2; r++)
559                {
560                    for(l=ai1; l<=ai2; l++)
561                    {
562                        i1_ = (bj1)-(aj1);
563                        v = 0.0;
564                        for(i_=aj1; i_<=aj2;i_++)
565                        {
566                            v += a[l,i_]*b[r,i_+i1_];
567                        }
568                        c[ci1+l-ai1,cj1+r-bi1] = c[ci1+l-ai1,cj1+r-bi1]+alpha*v;
569                    }
570                }
571                return;
572            }
573            else
574            {
575                for(l=ai1; l<=ai2; l++)
576                {
577                    for(r=bi1; r<=bi2; r++)
578                    {
579                        i1_ = (bj1)-(aj1);
580                        v = 0.0;
581                        for(i_=aj1; i_<=aj2;i_++)
582                        {
583                            v += a[l,i_]*b[r,i_+i1_];
584                        }
585                        c[ci1+l-ai1,cj1+r-bi1] = c[ci1+l-ai1,cj1+r-bi1]+alpha*v;
586                    }
587                }
588                return;
589            }
590        }
591       
592        //
593        // A'*B
594        //
595        if( transa & !transb )
596        {
597            for(l=aj1; l<=aj2; l++)
598            {
599                for(r=bi1; r<=bi2; r++)
600                {
601                    v = alpha*a[ai1+r-bi1,l];
602                    k = ci1+l-aj1;
603                    i1_ = (bj1) - (cj1);
604                    for(i_=cj1; i_<=cj2;i_++)
605                    {
606                        c[k,i_] = c[k,i_] + v*b[r,i_+i1_];
607                    }
608                }
609            }
610            return;
611        }
612       
613        //
614        // A'*B'
615        //
616        if( transa & transb )
617        {
618            if( arows*acols<brows*bcols )
619            {
620                for(r=bi1; r<=bi2; r++)
621                {
622                    for(i=1; i<=crows; i++)
623                    {
624                        work[i] = 0.0;
625                    }
626                    for(l=ai1; l<=ai2; l++)
627                    {
628                        v = alpha*b[r,bj1+l-ai1];
629                        k = cj1+r-bi1;
630                        i1_ = (aj1) - (1);
631                        for(i_=1; i_<=crows;i_++)
632                        {
633                            work[i_] = work[i_] + v*a[l,i_+i1_];
634                        }
635                    }
636                    i1_ = (1) - (ci1);
637                    for(i_=ci1; i_<=ci2;i_++)
638                    {
639                        c[i_,k] = c[i_,k] + work[i_+i1_];
640                    }
641                }
642                return;
643            }
644            else
645            {
646                for(l=aj1; l<=aj2; l++)
647                {
648                    k = ai2-ai1+1;
649                    i1_ = (ai1) - (1);
650                    for(i_=1; i_<=k;i_++)
651                    {
652                        work[i_] = a[i_+i1_,l];
653                    }
654                    for(r=bi1; r<=bi2; r++)
655                    {
656                        i1_ = (bj1)-(1);
657                        v = 0.0;
658                        for(i_=1; i_<=k;i_++)
659                        {
660                            v += work[i_]*b[r,i_+i1_];
661                        }
662                        c[ci1+l-aj1,cj1+r-bi1] = c[ci1+l-aj1,cj1+r-bi1]+alpha*v;
663                    }
664                }
665                return;
666            }
667        }
668    }
669}
Note: See TracBrowser for help on using the repository browser.