Solve A*X = B.

@param Matrix $B Right hand side
@return Matrix ... Solution if A is square, least squares solution otherwise
970123
public function solve($B)
    {
        if ($this->m == $this->n) {
            $LU = new PHPExcel_Shared_JAMA_LUDecomposition($this);
            return $LU->solve($B);
        } else {
            $QR = new PHPExcel_Shared_JAMA_QRDecomposition($this);
            return $QR->solve($B);
        }
    }
1302731
public DenseMatrix solve(DenseMatrix B) throws MatrixNotSPDException {
        if (notspd)
            throw new MatrixNotSPDException();
        if (n != B.numRows())
            throw new IllegalArgumentException("n != B.numRows()");

        intW info = new intW(0);
        if (upper)
            LAPACK.getInstance().dpotrs(UpLo.Upper.netlib(), Cu.numRows(),
                    B.numColumns(), Cu.getData(), Matrices.ld(Cu.numRows()),
                    B.getData(), Matrices.ld(Cu.numRows()), info);
        else
            LAPACK.getInstance().dpotrs(UpLo.Lower.netlib(), Cl.numRows(),
                    B.numColumns(), Cl.getData(), Matrices.ld(Cl.numRows()),
                    B.getData(), Matrices.ld(Cl.numRows()), info);

        if (info.val < 0)
            throw new IllegalArgumentException();

        return B;
    }
201353
public Matrix plus(final Matrix B)
    {
        final Matrix A = this;

        if ((B.m_rows != A.m_rows) || (B.m_columns != A.m_columns))
        {
            throw new GeometryException("Illegal matrix dimensions");
        }
        final Matrix C = new Matrix(m_rows, m_columns);

        for (int i = 0; i < m_rows; i++)
        {
            for (int j = 0; j < m_columns; j++)
            {
                C.m_data[i][j] = A.m_data[i][j] + B.m_data[i][j];
            }
        }
        return C;
    }
86313
@Override
    public void solve(DMatrixRMaj B, DMatrixRMaj X) {
        if( B.numRows != numRows )
            throw new IllegalArgumentException("Unexpected dimensions for X: X rows = "+X.numRows+" expected = "+numCols);
        X.reshape(numCols,B.numCols);

        int BnumCols = B.numCols;

        // solve each column one by one
        for( int colB = 0; colB < BnumCols; colB++ ) {

            // make a copy of this column in the vector
            for( int i = 0; i < numRows; i++ ) {
                a[i] = B.data[i*BnumCols + colB];
            }

            // Solve Qa=b
            // a = Q'b
            // a = Q_{n-1}...Q_2*Q_1*b
            //
            // Q_n*b = (I-gamma*u*u^T)*b = b - u*(gamma*U^T*b)
            for( int n = 0; n < numCols; n++ ) {
                u[n] = 1;
                double ub = a[n];
                // U^T*b
                for( int i = n+1; i < numRows; i++ ) {
                    ub += (u[i] = QR.unsafe_get(i,n))*a[i];
                }

                // gamma*U^T*b
                ub *= gammas[n];
 
                for( int i = n; i < numRows; i++ ) {
                    a[i] -= u[i]*ub;
                }
            }

            // solve for Rx = b using the standard upper triangular solver
            TriangularSolver_DDRM.solveU(QR.data,a,numCols);

            // save the results
            for( int i = 0; i < numCols; i++ ) {
                X.data[i*X.numCols+colB] = a[i];
            }
        }
    }
742438
public double[][] solve(double[][] B) {
    if(B.length != L.length) {
      throw new IllegalArgumentException(ERR_MATRIX_DIMENSIONS);
    }
    if(!isspd) {
      throw new ArithmeticException(ERR_MATRIX_NOT_SPD);
    }
    // Work on a copy!
    return solveLtransposed(solveL(copy(B)));
  }
500024
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];
		for (int i = 0; i < dim; i++) {
			double[] LI = LData[i];
			double sum = 0;
			for (int j = 0; j < i; j++) {
				sum += LI[j] * y[j];
			}
			y[i] = (b.getQuick(i) - sum) / LI[i];
		}

		// logger.debug("b: " + ArrayUtils.toString(b));
		// logger.debug("L.y: " + ArrayUtils.toString(getL().operate(y)));

		// Solve D.z = y
		final double[] z = new double[dim];
		for (int i = 0; i < dim; i++) {
			z[i] = y[i] / DData[i];
		}

		// Solve L[T].x = z
		final DoubleMatrix1D x = F1.make(dim);
		for (int i = dim - 1; i > -1; i--) {
			double sum = 0;
			for (int j = dim - 1; j > i; j--) {
				sum += LData[j][i] * x.getQuick(j);
			}
			x.setQuick(i, (z[i] - sum) / LData[i][i]);
		}

		if (this.rescaler != null) {
			// return ALG.mult(this.U, x);
			return ColtUtils.diagonalMatrixMult(this.U, x);
		} else {
			return x;
		}
	}
