Changeset 4068 for trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.5.0/ALGLIB-2.5.0/hblas.cs
- Timestamp:
- 07/22/10 00:44:01 (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.5.0/ALGLIB-2.5.0/hblas.cs
r3839 r4068 25 25 *************************************************************************/ 26 26 27 using System;28 27 29 namespace 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; 28 namespace alglib { 29 public class hblas { 30 public static void hermitianmatrixvectormultiply(ref AP.Complex[,] a, 31 bool isupper, 32 int i1, 33 int i2, 34 ref AP.Complex[] x, 35 AP.Complex alpha, 36 ref AP.Complex[] y) { 37 int i = 0; 38 int ba1 = 0; 39 int ba2 = 0; 40 int by1 = 0; 41 int by2 = 0; 42 int bx1 = 0; 43 int bx2 = 0; 44 int n = 0; 45 AP.Complex v = 0; 46 int i_ = 0; 47 int i1_ = 0; 52 48 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 } 49 n = i2 - i1 + 1; 50 if (n <= 0) { 51 return; 52 } 53 54 // 55 // Let A = L + D + U, where 56 // L is strictly lower triangular (main diagonal is zero) 57 // D is diagonal 58 // U is strictly upper triangular (main diagonal is zero) 59 // 60 // A*x = L*x + D*x + U*x 61 // 62 // Calculate D*x first 63 // 64 for (i = i1; i <= i2; i++) { 65 y[i - i1 + 1] = a[i, i] * x[i - i1 + 1]; 66 } 67 68 // 69 // Add L*x + U*x 70 // 71 if (isupper) { 72 for (i = i1; i <= i2 - 1; i++) { 73 74 // 75 // Add L*x to the result 76 // 77 v = x[i - i1 + 1]; 78 by1 = i - i1 + 2; 79 by2 = n; 80 ba1 = i + 1; 81 ba2 = i2; 82 i1_ = (ba1) - (by1); 83 for (i_ = by1; i_ <= by2; i_++) { 84 y[i_] = y[i_] + v * AP.Math.Conj(a[i, i_ + i1_]); 85 } 86 87 // 88 // Add U*x to the result 89 // 90 bx1 = i - i1 + 2; 91 bx2 = n; 92 ba1 = i + 1; 93 ba2 = i2; 94 i1_ = (ba1) - (bx1); 95 v = 0.0; 96 for (i_ = bx1; i_ <= bx2; i_++) { 97 v += x[i_] * a[i, i_ + i1_]; 98 } 99 y[i - i1 + 1] = y[i - i1 + 1] + v; 151 100 } 101 } else { 102 for (i = i1 + 1; i <= i2; i++) { 103 104 // 105 // Add L*x to the result 106 // 107 bx1 = 1; 108 bx2 = i - i1; 109 ba1 = i1; 110 ba2 = i - 1; 111 i1_ = (ba1) - (bx1); 112 v = 0.0; 113 for (i_ = bx1; i_ <= bx2; i_++) { 114 v += x[i_] * a[i, i_ + i1_]; 115 } 116 y[i - i1 + 1] = y[i - i1 + 1] + v; 117 118 // 119 // Add U*x to the result 120 // 121 v = x[i - i1 + 1]; 122 by1 = 1; 123 by2 = i - i1; 124 ba1 = i1; 125 ba2 = i - 1; 126 i1_ = (ba1) - (by1); 127 for (i_ = by1; i_ <= by2; i_++) { 128 y[i_] = y[i_] + v * AP.Math.Conj(a[i, i_ + i1_]); 129 } 130 } 131 } 132 for (i_ = 1; i_ <= n; i_++) { 133 y[i_] = alpha * y[i_]; 134 } 135 } 152 136 153 137 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; 138 public static void hermitianrank2update(ref AP.Complex[,] a, 139 bool isupper, 140 int i1, 141 int i2, 142 ref AP.Complex[] x, 143 ref AP.Complex[] y, 144 ref AP.Complex[] t, 145 AP.Complex alpha) { 146 int i = 0; 147 int tp1 = 0; 148 int tp2 = 0; 149 AP.Complex v = 0; 150 int i_ = 0; 151 int i1_ = 0; 169 152 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 } 153 if (isupper) { 154 for (i = i1; i <= i2; i++) { 155 tp1 = i + 1 - i1; 156 tp2 = i2 - i1 + 1; 157 v = alpha * x[i + 1 - i1]; 158 for (i_ = tp1; i_ <= tp2; i_++) { 159 t[i_] = v * AP.Math.Conj(y[i_]); 160 } 161 v = AP.Math.Conj(alpha) * y[i + 1 - i1]; 162 for (i_ = tp1; i_ <= tp2; i_++) { 163 t[i_] = t[i_] + v * AP.Math.Conj(x[i_]); 164 } 165 i1_ = (tp1) - (i); 166 for (i_ = i; i_ <= i2; i_++) { 167 a[i, i_] = a[i, i_] + t[i_ + i1_]; 168 } 216 169 } 170 } else { 171 for (i = i1; i <= i2; i++) { 172 tp1 = 1; 173 tp2 = i + 1 - i1; 174 v = alpha * x[i + 1 - i1]; 175 for (i_ = tp1; i_ <= tp2; i_++) { 176 t[i_] = v * AP.Math.Conj(y[i_]); 177 } 178 v = AP.Math.Conj(alpha) * y[i + 1 - i1]; 179 for (i_ = tp1; i_ <= tp2; i_++) { 180 t[i_] = t[i_] + v * AP.Math.Conj(x[i_]); 181 } 182 i1_ = (tp1) - (i1); 183 for (i_ = i1; i_ <= i; i_++) { 184 a[i, i_] = a[i, i_] + t[i_ + i1_]; 185 } 186 } 187 } 217 188 } 189 } 218 190 }
Note: See TracChangeset
for help on using the changeset viewer.
