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.
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];
            }
        }
    }
358801
private void extractSolutionForQ( DMatrix4x4 Q ) {
		DMatrixRMaj nv = new DMatrixRMaj(10,1);
		SingularOps_DDRM.nullVector(svd,true,nv);

		// Convert the solution into a fixed sized matrix because it's easier to read
		encodeQ(Q,nv.data);

		// diagonal elements must be positive because Q = [K*K' .. ; ... ]
		// If they are a mix of positive and negative that's bad and can't be fixed
		// All 3 are checked because technically they could be zero. Not a physically useful solution but...
		if( Q.a11 < 0 || Q.a22 < 0 || Q.a33 < 0 ){
			CommonOps_DDF4.scale(-1,Q);
		}
	}
134286
public void setMatrix(int i0, int i1, int[] c, Matrix X)
    {
        try
        {
            for (int i = i0; i <= i1; i++)
            {
                for (int j = 0; j < c.length; j++)
                {
                    A[i][c[j]] = X.get(i - i0, j);
                }
            }
        }
        catch (ArrayIndexOutOfBoundsException e)
        {
            throw new ArrayIndexOutOfBoundsException("Submatrix indices");
        }
    }
394707
public static DoubleMatrix absi(DoubleMatrix x) { 
		/*# mapfct('Math.abs') #*/
//RJPP-BEGIN------------------------------------------------------------
	   for (int i = 0; i < x.length; i++)
	      x.put(i, (double) Math.abs(x.get(i)));
	   return x;
//RJPP-END--------------------------------------------------------------
	}
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;
  }
306101
public static void addEquals( DMatrix5x5 a , DMatrix5x5 b ) {
        a.a11 += b.a11;
        a.a12 += b.a12;
        a.a13 += b.a13;
        a.a14 += b.a14;
        a.a15 += b.a15;
        a.a21 += b.a21;
        a.a22 += b.a22;
        a.a23 += b.a23;
        a.a24 += b.a24;
        a.a25 += b.a25;
        a.a31 += b.a31;
        a.a32 += b.a32;
        a.a33 += b.a33;
        a.a34 += b.a34;
        a.a35 += b.a35;
        a.a41 += b.a41;
        a.a42 += b.a42;
        a.a43 += b.a43;
        a.a44 += b.a44;
        a.a45 += b.a45;
        a.a51 += b.a51;
        a.a52 += b.a52;
        a.a53 += b.a53;
        a.a54 += b.a54;
        a.a55 += b.a55;
    }
1165685
private void solveInternalL() {
        // This takes advantage of the diagonal elements always being real numbers

        // solve L*y=b storing y in x
        TriangularSolver_ZDRM.solveL_diagReal(t, vv, n);

        // solve L^T*x=y
        TriangularSolver_ZDRM.solveConjTranL_diagReal(t, vv, n);
    }
1331107
@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(H.toArray()));
			log.debug("g: " + ArrayUtils.toString(g.toArray()));
			if (A != null) {
				log.debug("A: " + ArrayUtils.toString(A.toArray()));
			}
			if (h != null) {
				log.debug("h: " + ArrayUtils.toString(h.toArray()));
			}
		}
		
		//compose the full KKT matrix
		DoubleMatrix2D KKT = null;
		DoubleMatrix1D b = null;
		final DoubleMatrix2D HFull = ColtUtils.fillSubdiagonalSymmetricMatrix((SparseDoubleMatrix2D)this.H);
		
		if (this.A != null) {
			if(h!=null){
				//H.v + [A]T.w = -g
				//A.v = -h
				DoubleMatrix2D[][] parts = {
						{ HFull, this.AT },
						{ this.A, null } };
				if(HFull instanceof SparseDoubleMatrix2D && A instanceof SparseDoubleMatrix2D){
					KKT = DoubleFactory2D.sparse.compose(parts);
				}else{
					KKT = DoubleFactory2D.dense.compose(parts);
				}
				
				b = F1.append(g, h).assign(Mult.mult(-1));
			}else{
				//H.v + [A]T.w = -g
				DoubleMatrix2D[][] parts = {{ HFull, this.AT }};
				if(HFull instanceof SparseDoubleMatrix2D && A instanceof SparseDoubleMatrix2D){
					KKT = DoubleFactory2D.sparse.compose(parts);
				}else{
					KKT = DoubleFactory2D.dense.compose(parts);
				}
				b = ColtUtils.scalarMult(g, -1);
			}
		}else{
			KKT = HFull;
			b = ColtUtils.scalarMult(g, -1);
		}
		
		//factorization
		//LDLTPermutedFactorization KKTFact = new LDLTPermutedFactorization(KKT, new RequestFilter());
		QRSparseFactorization KKTFact = new QRSparseFactorization((SparseDoubleMatrix2D)KKT, new Matrix1NornRescaler());
		try{
			KKTFact.factorize();
		}catch(Exception e){
			log.error("singular KKT system");
			throw new Exception("singular KKT system");
		}
		
		DoubleMatrix1D x = KKTFact.solve(b);
		v = x.viewPart(0, H.rows());
		w = x.viewPart(H.rows()-1, x.size()-H.rows());

		// 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;// dim equals cols of A
		ret[1] = w;
		return ret;
	}
1269159
public double[] solveLInplace(double[] X) {
    final int n = L.length;
    X[0] /= L[0][0]; // First iteration, simplified.
    for(int k = 1; k < n; k++) {
      final double[] Lk = L[k];
      for(int i = 0; i < k; i++) {
        X[k] -= X[i] * Lk[i];
      }
      X[k] /= Lk[k];
    }
    return X;
  }