832854
public function augmentAbove(Matrix $b) : self
    {
        if ($this->m > 0 and $b->n() !== $this->n) {
            throw new DimensionalityMismatchException('Matrix A requires'
                . " $this->n columns but Matrix B has {$b->n()}.");
        }

        $b = array_merge($b->asArray(), $this->a);

        return self::quick($b);
    }
433358
public function equal($b)
    {
        switch (true) {
            case $b instanceof Vector:
                $c = $this->equalVector($b);
                break;

            case $b instanceof Matrix:
                $c = $this->equalMatrix($b);
                break;

            case is_int($b) or is_float($b):
                $c = $this->equalScalar($b);
                break;

            default:
                throw new InvalidArgumentException('Cannot compare vector'
                    . ' to a ' . gettype($b) . '.');
        }

        return $c;
    }
948293
public function solve(Vector $b) : ColumnVector
    {
        $k = $this->m - 1;

        [$l, $u, $p] = $this->lu();

        $pb = $p->dot($b);

        $y = [$pb[0] / ($l[0][0] ?: self::EPSILON)];

        for ($i = 1; $i < $this->m; $i++) {
            $sigma = 0;

            for ($j = 0; $j <= $i - 1; $j++) {
                $sigma += $l[$i][$j] * $y[$j];
            }

            $y[] = ($pb[$i] - $sigma) / $l[$i][$i];
        }

        $x = [];

        $x = [$k => $y[$k] / ($l[$k][$k] ?: self::EPSILON)];

        for ($i = $this->m - 2; $i >= 0; $i--) {
            $sigma = 0;

            for ($j = $i + 1; $j < $this->m; $j++) {
                $sigma += $u[$i][$j] * $x[$j];
            }

            $x[$i] = ($y[$i] - $sigma) / $u[$i][$i];
        }

        return ColumnVector::quick(array_reverse($x));
    }
614787
def solve(self, b):
        
        #b = self.permute_vec(b)
        y = []
        for bri in range(self.N):  # block row index
            for li in range(self.n):  # local row index
                s = 0.0
                for lci in range(li): # local column index
                    s += self.lu[bri][li, lci]*y[bri*self.n+lci]
                for di in range(1, self.ndiag+1):
                    if bri >= di:
                        # di:th sub diagonal (counted as distance from main diag)
                        ci = self.colbyrow[bri][li]
                        s += self.sub[di-1][(bri-di)*self.n+ci]*y[
                            (bri-di)*self.n + ci]
                y.append(b[bri*self.n+self.rowbycol[bri][li]]-s) # Doolittle: L[i, i] == 1
        x = [0]*len(y)
        for bri in range(self.N-1, -1, -1):
            for li in range(self.n - 1, -1, -1):
                s = 0.0
                for ci in range(li+1, self.n):
                    s += self.lu[bri][li, ci]*x[bri*self.n + ci]
                for di in range(1, self.ndiag+1):
                    if bri < self.N-di:
                        ci = self.colbyrow[bri][li]
                        s += self.sup[di-1][bri*self.n+ci]*x[(bri+di)*self.n + ci]
                x[bri*self.n+li] = (y[bri*self.n + li] - s)/self.lu[bri][li, li]
        return x
1021860
public function compute($a, $b)
    {
        $score = 2 * Product::scalar($a, $b);
        $squares = Product::scalar($a, $a) + Product::scalar($b, $b);
        $result = exp(-$this->gamma * ($squares - $score));

        return $result;
    }
792701
@Override
    public void solve(DMatrixRMaj B, DMatrixRMaj X) {
        if( B.numRows != numRows )
            throw new IllegalArgumentException("Unexpected dimensions for X: X rows = "+X.numRows+" expected = "+numCols);
        X.reshape(numCols,B.numCols);

        U = decomposer.getR(U,true);
        final double gammas[] = decomposer.getGammas();
        final double dataQR[] = QR.data;

        final int BnumCols = B.numCols;

        // solve each column one by one
        for( int colB = 0; colB < BnumCols; colB++ ) {

            // make a copy of this column in the vector
            for( int i = 0; i < numRows; i++ ) {
                a[i] = B.data[i*BnumCols + colB];
            }

            // Solve Qa=b
            // a = Q'b
            // a = Q_{n-1}...Q_2*Q_1*b
            //
            // Q_n*b = (I-gamma*u*u^T)*b = b - u*(gamma*U^T*b)
            for( int n = 0; n < numCols; n++ ) {
                int indexU = n*numRows + n + 1;

                double ub = a[n];
                // U^T*b
                for( int i = n+1; i < numRows; i++ , indexU++ ) {
                    ub += dataQR[indexU]*a[i];
                }

                // gamma*U^T*b
                ub *= gammas[n];

                a[n] -= ub;
                indexU = n*numRows + n + 1;
                for( int i = n+1; i < numRows; i++ , indexU++) {
                    a[i] -= dataQR[indexU]*ub;
                }
            }

            // solve for Rx = b using the standard upper triangular solver
            TriangularSolver_DDRM.solveU(U.data,a,numCols);

            // save the results
            for( int i = 0; i < numCols; i++ ) {
                X.data[i*X.numCols+colB] = a[i];
            }
        }
    }
