Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/ALGLIB/sblas.cs @ 2575

Last change on this file since 2575 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.6 KB
Line 
1/*************************************************************************
2Copyright (c) 1992-2007 The University of Tennessee.  All rights reserved.
3
4Contributors:
5    * Sergey Bochkanov (ALGLIB project). Translation from FORTRAN to
6      pseudocode.
7
8See subroutines comments for additional copyrights.
9
10>>> SOURCE LICENSE >>>
11This program is free software; you can redistribute it and/or modify
12it under the terms of the GNU General Public License as published by
13the Free Software Foundation (www.fsf.org); either version 2 of the
14License, or (at your option) any later version.
15
16This program is distributed in the hope that it will be useful,
17but WITHOUT ANY WARRANTY; without even the implied warranty of
18MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
19GNU General Public License for more details.
20
21A copy of the GNU General Public License is available at
22http://www.fsf.org/licensing/licenses
23
24>>> END OF LICENSE >>>
25*************************************************************************/
26
27using System;
28
29namespace alglib
30{
31    public class sblas
32    {
33        public static void symmetricmatrixvectormultiply(ref double[,] a,
34            bool isupper,
35            int i1,
36            int i2,
37            ref double[] x,
38            double alpha,
39            ref double[] y)
40        {
41            int i = 0;
42            int ba1 = 0;
43            int ba2 = 0;
44            int by1 = 0;
45            int by2 = 0;
46            int bx1 = 0;
47            int bx2 = 0;
48            int n = 0;
49            double v = 0;
50            int i_ = 0;
51            int i1_ = 0;
52
53            n = i2-i1+1;
54            if( n<=0 )
55            {
56                return;
57            }
58           
59            //
60            // Let A = L + D + U, where
61            //  L is strictly lower triangular (main diagonal is zero)
62            //  D is diagonal
63            //  U is strictly upper triangular (main diagonal is zero)
64            //
65            // A*x = L*x + D*x + U*x
66            //
67            // Calculate D*x first
68            //
69            for(i=i1; i<=i2; i++)
70            {
71                y[i-i1+1] = a[i,i]*x[i-i1+1];
72            }
73           
74            //
75            // Add L*x + U*x
76            //
77            if( isupper )
78            {
79                for(i=i1; i<=i2-1; i++)
80                {
81                   
82                    //
83                    // Add L*x to the result
84                    //
85                    v = x[i-i1+1];
86                    by1 = i-i1+2;
87                    by2 = n;
88                    ba1 = i+1;
89                    ba2 = i2;
90                    i1_ = (ba1) - (by1);
91                    for(i_=by1; i_<=by2;i_++)
92                    {
93                        y[i_] = y[i_] + v*a[i,i_+i1_];
94                    }
95                   
96                    //
97                    // Add U*x to the result
98                    //
99                    bx1 = i-i1+2;
100                    bx2 = n;
101                    ba1 = i+1;
102                    ba2 = i2;
103                    i1_ = (ba1)-(bx1);
104                    v = 0.0;
105                    for(i_=bx1; i_<=bx2;i_++)
106                    {
107                        v += x[i_]*a[i,i_+i1_];
108                    }
109                    y[i-i1+1] = y[i-i1+1]+v;
110                }
111            }
112            else
113            {
114                for(i=i1+1; i<=i2; i++)
115                {
116                   
117                    //
118                    // Add L*x to the result
119                    //
120                    bx1 = 1;
121                    bx2 = i-i1;
122                    ba1 = i1;
123                    ba2 = i-1;
124                    i1_ = (ba1)-(bx1);
125                    v = 0.0;
126                    for(i_=bx1; i_<=bx2;i_++)
127                    {
128                        v += x[i_]*a[i,i_+i1_];
129                    }
130                    y[i-i1+1] = y[i-i1+1]+v;
131                   
132                    //
133                    // Add U*x to the result
134                    //
135                    v = x[i-i1+1];
136                    by1 = 1;
137                    by2 = i-i1;
138                    ba1 = i1;
139                    ba2 = i-1;
140                    i1_ = (ba1) - (by1);
141                    for(i_=by1; i_<=by2;i_++)
142                    {
143                        y[i_] = y[i_] + v*a[i,i_+i1_];
144                    }
145                }
146            }
147            for(i_=1; i_<=n;i_++)
148            {
149                y[i_] = alpha*y[i_];
150            }
151        }
152
153
154        public static void symmetricrank2update(ref double[,] a,
155            bool isupper,
156            int i1,
157            int i2,
158            ref double[] x,
159            ref double[] y,
160            ref double[] t,
161            double alpha)
162        {
163            int i = 0;
164            int tp1 = 0;
165            int tp2 = 0;
166            double v = 0;
167            int i_ = 0;
168            int i1_ = 0;
169
170            if( isupper )
171            {
172                for(i=i1; i<=i2; i++)
173                {
174                    tp1 = i+1-i1;
175                    tp2 = i2-i1+1;
176                    v = x[i+1-i1];
177                    for(i_=tp1; i_<=tp2;i_++)
178                    {
179                        t[i_] = v*y[i_];
180                    }
181                    v = y[i+1-i1];
182                    for(i_=tp1; i_<=tp2;i_++)
183                    {
184                        t[i_] = t[i_] + v*x[i_];
185                    }
186                    for(i_=tp1; i_<=tp2;i_++)
187                    {
188                        t[i_] = alpha*t[i_];
189                    }
190                    i1_ = (tp1) - (i);
191                    for(i_=i; i_<=i2;i_++)
192                    {
193                        a[i,i_] = a[i,i_] + t[i_+i1_];
194                    }
195                }
196            }
197            else
198            {
199                for(i=i1; i<=i2; i++)
200                {
201                    tp1 = 1;
202                    tp2 = i+1-i1;
203                    v = x[i+1-i1];
204                    for(i_=tp1; i_<=tp2;i_++)
205                    {
206                        t[i_] = v*y[i_];
207                    }
208                    v = y[i+1-i1];
209                    for(i_=tp1; i_<=tp2;i_++)
210                    {
211                        t[i_] = t[i_] + v*x[i_];
212                    }
213                    for(i_=tp1; i_<=tp2;i_++)
214                    {
215                        t[i_] = alpha*t[i_];
216                    }
217                    i1_ = (tp1) - (i1);
218                    for(i_=i1; i_<=i;i_++)
219                    {
220                        a[i,i_] = a[i,i_] + t[i_+i1_];
221                    }
222                }
223            }
224        }
225    }
226}
Note: See TracBrowser for help on using the repository browser.