From 758e47c102500ec9eef32d9e1e76a8ba57a77231 Mon Sep 17 00:00:00 2001 From: Mathieu Dutour Sikiric Date: Sat, 23 Nov 2024 20:00:01 +0100 Subject: [PATCH] Some cleanup. --- src_matrix/MAT_Matrix.h | 24 +++++++++------ src_matrix/MAT_MatrixInt.h | 63 +++++++++++++++++++------------------- 2 files changed, 46 insertions(+), 41 deletions(-) diff --git a/src_matrix/MAT_Matrix.h b/src_matrix/MAT_Matrix.h index 25e37ca..50d5421 100644 --- a/src_matrix/MAT_Matrix.h +++ b/src_matrix/MAT_Matrix.h @@ -28,6 +28,10 @@ #define DEBUG_MAT_MATRIX #endif +#ifdef SANITY_CHECK +#define SANITY_CHECK_MAT_MATRIX +#endif + template struct is_mymatrix> { static const bool value = true; }; @@ -936,7 +940,7 @@ template MyVector VectorMatrix(MyVector const &eVect, MyMatrix 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"; @@ -957,7 +961,7 @@ template void AssignMatrixRow(MyMatrix &eMat, int const &iRow, MyVector 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"; @@ -1073,10 +1077,10 @@ void TMat_Inverse_destroy(MyMatrix &Input, MyMatrix &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}; @@ -1090,12 +1094,12 @@ void TMat_Inverse_destroy(MyMatrix &Input, MyMatrix &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); @@ -1110,7 +1114,7 @@ void TMat_Inverse_destroy(MyMatrix &Input, MyMatrix &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}; @@ -1137,7 +1141,7 @@ void TMat_Inverse_destroy(MyMatrix &Input, MyMatrix &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 } @@ -2125,7 +2129,7 @@ template struct SolutionMatRepetitive { public: SolutionMatRepetitive(MyMatrix 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"; @@ -2281,7 +2285,7 @@ bool IsSubspaceContained(MyMatrix const &M1, MyMatrix const &M2) { template bool TestEqualitySpannedSpaces(MyMatrix const &M1, MyMatrix 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}; diff --git a/src_matrix/MAT_MatrixInt.h b/src_matrix/MAT_MatrixInt.h index 5ea2b5b..ab115e4 100644 --- a/src_matrix/MAT_MatrixInt.h +++ b/src_matrix/MAT_MatrixInt.h @@ -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. @@ -84,7 +88,7 @@ template T Int_IndexLattice(MyMatrix 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}; @@ -138,8 +142,8 @@ template GCD_dot ComputeGcdDot(MyVector 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); } @@ -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) { @@ -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) { @@ -913,7 +917,7 @@ std::pair, MyMatrix> SmithNormalForm(MyMatrix const &M) { MyMatrix H = M; MyMatrix ROW = IdentityMat(nbRow); MyMatrix COL = IdentityMat(nbCol); -#ifdef DEBUG_MATRIX_INT +#ifdef SANITY_CHECK_MATRIX_INT auto check_consistency = [&](std::string const &mesg) -> void { MyMatrix eProd = ROW * M * COL; if (eProd != H) { @@ -968,7 +972,7 @@ std::pair, MyMatrix> SmithNormalForm(MyMatrix 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) + @@ -987,7 +991,7 @@ std::pair, MyMatrix> SmithNormalForm(MyMatrix 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) + @@ -1004,7 +1008,7 @@ std::pair, MyMatrix> SmithNormalForm(MyMatrix const &M) { MyMatrix Trans = TranspositionMatrix(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); @@ -1014,7 +1018,7 @@ std::pair, MyMatrix> SmithNormalForm(MyMatrix const &M) { MyMatrix Trans = TranspositionMatrix(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); @@ -1136,7 +1140,7 @@ template MyMatrix NullspaceIntTrMat(MyMatrix 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}; @@ -1255,7 +1259,7 @@ template MyMatrix NullspaceIntTrMat(MyMatrix 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++) { @@ -1319,14 +1323,14 @@ template MyMatrix ComplementToBasis(MyVector 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}; @@ -1370,7 +1374,7 @@ template MyMatrix ComplementToBasis(MyVector 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); @@ -1597,9 +1601,10 @@ std::optional> SolutionIntMat(MyMatrix 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}; @@ -1617,7 +1622,7 @@ std::optional> SolutionIntMat(MyMatrix 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}; @@ -1657,7 +1662,7 @@ CanSolIntMat ComputeCanonicalFormFastReduction(MyMatrix 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}; @@ -1693,7 +1698,7 @@ CanSolIntMat ComputeCanonicalFormFastReduction(MyMatrix 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}; @@ -1713,7 +1718,7 @@ CanSolIntMat ComputeCanonicalFormFastReduction(MyMatrix 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}; @@ -1932,8 +1937,7 @@ AffineBasisResult Kernel_ComputeAffineBasis(MyMatrix 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"; @@ -2185,7 +2189,7 @@ template MyMatrix GetZbasis(MyMatrix const &ListElement) { // std::cerr << "After TheRedMat construction\n"; MyMatrix 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"; @@ -2286,15 +2290,12 @@ template MyMatrix GetZbasis(MyMatrix 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 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 eLine = GetMatrixRow(TheBasis, iBas); @@ -2370,14 +2371,14 @@ MyMatrix IntersectionLattice_VectorSpace(MyMatrix const &Latt, V[i] = i + n_spa; MyMatrix Latt3 = SelectColumn(Latt2, V); MyMatrix 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 IntBasis = NSP * Latt; -#ifdef DEBUG_MATRIX_INT +#ifdef SANITY_CHECK_MATRIX_INT for (int i_s = 0; i_s < n_spa; i_s++) { MyVector v = GetMatrixRow(IntBasis, i_s); std::optional> opt1 = SolutionIntMat(Latt, v); @@ -2535,7 +2536,7 @@ T GetDenominatorQuotientSolution(MyMatrix const& Amat) { template MyMatrix IntegralSpaceSaturation(MyMatrix 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};