Solve A*X = B
@param  $B  A Matrix with as many rows as A and any number of columns.
@return  X so that L*U*X = B(piv,:)
@exception  IllegalArgumentException Matrix row dimensions must agree.
@exception  RuntimeException  Matrix is singular.
970521
public function solve($B)
    {
        if ($B->getRowDimension() == $this->m) {
            if ($this->isNonsingular()) {
                // Copy right hand side with pivoting
                $nx = $B->getColumnDimension();
                $X  = $B->getMatrix($this->piv, 0, $nx-1);
                // Solve L*Y = B(piv,:)
                for ($k = 0; $k < $this->n; $k++) {
                    for ($i = $k+1; $i < $this->n; $i++) {
                        for ($j = 0; $j < $nx; $j++) {
                            $X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k];
                        }
                    }
                }
                // Solve U*X = Y;
                for ($k = $this->n-1; $k >= 0; $k--) {
                    for ($j = 0; $j < $nx; $j++) {
                        $X->A[$k][$j] /= $this->LU[$k][$k];
                    }
                    for ($i = 0; $i < $k; $i++) {
                        for ($j = 0; $j < $nx; $j++) {
                            $X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k];
                        }
                    }
                }
                return $X;
            } else {
                trigger_error(MatrixSingularException, ERROR);
            }
        } else {
            trigger_error(MatrixSquareException, ERROR);
        }
    }
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);
        }
    }
1123995
public void solve(DenseMatrix B) {
        if (!full) {
            throw new IllegalStateException("This is not a FULL singular value decomposition.");
        }

        if (B.nrows() != m) {
            throw new IllegalArgumentException("Dimensions do not agree.");
        }

        double[] b = new double[m];
        double[] x = new double[n];
        int p = B.ncols();
        for (int j = 0; j < p; j++) {
            for (int i = 0; i < m; i++) {
                b[i] = B.get(i, j);
            }

            solve(b, x);
            for (int i = 0; i < n; i++) {
                B.set(i, j, x[i]);
            }
        }
    }
1318258
public DoubleMatrix solve(DoubleMatrix b)
    {
        DoubleMatrix x = getInstance(b.rows(), b.columns());
        solve(b, x);
        return x;
    }
319942
public IMatrix mul(IMatrix b) {
        IMatrix result = new IMatrix(rows, columns);
        mul(b, result);
        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];
            }
        }
    }
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;
    }
467630
public double[] solveInplace(double[] b) {
    if(b.length != m) {
      throw new IllegalArgumentException(ERR_MATRIX_DIMENSIONS);
    }
    if(!this.isNonsingular()) {
      throw new ArithmeticException(ERR_SINGULAR);
    }
    // Solve L*Y = B(piv,:)
    for(int k = 0; k < n; k++) {
      for(int i = k + 1; i < n; i++) {
        b[i] -= b[k] * LU[i][k];
      }
    }
    // Solve U*X = Y;
    for(int k = n - 1; k >= 0; k--) {
      b[k] /= LU[k][k];
      for(int i = 0; i < k; i++) {
        b[i] -= b[k] * LU[i][k];
      }
    }
    return b;
  }
