Solves for X using the QR decomposition.

@param B A matrix that is n by m.  Not modified.
@param X An n by m matrix where the solution is written to.  Modified.
486507
@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;

        Y.reshape(numRows,1, false);
        Z.reshape(numRows,1, false);

        // 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++ ) {
                Y.data[i] = B.get(i,colB);
            }

            // Solve Qa=b
            // a = Q'b
            CommonOps_DDRM.multTransA(Q,Y,Z);

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

            // save the results
            for( int i = 0; i < numCols; i++ ) {
                X.set(i,colB,Z.data[i]);
            }
        }
    }
784154
@Override
	public DoubleMatrix1D[] solve() throws Exception {

		final CholeskyFactorization HFact = new CholeskyFactorization(H, new Matrix1NornRescaler());
		boolean isHReducible = true;
		try{
			HFact.factorize();
		}catch(Exception e){
			isHReducible = false;
		}


        final DoubleMatrix1D v;// dim equals cols of A
        final DoubleMatrix1D w;// dim equals rank of A

		if (isHReducible) {
			// Solving KKT system via elimination
			DoubleMatrix1D HInvg = HFact.solve(g);
			if (A != null) {
				DoubleMatrix2D HInvAT = HFact.solve(AT);
				DoubleMatrix2D MenoSLower = ColtUtils.subdiagonalMultiply(A, HInvAT);
				DoubleMatrix1D AHInvg = ALG.mult(A, HInvg);
				
				CholeskyFactorization MSFact = new CholeskyFactorization(MenoSLower, new Matrix1NornRescaler());
				MSFact.factorize();
				if(h == null){
					w = MSFact.solve(ColtUtils.scalarMult(AHInvg, -1));
				}else{
					w = MSFact.solve(ColtUtils.add(h, AHInvg, -1));
				}
				
				v = HInvg.assign(ALG.mult(HInvAT, w), Functions.plus).assign(Mult.mult(-1));
			} else {
				w = null;
				v = HInvg.assign(Mult.mult(-1));
			}
		} else {
			// H is singular
			// Solving the full KKT system
			if(A!=null){
				DoubleMatrix1D[] fullSol =  this.solveAugmentedKKT();
				v = fullSol[0];
				w = fullSol[1];
			}else{
				//@TODO: try with rescaled H
				throw new Exception("KKT solution failed");
			}
		}

		// 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;
	}
926648
public static FloatMatrix solve(FloatMatrix A, FloatMatrix B) {
		A.assertSquare();
		FloatMatrix X = B.dup();
		int[] ipiv = new int[B.rows];
		SimpleBlas.gesv(A.dup(), ipiv, X);
		return X;
	}
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;
    }
637004
private void solve(int method) throws NullPointerException {
    // solution exists
    if(solved) {
      return;
    }

    // bring in reduced row echelon form
    if(!reducedRowEchelonForm) {
      reducedRowEchelonForm(method);
    }

    if(!isSolvable(method)) {
      if(LOG.isDebugging()) {
        LOG.debugFine("Equation system is not solvable!");
      }
      return;
    }

    // compute one special solution
    final int cols = coeff[0].length;
    int numbound = 0, numfree = 0;
    int[] boundIndices = new int[cols], freeIndices = new int[cols];
    x_0 = new double[cols];
    outer: for(int i = 0; i < coeff.length; i++) {
      for(int j = i; j < coeff[row[i]].length; j++) {
        if(coeff[row[i]][col[j]] == 1) {
          x_0[col[i]] = rhs[row[i]];
          boundIndices[numbound++] = col[i];
          continue outer;
        }
      }
      freeIndices[numfree++] = i;
    }

    StringBuilder msg = new StringBuilder();
    if(LOG.isDebugging()) {
      msg.append("\nSpecial solution x_0 = [").append(FormatUtil.format(x_0, ",", FormatUtil.NF4)).append(']') //
          .append("\nbound Indices ").append(FormatUtil.format(boundIndices, ",")) //
          .append("\nfree Indices ").append(FormatUtil.format(freeIndices, ","));
    }

    // compute solution space of homogeneous linear equation system
    Arrays.sort(boundIndices, 0, numbound);
    int freeIndex = 0;
    int boundIndex = 0;
    u = new double[cols][numfree];

    for(int j = 0; j < u[0].length; j++) {
      for(int i = 0; i < u.length; i++) {
        if(freeIndex < numfree && i == freeIndices[freeIndex]) {
          u[i][j] = 1;
        }
        else if(boundIndex < numbound && i == boundIndices[boundIndex]) {
          u[i][j] = -coeff[row[boundIndex]][freeIndices[freeIndex]];
          boundIndex++;
        }
      }
      freeIndex++;
      boundIndex = 0; // Restart
    }

    if(LOG.isDebugging()) {
      msg.append("\nU");
      for(double[] anU : u) {
        msg.append('\n').append(FormatUtil.format(anU, ",", FormatUtil.NF4));
      }
      LOG.debugFine(msg.toString());
    }

    solved = true;
  }
