Free cookie consent management tool by TermsFeed Policy Generator

source: branches/3.2/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.1.2.2591/ALGLIB-2.1.2.2591/hblas.cs @ 4539

Last change on this file since 4539 was 2645, checked in by mkommend, 15 years ago

extracted external libraries and adapted dependent plugins (ticket #837)

File size: 6.5 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 hblas
32    {
33        public static void hermitianmatrixvectormultiply(ref AP.Complex[,] a,
34            bool isupper,
35            int i1,
36            int i2,
37            ref AP.Complex[] x,
38            AP.Complex alpha,
39            ref AP.Complex[] 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            AP.Complex 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*AP.Math.Conj(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*AP.Math.Conj(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 hermitianrank2update(ref AP.Complex[,] a,
155            bool isupper,
156            int i1,
157            int i2,
158            ref AP.Complex[] x,
159            ref AP.Complex[] y,
160            ref AP.Complex[] t,
161            AP.Complex alpha)
162        {
163            int i = 0;
164            int tp1 = 0;
165            int tp2 = 0;
166            AP.Complex 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 = alpha*x[i+1-i1];
177                    for(i_=tp1; i_<=tp2;i_++)
178                    {
179                        t[i_] = v*AP.Math.Conj(y[i_]);
180                    }
181                    v = AP.Math.Conj(alpha)*y[i+1-i1];
182                    for(i_=tp1; i_<=tp2;i_++)
183                    {
184                        t[i_] = t[i_] + v*AP.Math.Conj(x[i_]);
185                    }
186                    i1_ = (tp1) - (i);
187                    for(i_=i; i_<=i2;i_++)
188                    {
189                        a[i,i_] = a[i,i_] + t[i_+i1_];
190                    }
191                }
192            }
193            else
194            {
195                for(i=i1; i<=i2; i++)
196                {
197                    tp1 = 1;
198                    tp2 = i+1-i1;
199                    v = alpha*x[i+1-i1];
200                    for(i_=tp1; i_<=tp2;i_++)
201                    {
202                        t[i_] = v*AP.Math.Conj(y[i_]);
203                    }
204                    v = AP.Math.Conj(alpha)*y[i+1-i1];
205                    for(i_=tp1; i_<=tp2;i_++)
206                    {
207                        t[i_] = t[i_] + v*AP.Math.Conj(x[i_]);
208                    }
209                    i1_ = (tp1) - (i1);
210                    for(i_=i1; i_<=i;i_++)
211                    {
212                        a[i,i_] = a[i,i_] + t[i_+i1_];
213                    }
214                }
215            }
216        }
217    }
218}
Note: See TracBrowser for help on using the repository browser.