Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/ALGLIB/blas.cs @ 2515

Last change on this file since 2515 was 2430, checked in by gkronber, 15 years ago

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