Solve A*X = b

@param b A column vector with as many rows as A.
@return X so that L*L^T*X = b
@exception IllegalArgumentException Matrix row dimensions must agree.
@exception RuntimeException Matrix is not symmetric positive definite.
149
public double[] solve(double[] b) {
    if(b.length != L.length) {
      throw new IllegalArgumentException(ERR_MATRIX_DIMENSIONS);
    }
    if(!isspd) {
      throw new ArithmeticException(ERR_MATRIX_NOT_SPD);
    }
    // Work on a copy!
    return solveLtransposed(solveLInplace(copy(b)));
  }
1269769
public static DoubleMatrix solveSymmetric(DoubleMatrix A, DoubleMatrix B) {
		A.assertSquare();
		DoubleMatrix X = B.dup();
		int[] ipiv = new int[B.rows];
		SimpleBlas.sysv('U', A.dup(), ipiv, X);
		return X;
	}
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));
            }
        }
    }
1131033
public function solve(Matrix $B)
    {
        if ($B->getRows() != $this->m) {
            throw MatrixException::notSquareMatrix();
        }

        if (!$this->isNonsingular()) {
            throw MatrixException::singularMatrix();
        }

        // Copy right hand side with pivoting
        $nx = $B->getColumns();
        $X  = $this->getSubMatrix($B->toArray(), $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[$i][$j] -= $X[$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[$k][$j] /= $this->LU[$k][$k];
            }
            for ($i = 0; $i < $k; ++$i) {
                for ($j = 0; $j < $nx; ++$j) {
                    $X[$i][$j] -= $X[$k][$j] * $this->LU[$i][$k];
                }
            }
        }
        return $X;
    }
213036
public void solve(double[] b, double[] x) {
        if (!full) {
            throw new IllegalStateException("This is not a FULL singular value decomposition.");
        }

        if (b.length != m || x.length != n) {
            throw new IllegalArgumentException("Dimensions do not agree.");
        }

        double[] tmp = new double[n];
        for (int j = 0; j < n; j++) {
            double r = 0.0;
            if (s[j] > tol) {
                for (int i = 0; i < m; i++) {
                    r += U.get(i, j) * b[i];
                }
                r /= s[j];
            }
            tmp[j] = r;
        }

        for (int j = 0; j < n; j++) {
            double r = 0.0;
            for (int jj = 0; jj < n; jj++) {
                r += V.get(j, jj) * tmp[jj];
            }
            x[j] = r;
        }
    }
1248419
public Matrix plus(Matrix B) {
        Matrix A = this;
        if (B.M != A.M || B.N != A.N) {
            throw new RuntimeException("Illegal matrix dimensions.");
        }
        Matrix C = new Matrix(M, N);
        for (int i = 0; i < M; i++) {
            for (int j = 0; j < N; j++) {
                C.data[i][j] = A.data[i][j] + B.data[i][j];
            }
        }
        return C;
    }
225751
public static final DoubleMatrix2D add(DoubleMatrix2D A, DoubleMatrix2D B){
		if(A.rows()!=B.rows() || A.columns()!=B.columns()){
			throw new IllegalArgumentException("wrong matrices dimensions");
		}
		DoubleMatrix2D ret = DoubleFactory2D.dense.make(A.rows(), A.columns());
		for(int i=0; i<ret.rows(); i++){
			for(int j=0; j<ret.columns(); j++){
				ret.setQuick(i, j, A.getQuick(i, j) + B.getQuick(i, j));
			}
		}

		return ret;
	}
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;
  }
1077707
public Matrix multiplyTranspose(final Matrix B)
    {
        Matrix C = new DenseMatrix(this.rows(), B.rows());
        multiplyTranspose(B, C);
        return C;
    }
563073
public void solve(DenseMatrix B) {
        if (B.nrows() != qr.nrows()) {
            throw new IllegalArgumentException(String.format("Row dimensions do not agree: A is %d x %d, but B is %d x %d", qr.nrows(), qr.nrows(), B.nrows(), B.ncols()));
        }

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

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

        // Compute Y = Q' * B
        for (int k = 0; k < n; k++) {
            for (int j = 0; j < nrhs; j++) {
                double s = 0.0;
                for (int i = k; i < m; i++) {
                    s += qr.get(i, k) * B.get(i, j);
                }
                s = -s / qr.get(k, k);
                for (int i = k; i < m; i++) {
                    B.add(i, j, s * qr.get(i, k));
                }
            }
        }

        // Solve R*X = Y;
        for (int k = n - 1; k >= 0; k--) {
            for (int j = 0; j < nrhs; j++) {
                B.set(k, j, B.get(k, j) / tau[k]);
            }
            
            for (int i = 0; i < k; i++) {
                for (int j = 0; j < nrhs; j++) {
                    B.sub(i, j, B.get(k, j) * qr.get(i, k));
                }
            }
        }
    }
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];
            }
        }
    }