775333
@Override
    public void solve(ZMatrixRMaj B, ZMatrixRMaj X) {
        if( X.numRows != numCols )
            throw new IllegalArgumentException("Unexpected dimensions for X: X rows = "+X.numRows+" expected = "+numCols);
        else if( B.numRows != numRows || B.numCols != X.numCols )
            throw new IllegalArgumentException("Unexpected dimensions for B");

        int BnumCols = B.numCols;
        
        // solve each column one by one
        for( int colB = 0; colB < BnumCols; colB++ ) {

            // make a copy of this column in the vector
            for( int i = 0; i < numRows; i++ ) {
                int indexB = (i*BnumCols + colB)*2;
                a.data[i*2]   = B.data[indexB];
                a.data[i*2+1] = B.data[indexB+1];
            }

            // Solve Qa=b
            // a = Q'b
            // a = Q_{n-1}...Q_2*Q_1*b
            //
            // Q_n*b = (I-gamma*u*u^T)*b = b - u*(gamma*U^T*b)
            for( int n = 0; n < numCols; n++ ) {
                double []u = QR[n];

                double realVV = u[n*2];
                double imagVV = u[n*2+1];

                u[n*2]   = 1;
                u[n*2+1] = 0;

                QrHelperFunctions_ZDRM.rank1UpdateMultR(a, u,0, gammas[n], 0, n, numRows, temp.data);

                u[n*2]   = realVV;
                u[n*2+1] = imagVV;
            }

            // solve for Rx = b using the standard upper triangular solver
            TriangularSolver_ZDRM.solveU(R.data, a.data, numCols);

            // save the results
            for( int i = 0; i < numCols; i++ ) {
                int indexB = (i*BnumCols + colB)*2;
                X.data[indexB]   = a.data[i*2];
                X.data[indexB+1] = a.data[i*2+1];
            }
        }
    }
675014
public function multiply($b)
    {
        switch (true) {
            case $b instanceof Vector:
                $c = $this->multiplyVector($b);
                break;

            case $b instanceof Matrix:
                $c = $this->multiplyMatrix($b);
                break;

            case is_int($b) or is_float($b):
                $c = $this->multiplyScalar($b);
                break;

            default:
                throw new InvalidArgumentException('Cannot multiply vector'
                    . ' with a ' . gettype($b) . '.');
        }

        return $c;
    }
600916
public function checkMatrixDimensions($B = null)
    {
        if (is_a($B, 'Matrix')) {
            if (($this->m == $B->m) && ($this->n == $B->n)) {
                return true;
            } else {
                trigger_error(MatrixDimensionException, ERROR);
            }
        } else {
            trigger_error(ArgumentTypeException, ERROR);
        }
    }
679722
public Matrix solve(final Matrix rhs)
    {
        if ((m_rows != m_columns) || (rhs.m_rows != m_columns) || (rhs.m_columns != 1))
        {
            throw new GeometryException("Illegal matrix dimensions");
        }
        // create copies of the data

        final Matrix A = new Matrix(this);

        final Matrix b = new Matrix(rhs);

        // Gaussian elimination with partial pivoting

        for (int i = 0; i < m_columns; i++)
        {
            // find pivot row and swap

            int max = i;

            for (int j = i + 1; j < m_columns; j++)
            {
                if (Math.abs(A.m_data[j][i]) > Math.abs(A.m_data[max][i]))
                {
                    max = j;
                }
            }
            A.swap(i, max);

            b.swap(i, max);

            // singular

            if (A.m_data[i][i] == 0.0)
            {
                throw new RuntimeException("Matrix is singular.");
            }
            // pivot within b

            for (int j = i + 1; j < m_columns; j++)
            {
                b.m_data[j][0] -= (b.m_data[i][0] * A.m_data[j][i]) / A.m_data[i][i];
            }
            // pivot within A

            for (int j = i + 1; j < m_columns; j++)
            {
                final double m = A.m_data[j][i] / A.m_data[i][i];

                for (int k = i + 1; k < m_columns; k++)
                {
                    A.m_data[j][k] -= A.m_data[i][k] * m;
                }
                A.m_data[j][i] = 0.0;
            }
        }
        // back substitution

        final Matrix x = new Matrix(m_columns, 1);

        for (int j = m_columns - 1; j >= 0; j--)
        {
            double t = 0.0;

            for (int k = j + 1; k < m_columns; k++)
            {
                t += A.m_data[j][k] * x.m_data[k][0];
            }
            x.m_data[j][0] = (b.m_data[j][0] - t) / A.m_data[j][j];
        }
        return x;
    }
766185
public function comp($a, $b)
    {
        return bccomp($this->bcNormalize($a), $this->bcNormalize($b));
    }