public DoubleMatrix1D solve ( DoubleMatrix1D b ) { if ( b . size ( ) != dim ) { log . error ( "wrong dimension of vector b: expected " + dim + ", actual " + b . size ( ) ) ; throw new RuntimeException ( "wrong dimension of vector b: expected " + dim + ", actual " + b . size ( ) ) ; } // with scaling, we must solve U.Q.U.z = U.b, after that we have x = U.z if ( this . rescaler != null ) { // b = ALG.mult(this.U, b); b = ColtUtils . diagonalMatrixMult ( this . U , b ) ; } // Solve L.y = b final double [ ] y = new double [ dim ] ; // L is upper left diagonal for ( int i = 0 ; i < diagonalLength ; i ++ ) { double LII = LData [ 0 ] [ i ] ; y [ i ] = b . getQuick ( i ) / LII ; } for ( int i = diagonalLength ; i < dim ; i ++ ) { double [ ] LI = LData [ i - diagonalLength + 1 ] ; double LII = LI [ i ] ; double sum = 0 ; for ( int j = 0 ; j < i ; j ++ ) { sum += LI [ j ] * y [ j ] ; } y [ i ] = ( b . getQuick ( i ) - sum ) / LII ; } // logger.debug("b: " + ArrayUtils.toString(b)); // logger.debug("L.y: " + ArrayUtils.toString(getL().operate(y))); // Solve L[T].x = y final DoubleMatrix1D x = F1 . make ( dim ) ; // for (int i = dim-1; i > -1; i--) { // double sum = 0; // int ll = Math.max(i, diagonalLength-1); // for (int j = dim-1; j > ll; j--) { // sum += LData[j][i] * x.getQuick(j); // // } // x.setQuick(i, (y[i] - sum)/ LData[i][i]); // } for ( int i = dim - 1 ; i > diagonalLength - 1 ; i -- ) { double sum = 0 ; for ( int j = dim - 1 ; j > i ; j -- ) { sum += LData [ j - diagonalLength + 1 ] [ i ] * x . getQuick ( j ) ; } x . setQuick ( i , ( y [ i ] - sum ) / LData [ i - diagonalLength + 1 ] [ i ] ) ; } for ( int i = diagonalLength - 1 ; i > - 1 ; i -- ) { double sum = 0 ; for ( int j = dim - 1 ; j > diagonalLength - 1 ; j -- ) { sum += LData [ j - diagonalLength + 1 ] [ i ] * x . getQuick ( j ) ; } x . setQuick ( i , ( y [ i ] - sum ) / LData [ 0 ] [ i ] ) ; } if ( this . rescaler != null ) { // return ALG.mult(this.U, x); return ColtUtils . diagonalMatrixMult ( this . U , x ) ; } else { return x ; } }
Solve Q . x = b