1 9 package org.jscience.mathematics.matrices; 10 11 import org.jscience.mathematics.numbers.Numeric; 12 import java.util.Comparator ; 13 14 import javolution.realtime.LocalContext; 15 import javolution.realtime.LocalReference; 16 import javolution.realtime.PoolContext; 17 import javolution.realtime.RealtimeObject; 18 19 43 public final class LUDecomposition<O extends Operable<O>> extends 44 RealtimeObject { 45 46 51 public static final Comparator <Operable> NUMERIC_COMPARATOR = new Comparator <Operable>() { 52 53 public int compare(Operable left, Operable right) { 54 if ((left instanceof Numeric) && (right instanceof Numeric)) 55 return ((Numeric) left).isLargerThan((Numeric) right) ? 1 : -1; 56 if (left.equals(left.plus(left))) return -1; 58 if (right.equals(right.plus(right))) return 1; 60 return 0; 61 } 62 }; 63 64 67 private static final LocalReference<Comparator <Operable>> PIVOT_COMPARATOR 68 = new LocalReference<Comparator <Operable>>( 69 NUMERIC_COMPARATOR); 70 71 74 private int _n; 75 76 79 private O[] _lu; 80 81 84 private O[] _tmp; 85 86 89 private final int[] _pivots; 90 91 94 private int _permutationCount; 95 96 101 private LUDecomposition(int dimension) { 102 _lu = (O[]) new Operable[dimension * dimension]; 103 _tmp = (O[]) new Operable[dimension * dimension]; 104 _pivots = new int[dimension]; 105 } 106 107 115 static <O extends Operable<O>> LUDecomposition<O> valueOf(Matrix<O> source) { 116 if (!source.isSquare()) 117 throw new MatrixException("Matrix is not square"); 118 int dimension = source.m; 119 LUDecomposition lu; 120 if (dimension <= 1 << 2) { 121 lu = FACTORY_2.object(); 122 } else if (dimension <= 1 << 3) { 123 lu = FACTORY_3.object(); 124 } else if (dimension <= 1 << 4) { 125 lu = FACTORY_4.object(); 126 } else if (dimension <= 1 << 5) { 127 lu = FACTORY_5.object(); 128 } else if (dimension <= 1 << 6) { 129 lu = FACTORY_6.object(); 130 } else if (dimension <= 1 << 7) { 131 lu = FACTORY_7.object(); 132 } else if (dimension <= 1 << 8) { 133 lu = FACTORY_8.object(); 134 } else if (dimension <= 1 << 9) { 135 lu = FACTORY_9.object(); 136 } else if (dimension <= 1 << 10) { 137 lu = FACTORY_10.object(); 138 } else if (dimension <= 1 << 11) { 139 lu = FACTORY_11.object(); 140 } else if (dimension <= 1 << 12) { 141 lu = FACTORY_12.object(); 142 } else if (dimension <= 1 << 13) { 143 lu = FACTORY_13.object(); 144 } else if (dimension <= 1 << 14) { 145 lu = FACTORY_14.object(); 146 } else if (dimension <= 1 << 15) { 147 lu = FACTORY_15.object(); 148 } else { 149 throw new UnsupportedOperationException ("Dimension : " + dimension 150 + " too large"); 151 } 152 lu._n = dimension; 153 lu._permutationCount = 0; 154 lu.construct(source); 155 return lu; 156 157 } 158 159 166 public void setPivotComparator(Comparator <Operable> cmp) { 167 PIVOT_COMPARATOR.set(cmp); 168 } 169 170 179 public Matrix<O> solve(Matrix<O> B) { 180 if (_n != B.m) 181 throw new MatrixException("Input vector has " + B.m 182 + " rows instead of " + _n); 183 184 final int n = B.n; 186 Matrix<O> X = Matrix.newInstance(_n, n); 187 final O[] x = X.o; 188 for (int i = 0; i < _n; i++) { 189 for (int j = 0; j < n; j++) { 190 x[i * n + j] = B.o[_pivots[i] * n + j]; 191 } 192 } 193 194 for (int k = 0; k < _n; k++) { 196 for (int i = k + 1; i < _n; i++) { 197 O luik = _lu[i * _n + k]; 198 for (int j = 0; j < n; j++) { 199 x[i * n + j] = x[i * n + j].plus(luik.times(x[k * n + j]) 200 .opposite()); 201 } 202 } 203 } 204 205 for (int k = _n - 1; k >= 0; k--) { 207 for (int j = 0; j < n; j++) { 208 x[k * n + j] = (_lu[k * _n + k].reciprocal()) 209 .times(x[k * n + j]); 210 } 211 for (int i = 0; i < k; i++) { 212 O luik = _lu[i * _n + k]; 213 for (int j = 0; j < n; j++) { 214 x[i * n + j] = x[i * n + j].plus(luik.times(x[k * n + j]) 215 .opposite()); 216 } 217 } 218 } 219 return X; 220 } 221 222 229 public Matrix<O> inverse() { 230 final int n = _n; 232 Matrix<O> R = Matrix.newInstance(n, n); 233 final O[] r = R.o; 234 for (int i = 0; i < n; i++) { 235 for (int j = i; j < n; j++) { 236 r[i * n + j] = _lu[i * n + j]; 237 } 238 } 239 for (int j = n - 1; j >= 0; j--) { 240 r[j * n + j] = r[j * n + j].reciprocal(); 241 for (int i = j - 1; i >= 0; i--) { 242 O sum = r[i * n + j].times(r[j * n + j]).opposite(); 245 for (int k = j - 1; k > i; k--) { 246 sum = sum.plus(r[i * n + k].times(r[k * n + j]) 247 .opposite()); 248 } 249 r[i * n + j] = (r[i * n + i].reciprocal()) 250 .times(sum); 251 r[i * n + j].move(ObjectSpace.OUTER); 252 } 256 } 257 for (int i = 0; i < n; i++) { 259 for (int j = n - 2; j >= 0; j--) { 260 for (int k = j + 1; k < n; k++) { 263 O lukj = _lu[k * n + j]; 264 if (r[i * n + j] != null) { 265 r[i * n + j] = r[i * n + j].plus(r[i * n + k] 266 .times(lukj).opposite()); 267 } else { 268 r[i * n + j] = r[i * n + k].times(lukj).opposite(); 269 } 270 } 271 r[i * n + j].move(ObjectSpace.OUTER); 272 } 276 } 277 for (int i = 0; i < n; i++) { 279 for (int j = 0; j < n; j++) { 280 _tmp[j] = r[i * n + j]; 281 } 282 for (int j = 0; j < n; j++) { 283 r[i * n + _pivots[j]] = _tmp[j]; 284 } 285 } 286 return R; 287 } 288 289 295 public O determinant() { 296 O product = _lu[0]; 297 for (int i = 1; i < _n; i++) { 298 product = product.times(_lu[i * _n + i]); 299 } 300 return ((_permutationCount & 1) == 0) ? product : product.opposite(); 301 } 302 303 311 private void construct(Matrix<O> source) { 312 System.arraycopy(source.o, 0, _lu, 0, _n * _n); 313 for (int i = 0; i < _n; i++) { 314 _pivots[i] = i; 315 } 316 317 Comparator <Operable> cmp = PIVOT_COMPARATOR.get(); 319 final int n = _n; 320 for (int k = 0; k < _n; k++) { 321 322 int pivot = k; 326 for (int i = k + 1; i < n; i++) { 327 if (cmp.compare(_lu[i * n + k], _lu[pivot * n + k]) > 0) { 328 pivot = i; 329 } 330 } 331 if (pivot != k) { for (int j = 0; j < n; j++) { 333 O tmp = _lu[pivot * n + j]; 334 _lu[pivot * n + j] = _lu[k * n + j]; 335 _lu[k * n + j] = tmp; 336 } 337 int j = _pivots[pivot]; 338 _pivots[pivot] = _pivots[k]; 339 _pivots[k] = j; 340 _permutationCount++; 341 } 342 343 O lukkInv = _lu[k * n + k].reciprocal(); 345 for (int i = k + 1; i < n; i++) { 346 _lu[i * n + k] = _lu[i * n + k].times(lukkInv); 349 for (int j = k + 1; j < n; j++) { 350 _lu[i * n + j] = _lu[i * n + j].plus(_lu[i * n + k].times( 351 _lu[k * n + j]).opposite()); 352 } 353 } 354 } 355 } 356 357 public boolean move(ObjectSpace os) { 359 if (super.move(os)) { 360 for (int i = _n * _n; i > 0;) { 361 _lu[--i].move(os); 362 } 363 return true; 364 } 365 return false; 366 } 367 368 371 private static final Factory<LUDecomposition> FACTORY_2 = new Factory<LUDecomposition>() { 372 protected LUDecomposition create() { 373 return new LUDecomposition(1 << 2); 374 } 375 }; 376 377 private static final Factory<LUDecomposition> FACTORY_3 = new Factory<LUDecomposition>() { 378 protected LUDecomposition create() { 379 return new LUDecomposition(1 << 3); 380 } 381 }; 382 383 private static final Factory<LUDecomposition> FACTORY_4 = new Factory<LUDecomposition>() { 384 protected LUDecomposition create() { 385 return new LUDecomposition(1 << 4); 386 } 387 }; 388 389 private static final Factory<LUDecomposition> FACTORY_5 = new Factory<LUDecomposition>() { 390 protected LUDecomposition create() { 391 return new LUDecomposition(1 << 5); 392 } 393 }; 394 395 private static final Factory<LUDecomposition> FACTORY_6 = new Factory<LUDecomposition>() { 396 protected LUDecomposition create() { 397 return new LUDecomposition(1 << 6); 398 } 399 }; 400 401 private static final Factory<LUDecomposition> FACTORY_7 = new Factory<LUDecomposition>() { 402 protected LUDecomposition create() { 403 return new LUDecomposition(1 << 7); 404 } 405 }; 406 407 private static final Factory<LUDecomposition> FACTORY_8 = new Factory<LUDecomposition>() { 408 protected LUDecomposition create() { 409 return new LUDecomposition(1 << 8); 410 } 411 }; 412 413 private static final Factory<LUDecomposition> FACTORY_9 = new Factory<LUDecomposition>() { 414 protected LUDecomposition create() { 415 return new LUDecomposition(1 << 9); 416 } 417 }; 418 419 private static final Factory<LUDecomposition> FACTORY_10 = new Factory<LUDecomposition>() { 420 protected LUDecomposition create() { 421 return new LUDecomposition(1 << 10); 422 } 423 }; 424 425 private static final Factory<LUDecomposition> FACTORY_11 = new Factory<LUDecomposition>() { 426 protected LUDecomposition create() { 427 return new LUDecomposition(1 << 11); 428 } 429 }; 430 431 private static final Factory<LUDecomposition> FACTORY_12 = new Factory<LUDecomposition>() { 432 protected LUDecomposition create() { 433 return new LUDecomposition(1 << 12); 434 } 435 }; 436 437 private static final Factory<LUDecomposition> FACTORY_13 = new Factory<LUDecomposition>() { 438 protected LUDecomposition create() { 439 return new LUDecomposition(1 << 13); 440 } 441 }; 442 443 private static final Factory<LUDecomposition> FACTORY_14 = new Factory<LUDecomposition>() { 444 protected LUDecomposition create() { 445 return new LUDecomposition(1 << 14); 446 } 447 }; 448 449 private static final Factory<LUDecomposition> FACTORY_15 = new Factory<LUDecomposition>() { 450 protected LUDecomposition create() { 451 return new LUDecomposition(1 << 15); 452 } 453 }; 454 455 } | Popular Tags |