Iteratively improve a solution to linear equations.

@param b right hand side of linear equations.
@param x a solution to linear equations.
546485
public void improve(double[] b, double[] x) {
        if (b.length != n || x.length != n) {
            throw new IllegalArgumentException(String.format("Row dimensions do not agree: A is %d x %d, but b is %d x 1 and x is %d x 1", n, n, b.length, x.length));
        }

        // Calculate the right-hand side, accumulating the residual
        // in higher precision.
        double[] r = b.clone();
        axpy(x, r, -1.0);

        // Solve for the error term.
        solve(r, r);

        // Subtract the error from the old solution.
        for (int i = 0; i < n; i++) {
            x[i] -= r[i];
        }
    }
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];
            }
        }
    }
747566
public Matrix add(Matrix B)
    {
        Matrix toReturn = getThisSideMatrix(B);
        toReturn.mutableAdd(1.0, B);
        return toReturn;
    }
73970
public void sub(IMatrix b, IMatrix result) {
        if ((b == null) || (rows != b.rows) || (columns != b.columns)) return;

        if ((result.rows != rows) || (result.columns != columns)) result.reshape(rows, columns);

        int i, j;
        for (i = 0; i < rows; i++)
            for (j = 0; j < columns; j++) {
                result.realmatrix[i][j] = realmatrix[i][j] - b.realmatrix[i][j];
                result.imagmatrix[i][j] = imagmatrix[i][j] - b.imagmatrix[i][j];
            }
    }
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;
        }
    }
460393
public ColumnVector solve(ColumnVector b, boolean improve)
        throws MatrixException
    {
        // Validate b's size.
        if (b.nRows != nRows) {
            throw new MatrixException(
                                MatrixException.INVALID_DIMENSIONS);
        }

        decompose();

        // Solve Ly = b for y by forward substitution.
        // Solve Ux = y for x by back substitution.
        ColumnVector y = forwardSubstitution(b);
        ColumnVector x = backSubstitution(y);

        // Improve and return x.
        if (improve) improve(b, x);
        return x;
    }
566987
public static void multiplyInPlace(double[] a, double b) {
    for (int i = 0; i < a.length; i++) {
      a[i] = a[i] * b;
    }
  }
802870
public void setMatrix(int[] r, int j0, int j1, Matrix X)
    {
        try
        {
            for (int i = 0; i < r.length; i++)
            {
                for (int j = j0; j <= j1; j++)
                {
                    A[r[i]][j] = X.get(i, j - j0);
                }
            }
        }
        catch (ArrayIndexOutOfBoundsException e)
        {
            throw new ArrayIndexOutOfBoundsException("Submatrix indices");
        }
    }
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()");
    }
119343
public static void addEquals( DMatrix5 a , DMatrix5 b ) {
        a.a1 += b.a1;
        a.a2 += b.a2;
        a.a3 += b.a3;
        a.a4 += b.a4;
        a.a5 += b.a5;
    }
494796
public void add(double x, double y)
    {
        if (count >= xarr.length)
        {
            grow();
        }
        xarr[count] = x;
        yarr[count++] = y;
        minX = Math.min(minX, x);
        minY = Math.min(minY, y);
        maxX = Math.max(maxX, x);
        maxY = Math.max(maxY, y);
    }
124673
public double[][] solve(double[][] B) {

        if (B.length != n) {
            throw new IllegalArgumentException("Matrix row dimensions must agree.");
        }
        if (!isspd) {
            throw new RuntimeException("Matrix is not symmetric positive definite.");
        }

        // Copy right hand side.
        double[][] X = Matrix.Copy(B);
        int nx = B[0].length;

        // Solve L*Y = B;
        for (int k = 0; k < n; k++) {
            for (int j = 0; j < nx; j++) {
                for (int i = 0; i < k; i++) {
                    X[k][j] -= X[i][j] * L[k][i];
                }
                X[k][j] /= L[k][k];
            }
        }

        // Solve L'*X = Y;
        for (int k = n - 1; k >= 0; k--) {
            for (int j = 0; j < nx; j++) {
                for (int i = k + 1; i < n; i++) {
                    X[k][j] -= X[i][j] * L[i][k];
                }
                X[k][j] /= L[k][k];
            }
        }

        return Matrix.SubMatrix(X, n, nx);
    }
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);
    }
1139657
public static void addEquals( DMatrix4 a , DMatrix4 b ) {
        a.a1 += b.a1;
        a.a2 += b.a2;
        a.a3 += b.a3;
        a.a4 += b.a4;
    }
1289063
public void add(double x) {
        Cell[] as; long b, v; int[] hc; Cell a; int n;
        if ((as = cells) != null ||
                !casBase(b = base,
                        Double.doubleToRawLongBits
                                (Double.longBitsToDouble(b) + x))) {
            boolean uncontended = true;
            if ((hc = threadHashCode.get()) == null ||
                    as == null || (n = as.length) < 1 ||
                    (a = as[(n - 1) & hc[0]]) == null ||
                    !(uncontended = a.cas(v = a.value,
                            Double.doubleToRawLongBits
                                    (Double.longBitsToDouble(v) + x))))
                retryUpdate(Double.doubleToRawLongBits(x), hc, uncontended);
        }
    }
698376
public static FloatMatrix solvePositive(FloatMatrix A, FloatMatrix B) {
		A.assertSquare();
		FloatMatrix X = B.dup();
		SimpleBlas.posv('U', A.dup(), X);
		return X;
	}
204025
protected void bicoUpdate(double[] x) {
		assert (!this.bufferPhase && this.numDimensions == x.length);
		// Starts with the global root node as the current root node
		ClusteringTreeNode r = this.root;
		int i = 1;
		while (true) {
			ClusteringTreeNode y = r.nearestChild(x);
			// Checks if the point can not be added to the current level
			if (r.hasNoChildren()
					|| y == null
					|| Metric.distanceSquared(x, y.getCenter()) > calcRSquared(i)) {
				// Creates a new node for the point and adds it to the current
				// root node
				r.addChild(new ClusteringTreeNode(x, new ClusteringFeature(x, calcR(i))));
				this.rootCount++;
				break;
			} else {
				// Checks if the point can be added to the nearest node without
				// exceeding the global threshold
				if (y.getClusteringFeature().calcKMeansCosts(y.getCenter(), x) <= this.T) {
					// Adds the point to the ClusteringFeature
					y.getClusteringFeature().add(1, x, Metric.distanceSquared(x));
					break;
				} else {
					// Navigates one level down in the tree
					r = y;
					i++;
				}
			}
		}
		// Checks if the number of nodes in the tree exceeds the maximum number
		if (this.rootCount > this.maxNumClusterFeatures) {
			rebuild();
		}
	}