614787
def solve(self, b):
        
        #b = self.permute_vec(b)
        y = []
        for bri in range(self.N):  # block row index
            for li in range(self.n):  # local row index
                s = 0.0
                for lci in range(li): # local column index
                    s += self.lu[bri][li, lci]*y[bri*self.n+lci]
                for di in range(1, self.ndiag+1):
                    if bri >= di:
                        # di:th sub diagonal (counted as distance from main diag)
                        ci = self.colbyrow[bri][li]
                        s += self.sub[di-1][(bri-di)*self.n+ci]*y[
                            (bri-di)*self.n + ci]
                y.append(b[bri*self.n+self.rowbycol[bri][li]]-s) # Doolittle: L[i, i] == 1
        x = [0]*len(y)
        for bri in range(self.N-1, -1, -1):
            for li in range(self.n - 1, -1, -1):
                s = 0.0
                for ci in range(li+1, self.n):
                    s += self.lu[bri][li, ci]*x[bri*self.n + ci]
                for di in range(1, self.ndiag+1):
                    if bri < self.N-di:
                        ci = self.colbyrow[bri][li]
                        s += self.sup[di-1][bri*self.n+ci]*x[(bri+di)*self.n + ci]
                x[bri*self.n+li] = (y[bri*self.n + li] - s)/self.lu[bri][li, li]
        return x
580596
public boolean eq(final Matrix B)
    {
        final Matrix A = this;

        if ((B.m_rows != A.m_rows) || (B.m_columns != A.m_columns))
        {
            return false;
        }
        for (int i = 0; i < m_rows; i++)
        {
            for (int j = 0; j < m_columns; j++)
            {
                if (A.m_data[i][j] != B.m_data[i][j])
                {
                    return false;
                }
            }
        }
        return true;
    }
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);
        }
    }
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]);
            }
        }
    }
871488
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());
		}

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

		// Solve L.y = b
		final double[] y = new double[dim];
		// L is upper left diagonal
		for (int i = 0; i < diagonalLength; i++) {
			double LII = LData[0][i];
			y[i] = b.getQuick(i) / LII;
		}
		for (int i = diagonalLength; i < dim; i++) {
			double[] LI = LData[i - diagonalLength + 1];
			double LII = LI[i];
			double sum = 0;
			for (int j = 0; j < i; j++) {
				sum += LI[j] * y[j];
			}
			y[i] = (b.getQuick(i) - sum) / LII;
		}

		// logger.debug("b: " + ArrayUtils.toString(b));
		// logger.debug("L.y: " + ArrayUtils.toString(getL().operate(y)));

		// Solve L[T].x = y
		final DoubleMatrix1D x = F1.make(dim);
		// for (int i = dim-1; i > -1; i--) {
		// double sum = 0;
		// int ll = Math.max(i, diagonalLength-1);
		// for (int j = dim-1; j > ll; j--) {
		// sum += LData[j][i] * x.getQuick(j);
		//
		// }
		// x.setQuick(i, (y[i] - sum)/ LData[i][i]);
		// }
		for (int i = dim - 1; i > diagonalLength - 1; i--) {
			double sum = 0;
			for (int j = dim - 1; j > i; j--) {
				sum += LData[j - diagonalLength + 1][i] * x.getQuick(j);

			}
			x.setQuick(i, (y[i] - sum) / LData[i - diagonalLength + 1][i]);
		}
		for (int i = diagonalLength - 1; i > -1; i--) {
			double sum = 0;
			for (int j = dim - 1; j > diagonalLength - 1; j--) {
				sum += LData[j - diagonalLength + 1][i] * x.getQuick(j);

			}
			x.setQuick(i, (y[i] - sum) / LData[0][i]);
		}

		if (this.rescaler != null) {
			// return ALG.mult(this.U, x);
			return ColtUtils.diagonalMatrixMult(this.U, x);
		} else {
			return x;
		}
	}
500024
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());
		}

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

		// Solve L.y = b
		final double[] y = new double[dim];
		for (int i = 0; i < dim; i++) {
			double[] LI = LData[i];
			double sum = 0;
			for (int j = 0; j < i; j++) {
				sum += LI[j] * y[j];
			}
			y[i] = (b.getQuick(i) - sum) / LI[i];
		}

		// logger.debug("b: " + ArrayUtils.toString(b));
		// logger.debug("L.y: " + ArrayUtils.toString(getL().operate(y)));

		// Solve D.z = y
		final double[] z = new double[dim];
		for (int i = 0; i < dim; i++) {
			z[i] = y[i] / DData[i];
		}

		// Solve L[T].x = z
		final DoubleMatrix1D x = F1.make(dim);
		for (int i = dim - 1; i > -1; i--) {
			double sum = 0;
			for (int j = dim - 1; j > i; j--) {
				sum += LData[j][i] * x.getQuick(j);
			}
			x.setQuick(i, (z[i] - sum) / LData[i][i]);
		}

		if (this.rescaler != null) {
			// return ALG.mult(this.U, x);
			return ColtUtils.diagonalMatrixMult(this.U, x);
		} else {
			return x;
		}
	}