From e380b4c07baf4846f9abc5f50c896c5540f42069 Mon Sep 17 00:00:00 2001
From: labourasse <labourassee@gmail.com>
Date: Mon, 18 Nov 2024 18:14:06 +0100
Subject: [PATCH] remove couts
---
src/algebra/EigenvalueSolver.cpp | 96 +++++++++++++++-----------------
1 file changed, 46 insertions(+), 50 deletions(-)
diff --git a/src/algebra/EigenvalueSolver.cpp b/src/algebra/EigenvalueSolver.cpp
index efaaa19a2..28ec7131d 100644
--- a/src/algebra/EigenvalueSolver.cpp
+++ b/src/algebra/EigenvalueSolver.cpp
@@ -216,10 +216,10 @@ EigenvalueSolver::rowReduce(const TinyMatrix<3, 3, T>& matrix, const double epsi
if (pivotRow == 3) {
continue;
}
- std::cout << "before: "
- << "i = " << i << " pivotRow = " << pivotRow << " " << redmat << "\n";
+ // std::cout << "before: "
+ // << "i = " << i << " pivotRow = " << pivotRow << " " << redmat << "\n";
redmat = swap(redmat, i, pivotRow);
- std::cout << "after: " << redmat << "\n";
+ // std::cout << "after: " << redmat << "\n";
double pivotValue = redmat(i, i);
for (size_t j = i; j < 3; ++j) {
@@ -244,8 +244,8 @@ EigenvalueSolver::findEigenvector(const TinyMatrix<3, 3>& A,
{
TinyMatrix<3, 3> eigenvectors = zero;
if (space_size == 3) {
- std::cout << "space_size = 3 "
- << "\n";
+ // std::cout << "space_size = 3 "
+ // << "\n";
return eigenvectors0;
}
TinyMatrix<3, 3> B = A;
@@ -254,19 +254,19 @@ EigenvalueSolver::findEigenvector(const TinyMatrix<3, 3>& A,
B(i, i) -= eigenvalue;
}
- std::cout << " Avant " << B << "\n"
- << "\n";
+ // std::cout << " Avant " << B << "\n"
+ // << "\n";
B = rowReduce(B, epsilon);
- std::cout << " Après " << B << "\n"
- << "\n";
+ // std::cout << " Après " << B << "\n"
+ // << "\n";
for (size_t i = 0; i < 3; ++i) {
bool isFreeVariableRow = true;
for (size_t j = 0; j < 3; ++j) {
if (std::fabs(B(i, j)) > epsilon) {
isFreeVariableRow = false;
- std::cout << " B(" << i << "," << j << ") = " << B(i, j) << "\n";
+ // std::cout << " B(" << i << "," << j << ") = " << B(i, j) << "\n";
break;
}
}
@@ -274,9 +274,9 @@ EigenvalueSolver::findEigenvector(const TinyMatrix<3, 3>& A,
if (isFreeVariableRow) {
TinyVector<3> eigenvector = zero;
eigenvector[i] = 1.0;
- std::cout << i << " is free "
- << "\n";
- std::cout << std::flush;
+ // std::cout << i << " is free "
+ // << "\n";
+ // std::cout << std::flush;
for (int k = i - 1; k >= 0; --k) {
double sum = 0.0;
for (size_t j = k + 1; j < 3; ++j) {
@@ -292,14 +292,14 @@ EigenvalueSolver::findEigenvector(const TinyMatrix<3, 3>& A,
}
}
}
- std::cout << " Before orthorgonalization " << B << "\n";
- std::cout << " Before orthorgonalization " << eigenvectors << "\n";
+ // std::cout << " Before orthorgonalization " << B << "\n";
+ // std::cout << " Before orthorgonalization " << eigenvectors << "\n";
if (space_size == 2) {
TinyVector<3> first = zero;
TinyVector<3> second = zero;
// TinyVector<3> new_second =zero;
- size_t rank1 = 0;
+ // size_t rank1 = 0;
size_t rank2 = 1;
for (size_t j = 0; j < 3; ++j) {
first[j] = eigenvectors(0, j);
@@ -308,7 +308,7 @@ EigenvalueSolver::findEigenvector(const TinyMatrix<3, 3>& A,
if (l2Norm(first) < 1e-3) {
for (size_t j = 0; j < 3; ++j) {
first[j] = eigenvectors(2, j);
- rank1 = 2;
+ // rank1 = 2;
}
}
if (l2Norm(second) < 1e-3) {
@@ -324,9 +324,9 @@ EigenvalueSolver::findEigenvector(const TinyMatrix<3, 3>& A,
for (size_t j = 0; j < 3; ++j) {
eigenvectors(rank2, j) = normal_vector[j];
}
- std::cout << " rank1 : " << rank1 << " rank2 " << rank2 << "\n";
+ // std::cout << " rank1 : " << rank1 << " rank2 " << rank2 << "\n";
}
- std::cout << "In findEigenvectors for eigenvalue " << eigenvalue << " eigenvectors " << eigenvectors << "\n";
+ // std::cout << "In findEigenvectors for eigenvalue " << eigenvalue << " eigenvectors " << eigenvectors << "\n";
return eigenvectors;
}
@@ -361,22 +361,21 @@ EigenvalueSolver::findEigenvectors(const TinyMatrix<3, 3>& A, TinyVector<3> eige
}
}
}
- std::cout << "space_size: " << space_size << "\n";
- std::cout << "eigenvalues: " << eigenvalues << "\n";
+ // std::cout << "space_size: " << space_size << "\n";
+ // std::cout << "eigenvalues: " << eigenvalues << "\n";
size_t it = 0;
while (count < 3 && it < 3) {
TinyMatrix<3> Try = findEigenvector(A, eigenvalues[count], space_size[count], epsilon);
- std::cout << eigenvalues[count] << " Frobenius norm try: " << frobeniusNorm(Try) << "\n";
- if (frobeniusNorm(Try) < 1e-3) {
- std::cout << " Error : no eigenvector corresponds to eigenvalue " << eigenvalues[count] << "\n";
- }
+ // std::cout << eigenvalues[count] << " Frobenius norm try: " << frobeniusNorm(Try) << "\n";
+ Assert(frobeniusNorm(Try) > 1e-3, " Error : no eigenvector corresponds to eigenvalue ");
+
TinyVector<3> eigenvector = zero;
int count2 = 0;
for (size_t i = 0; i < 3; ++i) {
for (size_t j = 0; j < 3; j++) {
eigenvector[j] = Try(i, j);
}
- std::cout << " count, i: " << count << ", " << i << " l2Norm(eigenvector) " << l2Norm(eigenvector) << "\n";
+ // std::cout << " count, i: " << count << ", " << i << " l2Norm(eigenvector) " << l2Norm(eigenvector) << "\n";
if (l2Norm(eigenvector) > 1e-3) {
for (size_t j = 0; j < 3; j++) {
eigenMatrix(count, j) = eigenvector[j];
@@ -385,13 +384,10 @@ EigenvalueSolver::findEigenvectors(const TinyMatrix<3, 3>& A, TinyVector<3> eige
count2 += 1;
}
}
- if (count2 != space_size[count - count2]) {
- std::cout << " error eigenvalue " << eigenvalues[count - 1] << " eigenvector space size " << count2
- << " is not what was expected " << space_size[count - 1] << "\n";
- }
+ Assert(count2 == space_size[count - count2], " eigenvector space size is not what was expected ");
it++;
}
- std ::cout << " eigenMatrix " << eigenMatrix << "\n";
+ // std ::cout << " eigenMatrix " << eigenMatrix << "\n";
return eigenMatrix;
}
@@ -400,22 +396,22 @@ EigenvalueSolver::findEigen(const TinyMatrix<3> A)
{
const double epsilon = 1.e-6; // * l2Norm(eigenvalues);
if (frobeniusNorm(A) < 1.e-15) {
- std::cout << " Frobenius norm 0 " << frobeniusNorm(A) << "\n";
+ // std::cout << " Frobenius norm 0 " << frobeniusNorm(A) << "\n";
return std::make_tuple(eigenvalues0, eigenvectors0);
}
TinyMatrix<3> C = 1. / frobeniusNorm(A) * A;
if (isDiagonal(C, epsilon)) {
- std::cout << "Matrix C " << C << " is diagonal "
- << "\n";
+ // std::cout << "Matrix C " << C << " is diagonal "
+ // << "\n";
const TinyVector<3> eigenvalues(A(0, 0), A(1, 1), A(2, 2));
TinyMatrix<3> Diag = zero;
for (size_t i = 0; i < 3; ++i) {
Diag(i, i) = eigenvalues[i];
}
TinyMatrix<3> B = eigenvectors0 * Diag * transpose(eigenvectors0);
- std::cout << "\n"
- << B << ", " << A << " Diff "
- << " A - B " << 1 / frobeniusNorm(A) * (A - B) << "\n";
+ // std::cout << "\n"
+ // << B << ", " << A << " Diff "
+ // << " A - B " << 1 / frobeniusNorm(A) * (A - B) << "\n";
return std::make_tuple(eigenvalues, eigenvectors0);
}
const TinyVector<3> eigenvalues = _findEigenValues(C, epsilon);
@@ -426,21 +422,21 @@ EigenvalueSolver::findEigen(const TinyMatrix<3> A)
Diag(i, i) = frobeniusNorm(A) * eigenvalues[i];
}
TinyMatrix<3> B = Eigenmatrix * Diag * transpose(Eigenmatrix);
- std::cout << "\n"
- << B << ", " << A << " Diff "
- << " A - B " << 1 / frobeniusNorm(A) * (A - B) << "\n";
+ // std::cout << "\n"
+ // << B << ", " << A << " Diff "
+ // << " A - B " << 1 / frobeniusNorm(A) * (A - B) << "\n";
return std::make_tuple(frobeniusNorm(A) * eigenvalues, Eigenmatrix);
}
TinyVector<3>
EigenvalueSolver::_findEigenValues(const TinyMatrix<3>& A, const double epsilon) const
{
- std::cout << " Matrix " << A << "\n";
+ // std::cout << " Matrix " << A << "\n";
TinyVector<3> eigenvalues(0., 0., 0.);
double b = -trace(A);
double c = 0.5 * (trace(A) * trace(A) - trace(A * A));
double d = -det(A);
Polynomial<3> P(d, c, b, 1.);
- std::cout << P << "\n";
+ // std::cout << P << "\n";
Polynomial<2> Q(c, b, 1.);
if (fabs(d) > 1e-11) {
eigenvalues[0] = _findFirstRoot(P);
@@ -448,9 +444,9 @@ EigenvalueSolver::_findEigenValues(const TinyMatrix<3>& A, const double epsilon)
Polynomial<3> S;
Polynomial<3> R;
divide(P, Q1, S, R);
- std::cout << "Q1 : " << Q1 << "\n";
- std::cout << "R : " << R << "\n";
- std::cout << "S : " << S << "\n";
+ // std::cout << "Q1 : " << Q1 << "\n";
+ // std::cout << "R : " << R << "\n";
+ // std::cout << "S : " << S << "\n";
Q.coefficients()[0] = S.coefficients()[0];
Q.coefficients()[1] = S.coefficients()[1];
Q.coefficients()[2] = S.coefficients()[2];
@@ -479,7 +475,7 @@ EigenvalueSolver::_findFirstRoot(Polynomial<3> P) const
{
double guess = -P.coefficients()[2] / 3.;
double old_guess = -P.coefficients()[2] / 3.;
- std::cout << " coefs P " << P.coefficients() << "\n";
+ // std::cout << " coefs P " << P.coefficients() << "\n";
// double delta = pow(2.*P.coefficients()[2],2) - 4*3.*P.coefficients()[3]*P.coefficients()[1];
Polynomial<2> Q = derivative(P);
Polynomial<1> R = derivative(Q);
@@ -489,13 +485,13 @@ EigenvalueSolver::_findFirstRoot(Polynomial<3> P) const
// guess -= P(guess) / Q(guess);
old_guess = guess;
guess -= 2 * P(guess) * Q(guess) / (2 * Q(guess) * Q(guess) - P(guess) * R(guess));
- std::cout << "guess = " << guess << "\n";
- std::cout << "g(guess) = " << Q(guess) << "\n";
+ // std::cout << "guess = " << guess << "\n";
+ // std::cout << "g(guess) = " << Q(guess) << "\n";
residu = P(guess);
- std::cout << "residu = " << residu << "\n";
+ // std::cout << "residu = " << residu << "\n";
iteration += 1;
} while ((fabs(residu) > 1.e-16) or (fabs(old_guess - guess) > 1.e-10));
- std::cout << " nb Newton iterations " << iteration << "\n";
+ // std::cout << " nb Newton iterations " << iteration << "\n";
return guess;
}
--
GitLab