| 14 | */ |
| 15 | public final class CholeskySolver implements LinearSolver { |
| 16 | @Override |
| 17 | public double[] solve(final double[][] ain, final double[] b) { |
| 18 | final int dim = b.length; |
| 19 | double sign = 1.0; |
| 20 | for(int i=0;i<dim;++i) { |
| 21 | if(Math.abs(ain[i][i])>0.0) { |
| 22 | if(ain[i][i]<0.0) { |
| 23 | sign = -1.0; |
| 24 | } |
| 25 | break; |
| 26 | } |
| 27 | } |
| 28 | final DoubleMatrix2D ma = new DenseDoubleMatrix2D(dim,dim); |
| 29 | for(int i=0;i<dim;++i) { |
| 30 | for(int j=0;j<dim;++j) { |
| 31 | ma.set(i,j,sign*ain[i][j]); |
| 32 | } |
| 33 | } |
| 34 | final DoubleMatrix2D mb = new DenseDoubleMatrix2D(dim,1); |
| 35 | for(int i=0;i<dim;++i) { |
| 36 | mb.set(i,0,sign*b[i]); |
| 37 | } |
| 38 | final CholeskyDecomposition decomp = new CholeskyDecomposition(ma); |
| 39 | final DoubleMatrix2D mx; |
| 40 | if(decomp.isSymmetricPositiveDefinite()) { |
| 41 | mx = decomp.solve(mb); |
| 42 | } else { |
| 43 | // fall back to direct |
| 44 | mx = Algebra.ZERO.solve(ma,mb); |
| 45 | } |
| 46 | final double[] x = new double[dim]; |
| 47 | for(int i=0;i<dim;++i) { |
| 48 | x[i] = mx.get(i,0); |
| 49 | } |
| 50 | return x; |
| 51 | } |
| 52 | } |