Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/ALGLIB/xblas.cs @ 2567

Last change on this file since 2567 was 2563, checked in by gkronber, 15 years ago

Updated ALGLIB to latest version. #751 (Plugin for for data-modeling with ANN (integrated into CEDMA))

File size: 6.1 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/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            int k = 0;
38            int ks = 0;
39            double mx = 0;
40            double v = 0;
41            double v1 = 0;
42            double v2 = 0;
43            double s = 0;
44            double ln2 = 0;
45            double chunk = 0;
46            double invchunk = 0;
47            bool allzeros = new bool();
48            int i_ = 0;
49
50           
51            //
52            // special cases:
53            // * N=0
54            // * N is too large to use integer arithmetics
55            //
56            if( n==0 )
57            {
58                r = 0;
59                rerr = 0;
60                return;
61            }
62            System.Diagnostics.Debug.Assert(n<536870912, "XDot: N is too large!");
63           
64            //
65            // Prepare
66            //
67            ln2 = Math.Log(2);
68           
69            //
70            // calculate pairwise products vector TEMP
71            // (relative precision of TEMP - almost full)
72            // find infinity-norm of products vector
73            //
74            mx = 0;
75            for(i=0; i<=n-1; i++)
76            {
77                v = a[i]*b[i];
78                temp[i] = v;
79                mx = Math.Max(mx, Math.Abs(v));
80            }
81            if( (double)(mx)==(double)(0) )
82            {
83                r = 0;
84                rerr = 0;
85                return;
86            }
87            rerr = mx*AP.Math.MachineEpsilon;
88           
89            //
90            // 1. find S such that 0.5<=S*MX<1
91            // 2. multiply TEMP by S, so task is normalized in some sense
92            // 3. S:=1/S so we can obtain original vector multiplying by S
93            //
94            k = (int)Math.Round(Math.Log(mx)/ln2);
95            s = xfastpow(2, -k);
96            while( (double)(s*mx)>=(double)(1) )
97            {
98                s = 0.5*s;
99            }
100            while( (double)(s*mx)<(double)(0.5) )
101            {
102                s = 2*s;
103            }
104            for(i_=0; i_<=n-1;i_++)
105            {
106                temp[i_] = s*temp[i_];
107            }
108            s = 1/s;
109           
110            //
111            // find Chunk=2^M such that N*Chunk<2^29
112            //
113            // we have chosen upper limit (2^29) with enough space left
114            // to tolerate possible problems with rounding and N's close
115            // to the limit, so we don't want to be very strict here.
116            //
117            k = (int)(Math.Log((double)(536870912)/(double)(n))/ln2);
118            chunk = xfastpow(2, k);
119            if( (double)(chunk)<(double)(2) )
120            {
121                chunk = 2;
122            }
123            invchunk = 1/chunk;
124           
125            //
126            // calculate result
127            //
128            r = 0;
129            for(i_=0; i_<=n-1;i_++)
130            {
131                temp[i_] = chunk*temp[i_];
132            }
133            while( true )
134            {
135                s = s*invchunk;
136                allzeros = true;
137                ks = 0;
138                for(i=0; i<=n-1; i++)
139                {
140                    v = temp[i];
141                    k = (int)(v);
142                    if( (double)(v)!=(double)(k) )
143                    {
144                        allzeros = false;
145                    }
146                    temp[i] = chunk*(v-k);
147                    ks = ks+k;
148                }
149                r = r+s*ks;
150                v = Math.Abs(r);
151                if( allzeros | (double)(s*n+mx)==(double)(mx) )
152                {
153                    break;
154                }
155            }
156           
157            //
158            // correct error
159            //
160            rerr = Math.Max(rerr, Math.Abs(r)*AP.Math.MachineEpsilon);
161        }
162
163
164        private static double xfastpow(double r,
165            int n)
166        {
167            double result = 0;
168
169            if( n>0 )
170            {
171                if( n%2==0 )
172                {
173                    result = AP.Math.Sqr(xfastpow(r, n/2));
174                }
175                else
176                {
177                    result = r*xfastpow(r, n-1);
178                }
179                return result;
180            }
181            if( n==0 )
182            {
183                result = 1;
184            }
185            if( n<0 )
186            {
187                result = xfastpow(1/r, -n);
188            }
189            return result;
190        }
191
192
193        private static double xfrac(double r)
194        {
195            double result = 0;
196            int i = 0;
197
198            if( (double)(r)==(double)(0) )
199            {
200                result = 0;
201                return result;
202            }
203            if( (double)(r)<(double)(0) )
204            {
205                result = -1;
206                r = -r;
207            }
208            else
209            {
210                result = 1;
211            }
212            result = result*(r-(int)Math.Floor(r));
213            return result;
214        }
215    }
216}
Note: See TracBrowser for help on using the repository browser.