Free cookie consent management tool by TermsFeed Policy Generator

source: branches/HeuristicLab.ALGLIB-2.5.0/ALGLIB-2.5.0/xblas.cs @ 5229

Last change on this file since 5229 was 3839, checked in by mkommend, 14 years ago

implemented first version of LR (ticket #1012)
File size: 10.0 KB
Line 
1
2using System;
3
4namespace alglib
5{
6    public class xblas
7    {
8        /*************************************************************************
9        More precise dot-product. Absolute error of  subroutine  result  is  about
10        1 ulp of max(MX,V), where:
11            MX = max( |a[i]*b[i]| )
12            V  = |(a,b)|
13
14        INPUT PARAMETERS
15            A       -   array[0..N-1], vector 1
16            B       -   array[0..N-1], vector 2
17            N       -   vectors length, N<2^29.
18            Temp    -   array[0..N-1], pre-allocated temporary storage
19
20        OUTPUT PARAMETERS
21            R       -   (A,B)
22            RErr    -   estimate of error. This estimate accounts for both  errors
23                        during  calculation  of  (A,B)  and  errors  introduced by
24                        rounding of A and B to fit in double (about 1 ulp).
25
26          -- ALGLIB --
27             Copyright 24.08.2009 by Bochkanov Sergey
28        *************************************************************************/
29        public static void xdot(ref double[] a,
30            ref double[] b,
31            int n,
32            ref double[] temp,
33            ref double r,
34            ref double rerr)
35        {
36            int i = 0;
37            double mx = 0;
38            double v = 0;
39
40           
41            //
42            // special cases:
43            // * N=0
44            //
45            if( n==0 )
46            {
47                r = 0;
48                rerr = 0;
49                return;
50            }
51            mx = 0;
52            for(i=0; i<=n-1; i++)
53            {
54                v = a[i]*b[i];
55                temp[i] = v;
56                mx = Math.Max(mx, Math.Abs(v));
57            }
58            if( (double)(mx)==(double)(0) )
59            {
60                r = 0;
61                rerr = 0;
62                return;
63            }
64            xsum(ref temp, mx, n, ref r, ref rerr);
65        }
66
67
68        /*************************************************************************
69        More precise complex dot-product. Absolute error of  subroutine  result is
70        about 1 ulp of max(MX,V), where:
71            MX = max( |a[i]*b[i]| )
72            V  = |(a,b)|
73
74        INPUT PARAMETERS
75            A       -   array[0..N-1], vector 1
76            B       -   array[0..N-1], vector 2
77            N       -   vectors length, N<2^29.
78            Temp    -   array[0..2*N-1], pre-allocated temporary storage
79
80        OUTPUT PARAMETERS
81            R       -   (A,B)
82            RErr    -   estimate of error. This estimate accounts for both  errors
83                        during  calculation  of  (A,B)  and  errors  introduced by
84                        rounding of A and B to fit in double (about 1 ulp).
85
86          -- ALGLIB --
87             Copyright 27.01.2010 by Bochkanov Sergey
88        *************************************************************************/
89        public static void xcdot(ref AP.Complex[] a,
90            ref AP.Complex[] b,
91            int n,
92            ref double[] temp,
93            ref AP.Complex r,
94            ref double rerr)
95        {
96            int i = 0;
97            double mx = 0;
98            double v = 0;
99            double rerrx = 0;
100            double rerry = 0;
101
102           
103            //
104            // special cases:
105            // * N=0
106            //
107            if( n==0 )
108            {
109                r = 0;
110                rerr = 0;
111                return;
112            }
113           
114            //
115            // calculate real part
116            //
117            mx = 0;
118            for(i=0; i<=n-1; i++)
119            {
120                v = a[i].x*b[i].x;
121                temp[2*i+0] = v;
122                mx = Math.Max(mx, Math.Abs(v));
123                v = -(a[i].y*b[i].y);
124                temp[2*i+1] = v;
125                mx = Math.Max(mx, Math.Abs(v));
126            }
127            if( (double)(mx)==(double)(0) )
128            {
129                r.x = 0;
130                rerrx = 0;
131            }
132            else
133            {
134                xsum(ref temp, mx, 2*n, ref r.x, ref rerrx);
135            }
136           
137            //
138            // calculate imaginary part
139            //
140            mx = 0;
141            for(i=0; i<=n-1; i++)
142            {
143                v = a[i].x*b[i].y;
144                temp[2*i+0] = v;
145                mx = Math.Max(mx, Math.Abs(v));
146                v = a[i].y*b[i].x;
147                temp[2*i+1] = v;
148                mx = Math.Max(mx, Math.Abs(v));
149            }
150            if( (double)(mx)==(double)(0) )
151            {
152                r.y = 0;
153                rerry = 0;
154            }
155            else
156            {
157                xsum(ref temp, mx, 2*n, ref r.y, ref rerry);
158            }
159           
160            //
161            // total error
162            //
163            if( (double)(rerrx)==(double)(0) & (double)(rerry)==(double)(0) )
164            {
165                rerr = 0;
166            }
167            else
168            {
169                rerr = Math.Max(rerrx, rerry)*Math.Sqrt(1+AP.Math.Sqr(Math.Min(rerrx, rerry)/Math.Max(rerrx, rerry)));
170            }
171        }
172
173
174        /*************************************************************************
175        Internal subroutine for extra-precise calculation of SUM(w[i]).
176
177        INPUT PARAMETERS:
178            W   -   array[0..N-1], values to be added
179                    W is modified during calculations.
180            MX  -   max(W[i])
181            N   -   array size
182           
183        OUTPUT PARAMETERS:
184            R   -   SUM(w[i])
185            RErr-   error estimate for R
186
187          -- ALGLIB --
188             Copyright 24.08.2009 by Bochkanov Sergey
189        *************************************************************************/
190        private static void xsum(ref double[] w,
191            double mx,
192            int n,
193            ref double r,
194            ref double rerr)
195        {
196            int i = 0;
197            int k = 0;
198            int ks = 0;
199            double v = 0;
200            double s = 0;
201            double ln2 = 0;
202            double chunk = 0;
203            double invchunk = 0;
204            bool allzeros = new bool();
205            int i_ = 0;
206
207           
208            //
209            // special cases:
210            // * N=0
211            // * N is too large to use integer arithmetics
212            //
213            if( n==0 )
214            {
215                r = 0;
216                rerr = 0;
217                return;
218            }
219            if( (double)(mx)==(double)(0) )
220            {
221                r = 0;
222                rerr = 0;
223                return;
224            }
225            System.Diagnostics.Debug.Assert(n<536870912, "XDot: N is too large!");
226           
227            //
228            // Prepare
229            //
230            ln2 = Math.Log(2);
231            rerr = mx*AP.Math.MachineEpsilon;
232           
233            //
234            // 1. find S such that 0.5<=S*MX<1
235            // 2. multiply W by S, so task is normalized in some sense
236            // 3. S:=1/S so we can obtain original vector multiplying by S
237            //
238            k = (int)Math.Round(Math.Log(mx)/ln2);
239            s = xfastpow(2, -k);
240            while( (double)(s*mx)>=(double)(1) )
241            {
242                s = 0.5*s;
243            }
244            while( (double)(s*mx)<(double)(0.5) )
245            {
246                s = 2*s;
247            }
248            for(i_=0; i_<=n-1;i_++)
249            {
250                w[i_] = s*w[i_];
251            }
252            s = 1/s;
253           
254            //
255            // find Chunk=2^M such that N*Chunk<2^29
256            //
257            // we have chosen upper limit (2^29) with enough space left
258            // to tolerate possible problems with rounding and N's close
259            // to the limit, so we don't want to be very strict here.
260            //
261            k = (int)(Math.Log((double)(536870912)/(double)(n))/ln2);
262            chunk = xfastpow(2, k);
263            if( (double)(chunk)<(double)(2) )
264            {
265                chunk = 2;
266            }
267            invchunk = 1/chunk;
268           
269            //
270            // calculate result
271            //
272            r = 0;
273            for(i_=0; i_<=n-1;i_++)
274            {
275                w[i_] = chunk*w[i_];
276            }
277            while( true )
278            {
279                s = s*invchunk;
280                allzeros = true;
281                ks = 0;
282                for(i=0; i<=n-1; i++)
283                {
284                    v = w[i];
285                    k = (int)(v);
286                    if( (double)(v)!=(double)(k) )
287                    {
288                        allzeros = false;
289                    }
290                    w[i] = chunk*(v-k);
291                    ks = ks+k;
292                }
293                r = r+s*ks;
294                v = Math.Abs(r);
295                if( allzeros | (double)(s*n+mx)==(double)(mx) )
296                {
297                    break;
298                }
299            }
300           
301            //
302            // correct error
303            //
304            rerr = Math.Max(rerr, Math.Abs(r)*AP.Math.MachineEpsilon);
305        }
306
307
308        /*************************************************************************
309        Fast Pow
310
311          -- ALGLIB --
312             Copyright 24.08.2009 by Bochkanov Sergey
313        *************************************************************************/
314        private static double xfastpow(double r,
315            int n)
316        {
317            double result = 0;
318
319            if( n>0 )
320            {
321                if( n%2==0 )
322                {
323                    result = AP.Math.Sqr(xfastpow(r, n/2));
324                }
325                else
326                {
327                    result = r*xfastpow(r, n-1);
328                }
329                return result;
330            }
331            if( n==0 )
332            {
333                result = 1;
334            }
335            if( n<0 )
336            {
337                result = xfastpow(1/r, -n);
338            }
339            return result;
340        }
341    }
342}
Note: See TracBrowser for help on using the repository browser.