Iteratively improve the solution x to machine accuracy.
@param b the right-hand side column vector
@param x the improved solution column vector
@throws matrix.MatrixException if failed to converge
349867
private void improve(ColumnVector b, ColumnVector x)
        throws MatrixException
    {
        // Find the largest x element.
        double largestX = 0;
        for (int r = 0; r < nRows; ++r) {
            double absX = Math.abs(x.values[r][0]);
            if (largestX < absX) largestX = absX;
        }

        // Is x already as good as possible?
        if (largestX == 0) return;

        ColumnVector residuals = new ColumnVector(nRows);

        // Iterate to improve x.
        for (int iter = 0; iter < MAX_ITER; ++iter) {

            // Compute residuals = b - Ax.
            // Must use double precision!
            for (int r = 0; r < nRows; ++r) {
                double dot   = 0;
                double  row[] = values[r];
                for (int c = 0; c < nRows; ++c) {
                    double elmt = at(r, c);
                    dot += elmt*x.at(c);        // dbl.prec. *
                }
                double value = b.at(r) - dot;   // dbl.prec. -
                residuals.set(r, (double) value);
            }

            // Solve Az = residuals for z.
            ColumnVector z = solve(residuals, false);

            // Set x = x + z.
            // Find largest the largest difference.
            double largestDiff = 0;
            for (int r = 0; r < nRows; ++r) {
                double oldX = x.at(r);
                x.set(r, oldX + z.at(r));

                double diff = Math.abs(x.at(r) - oldX);
                if (largestDiff < diff) largestDiff = diff;
            }

            // Is any further improvement possible?
            if (largestDiff < largestX*TOLERANCE) return;
        }

        // Failed to converge because A is nearly singular.
        throw new MatrixException(MatrixException.NO_CONVERGENCE);
    }
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];
            }
        }
    }
851509
public static void elementPower(double a , DMatrixD1 B , DMatrixD1 C ) {

        if( B.numRows != C.numRows || B.numCols != C.numCols ) {
            throw new MatrixDimensionException("All matrices must be the same shape");
        }

        int size = B.getNumElements();
        for( int i = 0; i < size; i++ ) {
            C.data[i] = Math.pow(a, B.data[i]);
        }
    }
1241069
public Matrix subtract(Matrix B, ExecutorService threadPool)
    {
        Matrix toReturn = getThisSideMatrix(B);
        toReturn.mutableSubtract(1.0, B, threadPool);
        return toReturn;
    }
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];
            }
        }
    }
927880
public void solve(DenseMatrix B) {
        if (B.nrows() != L.nrows()) {
            throw new IllegalArgumentException(String.format("Row dimensions do not agree: A is %d x %d, but B is %d x %d", L.nrows(), L.ncols(), B.nrows(), B.ncols()));
        }

        int n = B.nrows();
        int nrhs = B.ncols();

        // Solve L*Y = B;
        for (int k = 0; k < n; k++) {
            for (int j = 0; j < nrhs; j++) {
                for (int i = 0; i < k; i++) {
                    B.sub(k, j, B.get(i, j) * L.get(k, i));
                }
                B.div(k, j, L.get(k, k));
            }
        }

        // Solve L'*X = Y;
        for (int k = n - 1; k >= 0; k--) {
            for (int j = 0; j < nrhs; j++) {
                for (int i = k + 1; i < n; i++) {
                    B.sub(k, j, B.get(i, j) * L.get(i, k));
                }
                B.div(k, j, L.get(k, k));
            }
        }
    }
413676
public Matrix add(Matrix B, ExecutorService threadPool)
    {
        Matrix toReturn = getThisSideMatrix(B);
        toReturn.mutableAdd(1.0, B, threadPool);
        return toReturn;
    }
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;
	}
184818
public static void addEquals( DMatrix2x2 a , DMatrix2x2 b ) {
        a.a11 += b.a11;
        a.a12 += b.a12;
        a.a21 += b.a21;
        a.a22 += b.a22;
    }
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];
        }
    }
590725
public Matrix solve (Matrix B) {
      if (B.getRowDimension() != m) {
         throw new IllegalArgumentException("Matrix row dimensions must agree.");
      }
      if (!this.isNonsingular()) {
         throw new RuntimeException("Matrix is singular.");
      }

      // Copy right hand side with pivoting
      int nx = B.getColumnDimension();
      Matrix Xmat = B.getMatrix(piv,0,nx-1);
      double[][] X = Xmat.getArray();

      // Solve L*Y = B(piv,:)
      for (int k = 0; k < n; k++) {
         for (int i = k+1; i < n; i++) {
            for (int j = 0; j < nx; j++) {
               X[i][j] -= X[k][j]*LU[i][k];
            }
         }
      }
      // Solve U*X = Y;
      for (int k = n-1; k >= 0; k--) {
         for (int j = 0; j < nx; j++) {
            X[k][j] /= LU[k][k];
         }
         for (int i = 0; i < k; i++) {
            for (int j = 0; j < nx; j++) {
               X[i][j] -= X[k][j]*LU[i][k];
            }
         }
      }
      return Xmat;
   }
1287983
public double[][] solve(double[][] B) {
        if (B.length != m) {
            throw new IllegalArgumentException("Matrix row dimensions must agree.");
        }
        if (!this.isNonsingular()) {
            throw new RuntimeException("Matrix is singular.");
        }

        // Copy right hand side with pivoting
        int nx = B[0].length;
        double[][] X = Matrix.Submatrix(B, piv, 0, nx - 1);

        // Solve L*Y = B(piv,:)
        for (int k = 0; k < n; k++) {
            for (int i = k + 1; i < n; i++) {
                for (int j = 0; j < nx; j++) {
                    X[i][j] -= X[k][j] * LU[i][k];
                }
            }
        }
        // Solve U*X = Y;
        for (int k = n - 1; k >= 0; k--) {
            for (int j = 0; j < nx; j++) {
                X[k][j] /= LU[k][k];
            }
            for (int i = 0; i < k; i++) {
                for (int j = 0; j < nx; j++) {
                    X[i][j] -= X[k][j] * LU[i][k];
                }
            }
        }
        return X;
    }
956647
@Override
    public Matrix applyCommon(Matrix a, Matrix b) {
        int n = a.rows() * b.rows();
        int m = a.columns() * b.columns();
        int p = b.rows();
        int q = b.columns();

        Matrix result = a.blankOfShape(n, m);

        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                result.set(i, j, a.get(i / p, j / q) * b.get(i % p, j % q));
            }
        }

        return result;
    }
362156
protected void upgradeSolution( DMatrixRMaj X ) {
        DMatrixRMaj z = Y; // recycle Y

        // compute the z which will minimize the 2-norm of X
        // because of the identity matrix tacked onto the end 'A' should never be singular
        if( !internalSolver.setA(W) )
            throw new RuntimeException("This should never happen.  Is input NaN?");
        z.reshape(numCols-rank,1);
        internalSolver.solve(X, z);

        // compute X by tweaking the original
        CommonOps_DDRM.multAdd(-1, W, z, X);
    }
87239
public static void elementPower(DMatrixD1 A , double b, DMatrixD1 C ) {

        if( A.numRows != C.numRows || A.numCols != C.numCols ) {
            throw new MatrixDimensionException("All matrices must be the same shape");
        }

        int size = A.getNumElements();
        for( int i = 0; i < size; i++ ) {
            C.data[i] = Math.pow(A.data[i], b);
        }
    }
343008
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());
		}

		throw new RuntimeException("not yet implemented");
	}