759415
@Override
	public DoubleMatrix1D[] solve() throws Exception {

		DoubleMatrix1D v = null;// dim equals cols of A
		DoubleMatrix1D w = null;// dim equals rank of A
		
    if (log.isDebugEnabled()) {
			log.debug("H: " + ArrayUtils.toString(ColtUtils.fillSubdiagonalSymmetricMatrix((SparseDoubleMatrix2D)this.H).toArray()));
			log.debug("g: " + ArrayUtils.toString(g.toArray()));
			if (A!=null) {
				log.debug("A: " + ArrayUtils.toString(A.toArray()));
				log.debug("h: " + ArrayUtils.toString(h.toArray()));
			}
		}
    	MatrixRescaler rescaler = (avoidScaling)? null : new Matrix1NornRescaler();
		CholeskyUpperDiagonalFactorization HFact = new CholeskyUpperDiagonalFactorization((SparseDoubleMatrix2D)H, diagonalLength, rescaler);
		boolean isHReducible = true;
		try{
			HFact.factorize();
		}catch(Exception e){
			log.warn("warn", e);
			isHReducible = false;
		}

		if (isHReducible) {
			// Solving KKT system via elimination
			DoubleMatrix1D HInvg;
			HInvg = HFact.solve(g);
			
			if (A != null) {
				DoubleMatrix2D HInvAT;
				HInvAT = HFact.solve(AT);
				
//				double scaledResidualX = Utils.calculateScaledResidual(ColtUtils.fillSubdiagonalSymmetricMatrix(H), HInvAT, AT);
//				log.info("scaledResidualX: " + scaledResidualX);
				
				DoubleMatrix2D MenoSLower = ColtUtils.subdiagonalMultiply(A, HInvAT);
				log.debug("MenoS: " + ArrayUtils.toString(ColtUtils.fillSubdiagonalSymmetricMatrix(MenoSLower).toArray()));
				DoubleMatrix1D AHInvg = ALG.mult(A, HInvg);
				
				CholeskyFactorization MSFact = new CholeskyFactorization(MenoSLower, new Matrix1NornRescaler());
				//CholeskySparseFactorization MSFact = new CholeskySparseFactorization((SparseDoubleMatrix2D)MenoSLower, new RequestFilter());
				try{
					MSFact.factorize();
					if(h == null){
						w = MSFact.solve(ColtUtils.scalarMult(AHInvg, -1));
					}else{
						//w = MSFact.solve(h.copy().assign(AHInvg, Functions.minus));
						DoubleMatrix1D hmAHInvg = ColtUtils.add(h, AHInvg, -1);
						//log.debug("hmAHInvg: " + ArrayUtils.toString(hmAHInvg.toArray()));
						w = MSFact.solve(hmAHInvg);
					}
					
					v = HInvg.assign(ALG.mult(HInvAT, w), Functions.plus).assign(Mult.mult(-1));
				}catch(Exception e){
					log.warn("warn: "+ e.getMessage());
					log.debug("MenoS: " + ArrayUtils.toString(ColtUtils.fillSubdiagonalSymmetricMatrix(MenoSLower).toArray()));
					//is it a numeric issue? try solving the full KKT system
					try{
						//NOTE: it would be more appropriate to try solving the full KKT, but if the decomposition 
						//of the Shur complement of H (>0) in KKT fails it is certainty for a numerical issue and
						//the augmented KKT seems to be more able to recover from this situation
						//DoubleMatrix1D[] fullSol =  this.solveFullKKT();
						DoubleMatrix1D[] fullSol =  this.solveAugmentedKKT();
						v = fullSol[0];
						w = fullSol[1];
					}catch(Exception ex){
						log.error(ex.getMessage());
						throw new Exception("KKT solution failed");
					}
				}
				
			} else {
				//A==null
				w = null;
				v = HInvg.assign(Mult.mult(-1));
			}
		} else {
			// H not isReducible, try solving the augmented KKT system
			DoubleMatrix1D[] fullSol = this.solveAugmentedKKT();
			v = fullSol[0];
			w = fullSol[1];
		}

		// solution checking
		if (this.checkKKTSolutionAccuracy && !this.checkKKTSolutionAccuracy(v, w)) {
			log.error("KKT solution failed");
			throw new Exception("KKT solution failed");
		}

		DoubleMatrix1D[] ret = new DoubleMatrix1D[2];
		ret[0] = v;
		ret[1] = w;
		return ret;
	}
695698
public static function solve(\CArray $a, \CArray $b): \CArray
    {
        return parent::solve($a, $b);
    }
93854
public DoubleMatrix2D solve(DoubleMatrix2D B) {
		if (B.rows() != dim) {
			log.error("wrong dimension of vector b: expected " + dim + ", actual "	+ B.rows());
			throw new RuntimeException("wrong dimension of vector b: expected " + dim	+ ", actual " + B.rows());
		}

	// 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);
		}

		int nOfColumns = B.columns();

		// Solve LY = B (the same as L.Yc = Bc for each column Yc e Bc)
		final double[][] Y = new double[dim][nOfColumns];
		for (int i = 0; i < diagonalLength; i++) {
			double LII = LData[0][i];
			double[] YI = Y[i];
			DoubleMatrix1D BI = B.viewRow(i);
			for (int col = 0; col < nOfColumns; col++) {
				YI[col] = BI.getQuick(col) / LII;
			}
		}
		for (int i = diagonalLength; i < dim; i++) {
			double[] LI = LData[i - diagonalLength + 1];
			double LII = LI[i];
			double[] sum = new double[nOfColumns];
			for (int j = 0; j < i; j++) {
				double LIJ = LI[j];
				double[] YJ = Y[j];
				for (int col = 0; col < nOfColumns; col++) {
					sum[col] += LIJ * YJ[col];
				}
			}
			double[] YI = Y[i];
			DoubleMatrix1D BI = B.viewRow(i);
			for (int col = 0; col < nOfColumns; col++) {
				YI[col] = (BI.getQuick(col) - sum[col]) / LII;
			}
		}

		// Solve L[T].X = Y (the same as L[T].Xc = Yc for each column)
		final DoubleMatrix2D X = F2.make(dim, nOfColumns);
		for (int i = dim - 1; i > diagonalLength - 1; i--) {
			double LII = LData[i - diagonalLength + 1][i];
			double[] sum = new double[nOfColumns];
			for (int j = dim - 1; j > i; j--) {
				double[] LJ = LData[j - diagonalLength + 1];
				double LJI = LJ[i];
				DoubleMatrix1D XJ = X.viewRow(j);
				for (int col = 0; col < nOfColumns; col++) {
					sum[col] += LJI * XJ.getQuick(col);
				}
			}
			DoubleMatrix1D XI = X.viewRow(i);
			double[] YI = Y[i];
			for (int col = 0; col < nOfColumns; col++) {
				XI.setQuick(col, (YI[col] - sum[col]) / LII);
			}
		}
		for (int i = diagonalLength - 1; i > -1; i--) {
			double LII = LData[0][i];
			double[] sum = new double[nOfColumns];
			for (int j = dim - 1; j > diagonalLength - 1; j--) {
				double[] LJ = LData[j - diagonalLength + 1];
				double LJI = LJ[i];
				DoubleMatrix1D XJ = X.viewRow(j);
				for (int col = 0; col < nOfColumns; col++) {
					sum[col] += LJI * XJ.getQuick(col);
				}
			}
			DoubleMatrix1D XI = X.viewRow(i);
			double[] YI = Y[i];
			for (int col = 0; col < nOfColumns; col++) {
				XI.setQuick(col, (YI[col] - sum[col]) / LII);
			}
		}

		if (this.rescaler != null) {
			// return ALG.mult(this.U, X);
			return ColtUtils.diagonalMatrixMult(this.U, X);
		} else {
			return X;
		}
	}
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;
    }