1368873
public void setMatrix(int i0, int i1, int j0, int j1, Matrix X)
    {
        try
        {
            for (int i = i0; i <= i1; i++)
            {
                for (int j = j0; j <= j1; j++)
                {
                    A[i][j] = X.get(i - i0, j - j0);
                }
            }
        }
        catch (ArrayIndexOutOfBoundsException e)
        {
            throw new ArrayIndexOutOfBoundsException("Submatrix indices");
        }
    }
25465
@Override
    public void solve(ZMatrixRMaj B, ZMatrixRMaj X) {
        if( X.numRows != numCols )
            throw new IllegalArgumentException("Unexpected dimensions for X");
        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[i*2]   = B.data[indexB];
                a[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^H)*b = b - u*(gamma*U^H*b)
            for( int n = 0; n < numCols; n++ ) {
                u[n*2] = 1;
                u[n*2+1] = 0;

                double realUb = a[2*n];
                double imagUb = a[2*n+1];
                // U^H*b
                for( int i = n+1; i < numRows; i++ ) {
                    int indexQR = (i*QR.numCols+n)*2;
                    double realU = u[i*2]   = QR.data[indexQR];
                    double imagU = u[i*2+1] = QR.data[indexQR+1];

                    double realB = a[i*2];
                    double imagB = a[i*2+1];

                    realUb += realU*realB + imagU*imagB;
                    imagUb += realU*imagB - imagU*realB;
                }

                // gamma*U^H*b
                realUb *= gammas[n];
                imagUb *= gammas[n];

                for( int i = n; i < numRows; i++ ) {
                    double realU = u[i*2];
                    double imagU = u[i*2+1];

                    a[i*2  ] -= realU*realUb - imagU*imagUb;
                    a[i*2+1] -= realU*imagUb + imagU*realUb;
                }
            }

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

            // save the results
            for( int i = 0; i < numCols; i++ ) {
                int indexX = (i*X.numCols+colB)*2;

                X.data[indexX]   = a[i*2];
                X.data[indexX+1] = a[i*2+1];
            }
        }
    }
792246
protected DMatrixRMaj solveConstraintMatrix() {

		int rowAA = 0;
		for( int i = 0; i < numControl; i++ ) {
			for( int j = i+1; j < numControl; j++ ) {
				for( int k = j; k < numControl; k++ , rowAA++ ) {
					// x_{ii}*x_{jk} = x_{ik}*x_{ji}
					extractXaXb(getIndex(i, i), getIndex(j, k), XiiXjk);
					extractXaXb(getIndex(i, k), getIndex(j, i), XikXji);

					for( int l = 1; l <= AA.numCols; l++ ) {
						AA.set(rowAA,l-1,XikXji[l]-XiiXjk[l]);
					}
					yy.set(rowAA,XiiXjk[0]-XikXji[0]);
				}
			}
		}
//		AA.print();
		CommonOps_DDRM.solve(AA, yy, xx);

		return xx;
	}
928518
@Nonnull
  @ReturnsMutableCopy
  public Matrix solve (@Nonnull final Matrix aMatrix)
  {
    if (aMatrix.getRowDimension () != m_nRows)
      throw new IllegalArgumentException ("Matrix row dimensions must agree.");
    if (!isNonsingular ())
      throw new IllegalStateException ("Matrix is singular.");

    // Copy right hand side with pivoting
    final int nCols = aMatrix.getColumnDimension ();
    final Matrix aNewMatrix = aMatrix.getMatrix (m_aPivot, 0, nCols - 1);
    final double [] [] aNewArray = aNewMatrix.internalGetArray ();

    // Solve L*Y = B(piv,:)
    for (int k = 0; k < m_nCols; k++)
    {
      final double [] aNewk = aNewArray[k];
      for (int i = k + 1; i < m_nCols; i++)
      {
        final double [] aLUi = m_aLU[i];
        final double [] aNewi = aNewArray[i];
        for (int j = 0; j < nCols; j++)
          aNewi[j] -= aNewk[j] * aLUi[k];
      }
    }
    // Solve U*X = Y;
    for (int k = m_nCols - 1; k >= 0; k--)
    {
      final double [] aLUk = m_aLU[k];
      final double [] aNewk = aNewArray[k];
      for (int j = 0; j < nCols; j++)
        aNewk[j] /= aLUk[k];
      for (int i = 0; i < k; i++)
      {
        final double [] aLUi = m_aLU[i];
        final double [] aNewi = aNewArray[i];
        for (int j = 0; j < nCols; j++)
          aNewi[j] -= aNewk[j] * aLUi[k];
      }
    }
    return aNewMatrix;
  }
