Skip to content

Commit

Permalink
Some cleanup.
Browse files Browse the repository at this point in the history
  • Loading branch information
@MathieuDutSik
MathieuDutSik committed Nov 23, 2024
1 parent eae0bd8 commit 758e47c
Show file tree
Hide file tree
Showing 2 changed files with 46 additions and 41 deletions.
24 changes: 14 additions & 10 deletions src_matrix/MAT_Matrix.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,10 @@
#define DEBUG_MAT_MATRIX
#endif

#ifdef SANITY_CHECK
#define SANITY_CHECK_MAT_MATRIX
#endif

template <typename T> struct is_mymatrix<MyMatrix<T>> {
static const bool value = true;
};
Expand Down Expand Up @@ -936,7 +940,7 @@ template <typename T>
MyVector<T> VectorMatrix(MyVector<T> const &eVect, MyMatrix<T> const &eMat) {
int nbCol = eMat.cols();
int nbRow = eMat.rows();
#ifdef DEBUG_MAT_MATRIX
#ifdef SANITY_CHECK_MAT_MATRIX
int n = eVect.size();
if (n != nbRow) {
std::cerr << "n should be equal to nbRow\n";
Expand All @@ -957,7 +961,7 @@ template <typename T>
void AssignMatrixRow(MyMatrix<T> &eMat, int const &iRow,
MyVector<T> const &eVect) {
int nbCol = eMat.cols();
#ifdef DEBUG_MAT_MATRIX
#ifdef SANITY_CHECK_MAT_MATRIX
int n = eVect.size();
if (n != nbCol) {
std::cerr << "We should have eVect.size() = eMat.cols()\n";
Expand Down Expand Up @@ -1073,10 +1077,10 @@ void TMat_Inverse_destroy(MyMatrix<T> &Input, MyMatrix<T> &Output) {
int nbRow = Input.rows();
int nbCol = Input.cols();
T prov1;
#ifdef DEBUG_MAT_MATRIX
#ifdef DEBUG_MAT_MATRIX_DISABLE
std::cerr << "TMat_Inverse_destroy, step 1\n";
#endif
#ifdef DEBUG_MAT_MATRIX
#ifdef SANITY_CHECK_MAT_MATRIX
if (nbRow != nbCol) {
std::cerr << "Error on nbRow, nbCol in TMat_Inverse_destroy";
throw TerminalException{1};
Expand All @@ -1090,12 +1094,12 @@ void TMat_Inverse_destroy(MyMatrix<T> &Input, MyMatrix<T> &Output) {
prov1 = 0;
Output(iRow, iCol) = prov1;
}
#ifdef DEBUG_MAT_MATRIX
#ifdef DEBUG_MAT_MATRIX_DISABLE
std::cerr << "TMat_Inverse_destroy, step 2\n";
#endif
int iColFound;
for (iRow = 0; iRow < nbRow; iRow++) {
#ifdef DEBUG_MAT_MATRIX
#ifdef DEBUG_MAT_MATRIX_DISABLE
std::cerr << "iRow=" << iRow << "\n";
std::cerr << "Input=\n";
WriteMatrix(std::cerr, Input);
Expand All @@ -1110,7 +1114,7 @@ void TMat_Inverse_destroy(MyMatrix<T> &Input, MyMatrix<T> &Output) {
prov1 = 1 / prov1;
}
}
#ifdef DEBUG_MAT_MATRIX
#ifdef SANITY_CHECK_MAT_MATRIX
if (prov1 == 0) {
std::cerr << "Error during the computation of the matrix inverse\n";
throw TerminalException{1};
Expand All @@ -1137,7 +1141,7 @@ void TMat_Inverse_destroy(MyMatrix<T> &Input, MyMatrix<T> &Output) {
std::swap(Input(iRowB, iColFound), Input(iRowB, iRow));
}
}
#ifdef DEBUG_MAT_MATRIX
#ifdef DEBUG_MAT_MATRIX_DISABLE
std::cerr << "TMat_Inverse_destroy, step 3\n";
#endif
}
Expand Down Expand Up @@ -2125,7 +2129,7 @@ template <typename T> struct SolutionMatRepetitive {

public:
SolutionMatRepetitive(MyMatrix<T> const &_TheBasis) : TheBasis(_TheBasis) {
#ifdef DEBUG_MAT_MATRIX
#ifdef SANITY_CHECK_MAT_MATRIX
if (RankMat(TheBasis) != TheBasis.rows()) {
std::cerr << "RankMat(TheBasis)=" << RankMat(TheBasis) << "\n";
std::cerr << " TheBasis.rows()=" << TheBasis.rows() << "\n";
Expand Down Expand Up @@ -2281,7 +2285,7 @@ bool IsSubspaceContained(MyMatrix<T> const &M1, MyMatrix<T> const &M2) {

template <typename T>
bool TestEqualitySpannedSpaces(MyMatrix<T> const &M1, MyMatrix<T> const &M2) {
#ifdef DEBUG_MAT_MATRIX
#ifdef SANITY_CHECK_MAT_MATRIX
if (M1.cols() != M2.cols()) {
std::cerr << "That case is actually not allowed\n";
throw TerminalException{1};
Expand Down
63 changes: 32 additions & 31 deletions src_matrix/MAT_MatrixInt.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,10 @@
#define DEBUG_MATRIX_INT
#endif

#ifdef SANITY_CHECK
#define SANITY_CHECK_MATRIX_INT
#endif

// Now declarations of generic code.
// The code below generally requires the field T to be the ring (or fraction
// ring) of a factorial ring. Operations may work for fields and rings as well.
Expand Down Expand Up @@ -84,7 +88,7 @@ template <typename T> T Int_IndexLattice(MyMatrix<T> const &eMat) {
}
if (IsFirst)
return 0;
#ifdef DEBUG_MATRIX_INT
#ifdef SANITY_CHECK_MATRIX_INT
if (MinPivot == 0) {
std::cerr << "Clear error in the code of IndexLattice\n";
throw TerminalException{1};
Expand Down Expand Up @@ -138,8 +142,8 @@ template <typename T> GCD_dot<T> ComputeGcdDot(MyVector<T> const &x) {
V(v) = V(v) * res.a;
}
}
#ifdef DEBUG_MATRIX_INT
T sum = 0;
#ifdef SANITY_CHECK_MATRIX_INT
T sum(0);
for (int u = 0; u < siz; u++) {
sum += V(u) * x(u);
}
Expand Down Expand Up @@ -230,7 +234,7 @@ ComputePairGcd(T const &m, T const &n) {
Pmat(1, 1) = m / f;
// std::cerr << "Pmat=\n";
// WriteMatrix(std::cerr, Pmat);
#ifdef DEBUG_MATRIX_INT
#ifdef SANITY_CHECK_MATRIX_INT
T diff1 = f - Pmat(0, 0) * m - Pmat(1, 0) * n;
T diff2 = Pmat(0, 1) * m + Pmat(1, 1) * n;
if (diff1 != 0 || diff2 != 0) {
Expand Down Expand Up @@ -261,7 +265,7 @@ ComputePairGcd(T const &m, T const &n) {
Pmat(1, 0) = t;
Pmat(0, 1) = -n / eGCD;
Pmat(1, 1) = m / eGCD;
#ifdef DEBUG_MATRIX_INT
#ifdef SANITY_CHECK_MATRIX_INT
T diff1 = eGCD - Pmat(0, 0) * m - Pmat(1, 0) * n;
T diff2 = Pmat(0, 1) * m + Pmat(1, 1) * n;
if (diff1 != 0 || diff2 != 0) {
Expand Down Expand Up @@ -913,7 +917,7 @@ std::pair<MyMatrix<T>, MyMatrix<T>> SmithNormalForm(MyMatrix<T> const &M) {
MyMatrix<T> H = M;
MyMatrix<T> ROW = IdentityMat<T>(nbRow);
MyMatrix<T> COL = IdentityMat<T>(nbCol);
#ifdef DEBUG_MATRIX_INT
#ifdef SANITY_CHECK_MATRIX_INT
auto check_consistency = [&](std::string const &mesg) -> void {
MyMatrix<T> eProd = ROW * M * COL;
if (eProd != H) {
Expand Down Expand Up @@ -968,7 +972,7 @@ std::pair<MyMatrix<T>, MyMatrix<T>> SmithNormalForm(MyMatrix<T> const &M) {
T TheQ = QuoInt(eVal, ThePivot);
H.row(iRow) -= TheQ * H.row(iRowF);
ROW.row(iRow) -= TheQ * ROW.row(iRowF);
#ifdef DEBUG_MATRIX_INT
#ifdef SANITY_CHECK_MATRIX_INT
mesg = "1 : Error_at iRowF=" + std::to_string(iRowF) +
" iColF=" + std::to_string(iColF) +
" iRpw=" + std::to_string(iRow) +
Expand All @@ -987,7 +991,7 @@ std::pair<MyMatrix<T>, MyMatrix<T>> SmithNormalForm(MyMatrix<T> const &M) {
T TheQ = QuoInt(eVal, ThePivot);
H.col(iCol) -= TheQ * H.col(iColF);
COL.col(iCol) -= TheQ * COL.col(iColF);
#ifdef DEBUG_MATRIX_INT
#ifdef SANITY_CHECK_MATRIX_INT
mesg = "2 : Error_at iRowF=" + std::to_string(iRowF) +
" iColF=" + std::to_string(iColF) +
" iCol=" + std::to_string(iCol) +
Expand All @@ -1004,7 +1008,7 @@ std::pair<MyMatrix<T>, MyMatrix<T>> SmithNormalForm(MyMatrix<T> const &M) {
MyMatrix<T> Trans = TranspositionMatrix<T>(nbRow, posDone, iRowF);
ROW = Trans * ROW;
H = Trans * H;
#ifdef DEBUG_MATRIX_INT
#ifdef SANITY_CHECK_MATRIX_INT
mesg = "3 : Error_at iRowF=" + std::to_string(iRowF) +
" posDone=" + std::to_string(posDone);
check_consistency(mesg);
Expand All @@ -1014,7 +1018,7 @@ std::pair<MyMatrix<T>, MyMatrix<T>> SmithNormalForm(MyMatrix<T> const &M) {
MyMatrix<T> Trans = TranspositionMatrix<T>(nbCol, posDone, iColF);
COL = COL * Trans;
H = H * Trans;
#ifdef DEBUG_MATRIX_INT
#ifdef SANITY_CHECK_MATRIX_INT
mesg = "4 : Error_at iColF=" + std::to_string(iColF) +
" posDone=" + std::to_string(posDone);
check_consistency(mesg);
Expand Down Expand Up @@ -1136,7 +1140,7 @@ template <typename T> MyMatrix<T> NullspaceIntTrMat(MyMatrix<T> const &eMat) {
nbFound++;
}
}
#ifdef DEBUG_MATRIX_INT
#ifdef SANITY_CHECK_MATRIX_INT
if (nbFound == 0) {
std::cerr << "The column is zero. No work possible\n";
throw TerminalException{1};
Expand Down Expand Up @@ -1255,7 +1259,7 @@ template <typename T> MyMatrix<T> NullspaceIntTrMat(MyMatrix<T> const &eMat) {
idx++;
}
}
#ifdef DEBUG_MATRIX_INT
#ifdef SANITY_CHECK_MATRIX_INT
size_t nbRow = eMat.rows();
for (size_t iVect = 0; iVect < dimSpace; iVect++)
for (size_t iRow = 0; iRow < nbRow; iRow++) {
Expand Down Expand Up @@ -1319,14 +1323,14 @@ template <typename T> MyMatrix<T> ComplementToBasis(MyVector<T> const &TheV) {
}
}
}
#ifdef DEBUG_MATRIX_INT
#ifdef SANITY_CHECK_MATRIX_INT
if (idxSelect == -1) {
std::cerr << "Inconsistency in computation of value\n";
throw TerminalException{1};
}
#endif
if (nbDiffZero == 1) {
#ifdef DEBUG_MATRIX_INT
#ifdef SANITY_CHECK_MATRIX_INT
if (AbsVal != 1) {
std::cerr << "Wrong value for AbsVal\n";
throw TerminalException{1};
Expand Down Expand Up @@ -1370,7 +1374,7 @@ template <typename T> MyMatrix<T> ComplementToBasis(MyVector<T> const &TheV) {
T eVal = TheVcopy(i2) + eCoeff * TheVcopy(i1);
TheVcopy(i2) = eVal;
}
#ifdef DEBUG_MATRIX_INT
#ifdef SANITY_CHECK_MATRIX_INT
if (!TestEquality(TheVcopy, TheV)) {
std::cerr << "TheVcopy =";
WriteVectorNoDim(std::cerr, TheVcopy);
Expand Down Expand Up @@ -1597,9 +1601,10 @@ std::optional<MyVector<T>> SolutionIntMat(MyMatrix<T> const &TheMat,
}
}
}
if (nbDiff == 1 || nbDiff == 0)
if (nbDiff == 1 || nbDiff == 0) {
break;
#ifdef DEBUG_MATRIX_INT
}
#ifdef SANITY_CHECK_MATRIX_INT
if (MinValue == 0) {
std::cerr << "MinValue should not be zero\n";
throw TerminalException{1};
Expand All @@ -1617,7 +1622,7 @@ std::optional<MyVector<T>> SolutionIntMat(MyMatrix<T> const &TheMat,
}
}
if (nbDiff == 1) {
#ifdef DEBUG_MATRIX_INT
#ifdef SANITY_CHECK_MATRIX_INT
if (iVectFound == -1) {
std::cerr << "Clear error in the program\n";
throw TerminalException{1};
Expand Down Expand Up @@ -1657,7 +1662,7 @@ CanSolIntMat<T> ComputeCanonicalFormFastReduction(MyMatrix<T> const &TheMat) {
int nbDiff;
int nbRow = TheMat.rows();
int nbCol = TheMat.cols();
#ifdef DEBUG_MATRIX_INT
#ifdef SANITY_CHECK_MATRIX_INT
if (nbRow == 0) {
std::cerr << "Need to write the code here\n";
throw TerminalException{1};
Expand Down Expand Up @@ -1693,7 +1698,7 @@ CanSolIntMat<T> ComputeCanonicalFormFastReduction(MyMatrix<T> const &TheMat) {
}
if (nbDiff == 1 || nbDiff == 0)
break;
#ifdef DEBUG_MATRIX_INT
#ifdef SANITY_CHECK_MATRIX_INT
if (MinValue == 0) {
std::cerr << "MinValue should not be zero\n";
throw TerminalException{1};
Expand All @@ -1713,7 +1718,7 @@ CanSolIntMat<T> ComputeCanonicalFormFastReduction(MyMatrix<T> const &TheMat) {
int eVal;
if (nbDiff == 1) {
eVal = iVectFound;
#ifdef DEBUG_MATRIX_INT
#ifdef SANITY_CHECK_MATRIX_INT
if (iVectFound == -1) {
std::cerr << "Clear error in the program\n";
throw TerminalException{1};
Expand Down Expand Up @@ -1932,8 +1937,7 @@ AffineBasisResult Kernel_ComputeAffineBasis(MyMatrix<T> const &EXT) {
if (eCol == miss_val && EXTwork(iVect, iCol) != 0 &&
ColumnStatus[iCol] == 1)
eCol = iCol;
#ifdef DEBUG_MATRIX_INT
std::cerr << "eCol=" << eCol << "\n";
#ifdef SANITY_CHECK_MATRIX_INT
if (eCol == miss_val) {
std::cerr << "This should not be selected\n";
std::cerr << "nbCol=" << nbCol << "\n";
Expand Down Expand Up @@ -2185,7 +2189,7 @@ template <typename T> MyMatrix<T> GetZbasis(MyMatrix<T> const &ListElement) {
// std::cerr << "After TheRedMat construction\n";
MyMatrix<T> NSP = NullspaceIntMat(TheRedMat);
// std::cerr << "We have NSP\n";
#ifdef DEBUG_MATRIX_INT
#ifdef SANITY_CHECK_MATRIX_INT
if (NSP.rows() != 1) {
std::cerr << "|NSP|=" << NSP.rows() << " when it should be 1\n";
std::cerr << "TheRedMat:\n";
Expand Down Expand Up @@ -2286,15 +2290,12 @@ template <typename T> MyMatrix<T> GetZbasis(MyMatrix<T> const &ListElement) {
};
fComputeSpeed();
int nbElt = ListElement.rows();
// std::cerr << "nbElt=" << nbElt << "\n";
for (int iElt = 0; iElt < nbElt; iElt++) {
// std::cerr << "iElt=" << iElt << "\n";
MyVector<T> eElt = GetMatrixRow(ListElement, iElt);
fInsert(eElt);
// std::cerr << "After fInsert\n";
}

#ifdef DEBUG_MATRIX_INT
#ifdef SANITY_CHECK_MATRIX_INT
int DimSpace = TheBasis.rows();
for (int iBas = 0; iBas < DimSpace; iBas++) {
MyVector<T> eLine = GetMatrixRow(TheBasis, iBas);
Expand Down Expand Up @@ -2370,14 +2371,14 @@ MyMatrix<T> IntersectionLattice_VectorSpace(MyMatrix<T> const &Latt,
V[i] = i + n_spa;
MyMatrix<T> Latt3 = SelectColumn(Latt2, V);
MyMatrix<T> NSP = NullspaceIntMat(Latt3);
#ifdef DEBUG_MATRIX_INT
#ifdef SANITY_CHECK_MATRIX_INT
if (!IsIntegralMatrix(NSP)) {
std::cerr << "NSP should be integral\n";
throw TerminalException{1};
}
#endif
MyMatrix<T> IntBasis = NSP * Latt;
#ifdef DEBUG_MATRIX_INT
#ifdef SANITY_CHECK_MATRIX_INT
for (int i_s = 0; i_s < n_spa; i_s++) {
MyVector<T> v = GetMatrixRow(IntBasis, i_s);
std::optional<MyVector<T>> opt1 = SolutionIntMat(Latt, v);
Expand Down Expand Up @@ -2535,7 +2536,7 @@ T GetDenominatorQuotientSolution(MyMatrix<T> const& Amat) {
template<typename T>
MyMatrix<T> IntegralSpaceSaturation(MyMatrix<T> const& TheSpace) {
int nbRow = TheSpace.rows();
#ifdef DEBUG_MATRIX_INT
#ifdef SANITY_CHECK_MATRIX_INT
if (nbRow != RankMat(TheSpace)) {
std::cerr << "TheSpace should be full dimensional\n";
throw TerminalException{1};
Expand Down

0 comments on commit 758e47c

Please sign in to comment.