679998
public function subtract($b)
    {
        switch (true) {
            case $b instanceof Matrix:
                $c = $this->subtractMatrix($b);
                break;

            case $b instanceof ColumnVector:
                $c = $this->subtractColumnVector($b);
                break;
                
            case $b instanceof Vector:
                $c = $this->subtractVector($b);
                break;

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

            default:
                throw new InvalidArgumentException('Cannot subtract a '
                    . gettype($b) . ' from a matrix.');
        }

        return $c;
    }
60985
public static function calculateDeterminant($matrix) : float {
        // The ´rows´ args is the only needed one. Determinants can only be calculated
        // from a square matrix. Hence we use the row as size for both row and column.
        $determinant = 0;
        $rcCount     = count($matrix);
        if ($rcCount !== count($matrix[0])) {
            throw new InvalidArgumentException(
                "Can only calculate determinant from square matrix (same amount of rows as columns)."
            );
        }

        // If only one element, return the value from it.
        if ($rcCount === 1) {
            return $matrix[0][0];
        }

        $sign = 1;
        // Go through each element in first row.
        for ($counter=0;$counter<$rcCount;$counter++) {
            // Create a cofactor sub-matrix ($rcCount-1 in size).
            $subMatrix = self::getSubMatrix($matrix, 0, $counter);
            // Calculate the determinant of the sub-matrix.
            $determinant += $sign * $matrix[0][$counter] * self::calculateDeterminant($subMatrix);
            $sign         = -$sign;
        }

        return $determinant;
    }
1269832
public function outer(Vector $b) : Matrix
    {
        $n = $b->n();

        $c = [];

        foreach ($this->a as $i => $value) {
            $temp = [];

            for ($j = 0; $j < $n; $j++) {
                $temp[] = $value * $b[$j];
            }

            $c[] = $temp;
        }

        return Matrix::quick($c);
    }
651471
public Matrix solve (Matrix B) {
      return (m == n ? (new LUDecomposition(this)).solve(B) :
                       (new QRDecomposition(this)).solve(B));
   }
111244
public double[] solve(double[] b, double[] x) {
        if (b.length != n) {
            throw new IllegalArgumentException(String.format("Row dimensions do not agree: A is %d x %d, but b is %d x 1", n, n, b.length));
        }

        if (b.length != x.length) {
            throw new IllegalArgumentException("b and x dimensions do not agree.");
        }

        if (au == null) {
            decompose();
        }
        System.arraycopy(b, 0, x, 0, n);

        int mm = m1 + m2 + 1;
        int l = m1;

        for (int k = 0; k < n; k++) {
            int j = index[k] - 1;

            if (j != k) {
                // swap x[k] and x[j]
                double swap = x[k];
                x[k] = x[j];
                x[j] = swap;
            }

            if (l < n) {
                l++;
            }

            for (j = k + 1; j < l; j++) {
                x[j] -= al[k][j - k - 1] * x[k];
            }
        }

        l = 1;
        for (int i = n - 1; i >= 0; i--) {
            double dum = x[i];
            for (int k = 1; k < l; k++) {
                dum -= au[i][k] * x[k + i];
            }
            x[i] = dum / au[i][0];
            if (l < mm) {
                l++;
            }
        }

        return x;
    }