1318258
public DoubleMatrix solve(DoubleMatrix b)
    {
        DoubleMatrix x = getInstance(b.rows(), b.columns());
        solve(b, x);
        return x;
    }
1253729
public boolean eq(Matrix B) {
        Matrix A = this;
        if (B.M != A.M || B.N != A.N) {
            throw new RuntimeException("Illegal matrix dimensions.");
        }
        for (int i = 0; i < M; i++) {
            for (int j = 0; j < N; j++) {
                if (A.data[i][j] != B.data[i][j]) {
                    return false;
                }
            }
        }
        return true;
    }
1316486
public void decompose( SimpleMatrix A ) {

        this.QR = A.copy();

        int N = Math.min(A.numCols(),A.numRows());
        gammas = new double[ A.numCols() ];

        for( int i = 0; i < N; i++ ) {
            // use extract matrix to get the column that is to be zeroed
            SimpleMatrix v = QR.extractMatrix(i, END,i,i+1);
            double max = v.elementMaxAbs();

            if( max > 0 && v.getNumElements() > 1 ) {
                // normalize to reduce overflow issues
                v = v.divide(max);

                // compute the magnitude of the vector
                double tau = v.normF();

                if( v.get(0) < 0 )
                    tau *= -1.0;

                double u_0 = v.get(0) + tau;
                double gamma = u_0/tau;

                v = v.divide(u_0);
                v.set(0,1.0);

                // extract the submatrix of A which is being operated on
                SimpleMatrix A_small = QR.extractMatrix(i,END,i,END);

                // A = (I - &gamma;*u*u<sup>T</sup>)A
                A_small = A_small.plus(-gamma,v.mult(v.transpose()).mult(A_small));

                // save the results
                QR.insertIntoThis(i,i,A_small);
                QR.insertIntoThis(i+1,i,v.extractMatrix(1,END,0,1));

                // Alternatively, the two lines above can be replaced with in-place equations
                // READ THE JAVADOC TO UNDERSTAND HOW THIS WORKS!
//                QR.equation("QR(i:,i:) = A","QR",i,"i",A_small,"A");
//                QR.equation("QR((i+1):,i) = v(1:,0)","QR",i,"i",v,"v");

                // save gamma for recomputing Q later on
                gammas[i] = gamma;
            }
        }
    }
1331224
public static double[] solve(double[][] A, double[] b) {
    final int rows = A.length, cols = A[0].length;
    return rows == cols //
        ? (new LUDecomposition(A, rows, cols)).solve(b) //
        : (new QRDecomposition(A, rows, cols)).solve(b);
  }
1324521
public Vec solve(Vec b)
    {
        //Solve  A x = L L^T x = b, for x 
        
        //First solve L y = b
        Vec y = forwardSub(L, b);
        //Sole L^T x = y
        Vec x = backSub(L, y);
        
        return x;
    }
1232964
public void _solveVectorInternal( double []vv )
    {
        // Solve L*Y = B
        int ii = 0;

        for( int i = 0; i < n; i++ ) {
            int ip = indx[i];
            double sum = vv[ip];
            vv[ip] = vv[i];
            if( ii != 0 ) {
//                for( int j = ii-1; j < i; j++ )
//                    sum -= dataLU[i* n +j]*vv[j];
                int index = i*n + ii-1;
                for( int j = ii-1; j < i; j++ )
                    sum -= dataLU[index++]*vv[j];
            } else if( sum != 0.0 ) {
                ii=i+1;
            }
            vv[i] = sum;
        }

        // Solve U*X = Y;
        TriangularSolver_DDRM.solveU(dataLU,vv,n);
    }
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);
        }
    }
939018
public static void subtractEquals(DMatrixD1 a, DMatrixD1 b)
    {
        if( a.numCols != b.numCols || a.numRows != b.numRows ) {
            throw new MatrixDimensionException("The 'a' and 'b' matrices do not have compatible dimensions");
        }

        final int length = a.getNumElements();

        for( int i = 0; i < length; i++ ) {
            a.data[i] -= b.data[i];
        }
    }