81780
@Override
    public void solve(DMatrixRBlock B, DMatrixRBlock X) {
        if( B.blockLength != blockLength )
            throw new IllegalArgumentException("Unexpected blocklength in B.");

        DSubmatrixD1 L = new DSubmatrixD1(decomposer.getT(null));

        if( X != null ) {
            if( X.blockLength != blockLength )
                throw new IllegalArgumentException("Unexpected blocklength in X.");
            if( X.numRows != L.col1 ) throw new IllegalArgumentException("Not enough rows in X");
        }
        
        if( B.numRows != L.col1 ) throw new IllegalArgumentException("Not enough rows in B");

        //  L * L^T*X = B

        // Solve for Y:  L*Y = B
        TriangularSolver_DDRB.solve(blockLength,false,L,new DSubmatrixD1(B),false);

        // L^T * X = Y
        TriangularSolver_DDRB.solve(blockLength,false,L,new DSubmatrixD1(B),true);

        if( X != null ) {
            // copy the solution from B into X
            MatrixOps_DDRB.extractAligned(B,X);
        }

    }
852542
public void solve(double[] b, double[] x) {
        if (b.length != qr.nrows()) {
            throw new IllegalArgumentException(String.format("Row dimensions do not agree: A is %d x %d, but B is %d x 1", qr.nrows(), qr.nrows(), b.length));
        }

        if (x.length != qr.ncols()) {
            throw new IllegalArgumentException("A and x dimensions don't match.");
        }

        if (singular) {
            throw new RuntimeException("Matrix is rank deficient.");
        }

        double[] B = b.clone();
        solve(Matrix.newInstance(B));
        System.arraycopy(B, 0, x, 0, x.length);
    }
947061
public static void solveL(DMatrixSparseCSC L , double []x )
    {
        final int N = L.numCols;

        int idx0 = L.col_idx[0];
        for (int col = 0; col < N; col++) {

            int idx1 = L.col_idx[col+1];
            double x_j = x[col] /= L.nz_values[idx0];

            for (int i = idx0+1; i < idx1; i++) {
                int row = L.nz_rows[i];
                x[row] -=  L.nz_values[i]*x_j;
            }

            idx0 = idx1;
        }
    }
1289696
@Override
	public DoubleMatrix1D[] solve() throws Exception {

		if (A == null) {
			throw new IllegalStateException("Matrix A cannot be null");
		}

		DoubleMatrix1D v = null;// dim equals cols of A
		DoubleMatrix1D w = null;// dim equals rank of A

		if (log.isDebugEnabled()) {
			log.debug("H: " + ArrayUtils.toString(H.toArray()));
			log.debug("g: " + ArrayUtils.toString(g.toArray()));
			log.debug("A: " + ArrayUtils.toString(A.toArray()));
			if (h != null) {
				log.debug("h: " + ArrayUtils.toString(h.toArray()));
			}
		}

		// augmentation
		final DoubleMatrix2D HAugm = ColtUtils.subdiagonalMultiply(AT, A);// H + ATQA
		HAugm.forEachNonZero(new IntIntDoubleFunction() {
			@Override
			public double apply(int i, int j, double HAugmij) {
				return s * HAugmij;
			}
		});
		H.forEachNonZero(new IntIntDoubleFunction() {
			@Override
			public double apply(int i, int j, double Hij) {
				if (i + 1 > j) {
					// the subdiagonal elements
					HAugm.setQuick(i, j, Hij + HAugm.getQuick(i, j));
				}
				return Hij;
			}
		});

		DoubleMatrix1D gAugm = null;// g + ATQh
		if (h != null) {
			DoubleMatrix1D ATQh = ALG.mult(AT, ColtUtils.diagonalMatrixMult(F1.make(A.rows(), 1), h));
			DoubleMatrix1D gATQh = ColtUtils.add(g, ATQh, defaultScalar);
			gAugm = gATQh;
		}else{
			gAugm = g.copy();
		}

		// solving the augmented system
		CholeskyFactorization HFact = new CholeskyFactorization(HAugm, (MatrixRescaler) new Matrix1NornRescaler());
		try {
			HFact.factorize();
		} catch (Exception e) {
			// The KKT matrix is nonsingular if and only if H + ATQA > 0 for some Q > 0
			log.error("singular KKT system");
			throw new Exception("singular KKT system");
		}

		// Solving KKT system via elimination
		DoubleMatrix1D HInvg = HFact.solve(gAugm);
		DoubleMatrix2D HInvAT = HFact.solve(AT);
		DoubleMatrix2D MenoSLower = ColtUtils.subdiagonalMultiply(A, HInvAT);
		DoubleMatrix1D AHInvg = ALG.mult(A, HInvg);

		CholeskyFactorization MSFact = new CholeskyFactorization(MenoSLower, (MatrixRescaler) 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));

		// 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;
	}
560
def _initialize_b_x(self):
        
        """
        # Get vector(s) compatible with the matrix,
        # i.e., with the same parallel layout.
        self.petsc_x, self.petsc_b = self.petsc_A.getVecs()

        #  Set the solution vector to zeros.
        self.petsc_x.set(0)

        # Define the petsc rhs vector from the numpy one.
        # If the rhs is defined by blocks, use this:
        PETSc.Vec.setValuesBlocked(self.petsc_b, [sp.arange(self.m)], self.b)
651471
public Matrix solve (Matrix B) {
      return (m == n ? (new LUDecomposition(this)).solve(B) :
                       (new QRDecomposition(this)).solve(B));
   }
899971
protected void checkSizes(Matrix A, Vector b, Vector x) {
        if (!A.isSquare())
            throw new IllegalArgumentException("!A.isSquare()");
        if (b.size() != A.numRows())
            throw new IllegalArgumentException("b.size() != A.numRows()");
        if (b.size() != x.size())
            throw new IllegalArgumentException("b.size() != x.size()");
    }
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;
    }