Skip to content

Commit

Permalink
Address some bugs, more remain.
Browse files Browse the repository at this point in the history
  • Loading branch information
@MathieuDutSik
MathieuDutSik committed Sep 8, 2023
1 parent c2c65cc commit 353e25b
Show file tree
Hide file tree
Showing 2 changed files with 43 additions and 3 deletions.
13 changes: 12 additions & 1 deletion src_matrix/MAT_Matrix_SubsetSolver.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,7 @@ struct SubsetRankOneSolver_Acceleration {
//
// Faster modular version of EXT_red
//
std::cerr << "SubsetRankOneSolver_Acceleration, constructor, step 1\n";
max_bits = 0;
EXT_fast = MyMatrix<Tfast>(nbRow, nbCol);
EXT_lift = MyMatrix<Tlift>(nbRow, nbCol);
Expand All @@ -69,13 +70,16 @@ struct SubsetRankOneSolver_Acceleration {
}
try_int = (max_bits <= 30);
max_bits += get_bit(static_cast<int64_t>(nbCol));
std::cerr << "SubsetRankOneSolver_Acceleration, constructor, step 2\n";
}

MyVector<Tint> GetKernelVector(Face const& sInc) {
std::cerr << "SubsetRankOneSolver_Acceleration, GetKernelVector, step 1\n";
size_t nb = sInc.count();
MyVector<Tint> Vkernel(nbCol);
bool failed_int = false;
if (try_int) {
std::cerr << "SubsetRankOneSolver_Acceleration, GetKernelVector, step 2\n";
boost::dynamic_bitset<>::size_type jRow = sInc.find_first();
auto f = [&](MyMatrix<Tfast> &M, size_t eRank,
[[maybe_unused]] size_t iRow) -> void {
Expand All @@ -84,6 +88,7 @@ struct SubsetRankOneSolver_Acceleration {
};
MyVector<Tfast> Vzero_Tfast =
NullspaceTrMatTargetOne_Kernel<Tfast, decltype(f)>(nb, nbCol, f);
std::cerr << "SubsetRankOneSolver_Acceleration, GetKernelVector, step 3\n";
// check result at full precision in case of overflows
bool allzero = true;
for (int iCol = 0; iCol < nbCol; iCol++) {
Expand All @@ -92,23 +97,28 @@ struct SubsetRankOneSolver_Acceleration {
break;
}
}
std::cerr << "SubsetRankOneSolver_Acceleration, GetKernelVector, step 4\n";
if (allzero) {
failed_int = true;
} else {
std::cerr << "SubsetRankOneSolver_Acceleration, GetKernelVector, step 5\n";
MyVector<Tlift> VZ_lift(nbCol);
// reconstruct the vector
size_t max_bits_NSP = 0;
std::cerr << "nbCol=" << nbCol << " |Vzero_Tfast|=" << Vzero_Tfast.size() << " |lifts|=" << lifts.size() << "\n";
lifts[0] = Vzero_Tfast(0, 0).rational_lift();
Tlift lcm = lifts[0].second;
for (int iCol = 1; iCol < nbCol; iCol++) {
lifts[iCol] = Vzero_Tfast(0, iCol).rational_lift();
lifts[iCol] = Vzero_Tfast(iCol).rational_lift();
lcm = LCMpair(lcm, lifts[iCol].second);
}
std::cerr << "SubsetRankOneSolver_Acceleration, GetKernelVector, step 6\n";
for (int iCol = 0; iCol < nbCol; iCol++) {
VZ_lift(iCol) = lifts[iCol].first * (lcm / lifts[iCol].second);
Vkernel(iCol) = UniversalScalarConversion<Tint,Tlift>(VZ_lift(iCol));
max_bits_NSP = std::max(max_bits_NSP, get_bit(VZ_lift(iCol)));
}
std::cerr << "SubsetRankOneSolver_Acceleration, GetKernelVector, step 7\n";
// check if elements are small enough to do computation in
if (max_bits + max_bits_NSP <= 60) {
// check if part of kernel
Expand All @@ -124,6 +134,7 @@ struct SubsetRankOneSolver_Acceleration {
}
jRow = sInc.find_next(jRow);
}
std::cerr << "SubsetRankOneSolver_Acceleration, GetKernelVector, step 8\n";
} else {
failed_int = true;
}
Expand Down
33 changes: 31 additions & 2 deletions src_matrix/MatrixSubsetSolver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@ template<typename T>
inline typename std::enable_if<!has_reduction_subset_solver<T>::value, void>::type process(std::string const& FileI, std::string const& FileO) {
MyMatrix<T> TheMat = ReadMatrixFile<T>(FileI);
int n_row = TheMat.rows();
int n_col = TheMat.cols();
// This is a test case
SubsetRankOneSolver_Field<T> solver(TheMat);
for (int i_row=0; i_row<n_row; i_row++) {
Expand All @@ -24,7 +25,20 @@ inline typename std::enable_if<!has_reduction_subset_solver<T>::value, void>::ty
}
}
MyVector<T> V = solver.GetKernelVector(f);
std::cerr << "i_row=" << i_row << " V=" << StringVector(V) << "\n";
for (int j_row=0; j_row<n_row; j_row++) {
if (j_row != i_row) {
T scal(0);
for (int i_col=0; i_col<n_col; i_col++) {
scal += V(i_col) * TheMat(j_row,i_col);
}
if (scal != 0) {
std::cerr << "j_row=" << j_row << " V=" << StringVector(V) << "\n";
std::cerr << "Wrong scalar product\n";
throw TerminalException{1};
}
}
}
std::cerr << "SCH A : i_row=" << i_row << " V=" << StringVector(V) << "\n";
}
}

Expand All @@ -33,6 +47,7 @@ inline typename std::enable_if<has_reduction_subset_solver<T>::value, void>::typ
using Tint = typename SubsetRankOneSolver_Acceleration<T>::Tint;
MyMatrix<Tint> TheMat = ReadMatrixFile<Tint>(FileI);
int n_row = TheMat.rows();
int n_col = TheMat.cols();
// This is a test case
SubsetRankOneSolver_Acceleration<T> solver(TheMat);
for (int i_row=0; i_row<n_row; i_row++) {
Expand All @@ -43,7 +58,20 @@ inline typename std::enable_if<has_reduction_subset_solver<T>::value, void>::typ
}
}
MyVector<Tint> V = solver.GetKernelVector(f);
std::cerr << "i_row=" << i_row << " V=" << StringVector(V) << "\n";
for (int j_row=0; j_row<n_row; j_row++) {
if (j_row != i_row) {
Tint scal(0);
for (int i_col=0; i_col<n_col; i_col++) {
scal += V(i_col) * TheMat(j_row,i_col);
}
if (scal != 0) {
std::cerr << "j_row=" << j_row << " V=" << StringVector(V) << "\n";
std::cerr << "Wrong scalar product\n";
throw TerminalException{1};
}
}
}
std::cerr << "SCH B : i_row=" << i_row << " V=" << StringVector(V) << "\n";
}
}

Expand Down Expand Up @@ -73,6 +101,7 @@ int main(int argc, char *argv[]) {
return process<T>(FileI, FileO);
}
std::cerr << "Failed to find a matching type\n";
std::cerr << "Allowed types are mpq_class, cpp_rational, safe_rational\n";
throw TerminalException{1};
};
f();
Expand Down

0 comments on commit 353e25b

Please sign in to comment.