Function getRhoJ

Applications/HHL_Algorithm/BasicFunctionByOrigin.cpp:163–201

}

int getRhoJ(VectorXd& r, double dd, vector<int>& No0Ind, vector<int>& SpIni, MatrixXd& Ae) {

	int size = r.size();
	int rows = Ae.rows();

	VectorXd r_all(rows);
	vector<int> Jset;
	//vector<int> JsetNew;

	for (int i = 0; i < size; i++) {
		if (r(i) != 0) {
			Jset.push_back(No0Ind[i]);
		}
	}

	for (auto it = SpIni.begin(); it != SpIni.end(); ++it) {
		Jset.erase(std::remove(Jset.begin(), Jset.end(), *it), Jset.end());
	}

	MatrixXd sA(rows, No0Ind.size());
	for (int i = 0; i < No0Ind.size(); i++) {
		sA.col(i) = Ae.col(No0Ind[i]);
	}
	//vector<double> rho_j;

	double min_rho = 10;
	int J_add=0;
	for (int j = 0; j < Jset.size(); j++) {
		
		auto up = pow(r.dot(sA.row(Jset[j])), 2);
		auto down = pow(sA.row(Jset[j]).norm(), 2);
		auto rho = dd * dd - up / down;
		if (rho < min_rho) {
			min_rho = rho;
			J_add = j;
		}
	}
	return Jset[J_add];
}

int getRhoJ(VectorXd& r, double dd, vector<int>& No0Ind, vector<int>& SpIni, SparseQMatrix<double>& Ae) {