From fe4ec7d56293763c94d555803acf2267379a2e91 Mon Sep 17 00:00:00 2001 From: maekclena Date: Wed, 8 Jul 2026 20:19:33 +0200 Subject: [PATCH] perf: keep Eigen matrices to reuse memory --- include/OpenABF/ABF.hpp | 34 +++++----- include/OpenABF/ABFPlusPlus.hpp | 52 ++++++++-------- include/OpenABF/detail/LSCMSystem.hpp | 3 - single_include/OpenABF/OpenABF.hpp | 89 ++++++++++++++------------- 4 files changed, 92 insertions(+), 86 deletions(-) diff --git a/include/OpenABF/ABF.hpp b/include/OpenABF/ABF.hpp index 98f5e27..6ad26f8 100644 --- a/include/OpenABF/ABF.hpp +++ b/include/OpenABF/ABF.hpp @@ -307,23 +307,30 @@ class ABF } } + // Typedefs + using Triplet = Eigen::Triplet; + using SparseMatrix = Eigen::SparseMatrix; + using DenseVector = Eigen::Matrix; + + // Helpful parameters + auto vIntCnt = mesh->num_vertices_interior(); + auto edgeCnt = mesh->num_edges(); + auto faceCnt = mesh->num_faces(); + + std::vector triplets; + auto Asize = edgeCnt + faceCnt + 2 * vIntCnt; + SparseMatrix b(Asize, 1); + SparseMatrix A(Asize, Asize); + DenseVector delta(Asize, 1); + while (gradient > gradThreshold and gradDelta > gradThreshold and iters < maxIters) { if (std::isnan(gradient) or std::isinf(gradient)) { throw MeshException("Mesh gradient cannot be computed"); } - // Typedefs - using Triplet = Eigen::Triplet; - using SparseMatrix = Eigen::SparseMatrix; - using DenseVector = Eigen::Matrix; - - // Helpful parameters - auto vIntCnt = mesh->num_vertices_interior(); - auto edgeCnt = mesh->num_edges(); - auto faceCnt = mesh->num_faces(); //// RHS //// // b1 = -alpha gradient - std::vector triplets; + triplets.clear(); std::size_t idx{0}; for (const auto& e : mesh->edges()) { triplets.emplace_back(idx, 0, -AlphaGrad(e)); @@ -342,8 +349,6 @@ class ABF triplets.emplace_back(vIntCnt + idx, 0, -LenGrad(v)); ++idx; } - SparseMatrix b(edgeCnt + faceCnt + 2 * vIntCnt, 1); - b.reserve(triplets.size()); b.setFromTriplets(triplets.begin(), triplets.end()); ///// LHS ///// @@ -390,9 +395,6 @@ class ABF } ++idx; } - auto Asize = edgeCnt + faceCnt + 2 * vIntCnt; - SparseMatrix A(Asize, Asize); - A.reserve(triplets.size()); A.setFromTriplets(triplets.begin(), triplets.end()); A.makeCompressed(); @@ -401,7 +403,7 @@ class ABF if (solver.info() != Eigen::ComputationInfo::Success) { throw SolverException("ABF: Failed to solve A"); } - DenseVector delta = solver.solve(b); + delta = solver.solve(b); if (solver.info() != Eigen::ComputationInfo::Success) { throw SolverException("ABF: Failed to solve b"); } diff --git a/include/OpenABF/ABFPlusPlus.hpp b/include/OpenABF/ABFPlusPlus.hpp index 6999d9a..1612cc2 100644 --- a/include/OpenABF/ABFPlusPlus.hpp +++ b/include/OpenABF/ABFPlusPlus.hpp @@ -120,29 +120,39 @@ class ABFPlusPlus } } + // Typedefs + using Triplet = Eigen::Triplet; + using SparseMatrix = Eigen::SparseMatrix; + using DenseVector = Eigen::Matrix; + + // Helpful parameters + auto vIntCnt = mesh->num_vertices_interior(); + auto edgeCnt = mesh->num_edges(); + auto faceCnt = mesh->num_faces(); + + std::vector triplets; + SparseMatrix b1(edgeCnt, 1); + SparseMatrix b2(faceCnt + 2 * vIntCnt, 1); + SparseMatrix J(faceCnt + 2 * vIntCnt, 3 * faceCnt); + SparseMatrix LambdaInv(edgeCnt, edgeCnt); + SparseMatrix LambdaStarInv(faceCnt, faceCnt); + SparseMatrix A(2 * vIntCnt, 2 * vIntCnt); + SparseMatrix b(2 * vIntCnt, 1); + DenseVector deltaLambda(faceCnt + 2 * vIntCnt, 1); + DenseVector deltaAlpha(edgeCnt, 1); + while (gradient > gradThreshold and gradDelta > gradThreshold and iters < maxIters) { if (std::isnan(gradient) or std::isinf(gradient)) { throw MeshException("Mesh gradient cannot be computed"); } - // Typedefs - using Triplet = Eigen::Triplet; - using SparseMatrix = Eigen::SparseMatrix; - using DenseVector = Eigen::Matrix; - - // Helpful parameters - auto vIntCnt = mesh->num_vertices_interior(); - auto edgeCnt = mesh->num_edges(); - auto faceCnt = mesh->num_faces(); // b1 = -alpha gradient - std::vector triplets; + triplets.clear(); std::size_t idx{0}; for (const auto& e : mesh->edges()) { triplets.emplace_back(idx, 0, -AlphaGrad(e)); ++idx; } - SparseMatrix b1(edgeCnt, 1); - b1.reserve(triplets.size()); b1.setFromTriplets(triplets.begin(), triplets.end()); // b2 = -lambda gradient @@ -159,8 +169,6 @@ class ABFPlusPlus triplets.emplace_back(vIntCnt + idx, 0, -LenGrad(v)); idx++; } - SparseMatrix b2(faceCnt + 2 * vIntCnt, 1); - b2.reserve(triplets.size()); b2.setFromTriplets(triplets.begin(), triplets.end()); // Compute J1 + J2 @@ -187,8 +195,6 @@ class ABFPlusPlus } ++idx; } - SparseMatrix J(faceCnt + 2 * vIntCnt, 3 * faceCnt); - J.reserve(triplets.size()); J.setFromTriplets(triplets.begin(), triplets.end()); // Lambda = diag(2/w) @@ -200,15 +206,13 @@ class ABFPlusPlus triplets.emplace_back(idx, idx, T(1) / (2 * e->weight)); ++idx; } - SparseMatrix LambdaInv(edgeCnt, edgeCnt); - LambdaInv.reserve(edgeCnt); LambdaInv.setFromTriplets(triplets.begin(), triplets.end()); // solve Eq. 16 auto bstar = J * LambdaInv * b1 - b2; auto JLiJt = J * LambdaInv * J.transpose(); - SparseMatrix LambdaStarInv = JLiJt.block(0, 0, faceCnt, faceCnt); + LambdaStarInv = JLiJt.block(0, 0, faceCnt, faceCnt); for (int k = 0; k < LambdaStarInv.outerSize(); ++k) { for (typename SparseMatrix::InnerIterator it(LambdaStarInv, k); it; ++it) { it.valueRef() = T(1) / it.value(); @@ -221,8 +225,8 @@ class ABFPlusPlus auto bstar2 = bstar.block(faceCnt, 0, 2 * vIntCnt, 1); // (J* Lam*^-1 J*^t - J**) delta_lambda_2 = J* Lam*^-1 b*_1 - b*_2 - SparseMatrix A = Jstar * LambdaStarInv * JstarT - Jstar2; - SparseMatrix b = Jstar * LambdaStarInv * bstar1 - bstar2; + A = Jstar * LambdaStarInv * JstarT - Jstar2; + b = Jstar * LambdaStarInv * bstar1 - bstar2; A.makeCompressed(); Solver solver; solver.compute(A); @@ -238,11 +242,11 @@ class ABFPlusPlus auto deltaLambda1 = LambdaStarInv * (bstar1 - JstarT * deltaLambda2); // Construct deltaLambda - DenseVector deltaLambda(deltaLambda1.rows() + deltaLambda2.rows(), 1); - deltaLambda << DenseVector(deltaLambda1), DenseVector(deltaLambda2); + deltaLambda.topRows(faceCnt) = deltaLambda1; + deltaLambda.bottomRows(2 * vIntCnt) = deltaLambda2; // Compute Eq. 10 -> delta_alpha - DenseVector deltaAlpha = LambdaInv * (b1 - J.transpose() * deltaLambda); + deltaAlpha = LambdaInv * (b1 - J.transpose() * deltaLambda); // lambda += delta_lambda for (auto& f : mesh->faces()) { diff --git a/include/OpenABF/detail/LSCMSystem.hpp b/include/OpenABF/detail/LSCMSystem.hpp index cb6188a..b66d7dc 100644 --- a/include/OpenABF/detail/LSCMSystem.hpp +++ b/include/OpenABF/detail/LSCMSystem.hpp @@ -187,7 +187,6 @@ auto BuildSystem(const typename MeshType::Pointer& mesh, const PinMap& pins) tripletsB.emplace_back(2 * s + 1, 0, uv[1]); } SparseMatrix bFixed(2 * numFixed, 1); - bFixed.reserve(tripletsB.size()); bFixed.setFromTriplets(tripletsB.begin(), tripletsB.end()); // Permutation for free vertices: maps mesh vertex idx → row-pair slot in A. @@ -269,11 +268,9 @@ auto BuildSystem(const typename MeshType::Pointer& mesh, const PinMap& pins) } SparseMatrix A(2 * numFaces, 2 * numFree); - A.reserve(tripletsA.size()); A.setFromTriplets(tripletsA.begin(), tripletsA.end()); SparseMatrix bFree(2 * numFaces, 2 * numFixed); - bFree.reserve(tripletsB.size()); bFree.setFromTriplets(tripletsB.begin(), tripletsB.end()); SparseMatrix b = bFree * bFixed * T(-1); diff --git a/single_include/OpenABF/OpenABF.hpp b/single_include/OpenABF/OpenABF.hpp index b21d488..15755f1 100644 --- a/single_include/OpenABF/OpenABF.hpp +++ b/single_include/OpenABF/OpenABF.hpp @@ -2615,23 +2615,30 @@ class ABF } } + // Typedefs + using Triplet = Eigen::Triplet; + using SparseMatrix = Eigen::SparseMatrix; + using DenseVector = Eigen::Matrix; + + // Helpful parameters + auto vIntCnt = mesh->num_vertices_interior(); + auto edgeCnt = mesh->num_edges(); + auto faceCnt = mesh->num_faces(); + + std::vector triplets; + auto Asize = edgeCnt + faceCnt + 2 * vIntCnt; + SparseMatrix b(Asize, 1); + SparseMatrix A(Asize, Asize); + DenseVector delta(Asize, 1); + while (gradient > gradThreshold and gradDelta > gradThreshold and iters < maxIters) { if (std::isnan(gradient) or std::isinf(gradient)) { throw MeshException("Mesh gradient cannot be computed"); } - // Typedefs - using Triplet = Eigen::Triplet; - using SparseMatrix = Eigen::SparseMatrix; - using DenseVector = Eigen::Matrix; - - // Helpful parameters - auto vIntCnt = mesh->num_vertices_interior(); - auto edgeCnt = mesh->num_edges(); - auto faceCnt = mesh->num_faces(); //// RHS //// // b1 = -alpha gradient - std::vector triplets; + triplets.clear(); std::size_t idx{0}; for (const auto& e : mesh->edges()) { triplets.emplace_back(idx, 0, -AlphaGrad(e)); @@ -2650,8 +2657,6 @@ class ABF triplets.emplace_back(vIntCnt + idx, 0, -LenGrad(v)); ++idx; } - SparseMatrix b(edgeCnt + faceCnt + 2 * vIntCnt, 1); - b.reserve(triplets.size()); b.setFromTriplets(triplets.begin(), triplets.end()); ///// LHS ///// @@ -2698,9 +2703,6 @@ class ABF } ++idx; } - auto Asize = edgeCnt + faceCnt + 2 * vIntCnt; - SparseMatrix A(Asize, Asize); - A.reserve(triplets.size()); A.setFromTriplets(triplets.begin(), triplets.end()); A.makeCompressed(); @@ -2709,7 +2711,7 @@ class ABF if (solver.info() != Eigen::ComputationInfo::Success) { throw SolverException("ABF: Failed to solve A"); } - DenseVector delta = solver.solve(b); + delta = solver.solve(b); if (solver.info() != Eigen::ComputationInfo::Success) { throw SolverException("ABF: Failed to solve b"); } @@ -2891,29 +2893,39 @@ class ABFPlusPlus } } + // Typedefs + using Triplet = Eigen::Triplet; + using SparseMatrix = Eigen::SparseMatrix; + using DenseVector = Eigen::Matrix; + + // Helpful parameters + auto vIntCnt = mesh->num_vertices_interior(); + auto edgeCnt = mesh->num_edges(); + auto faceCnt = mesh->num_faces(); + + std::vector triplets; + SparseMatrix b1(edgeCnt, 1); + SparseMatrix b2(faceCnt + 2 * vIntCnt, 1); + SparseMatrix J(faceCnt + 2 * vIntCnt, 3 * faceCnt); + SparseMatrix LambdaInv(edgeCnt, edgeCnt); + SparseMatrix LambdaStarInv(faceCnt, faceCnt); + SparseMatrix A(2 * vIntCnt, 2 * vIntCnt); + SparseMatrix b(2 * vIntCnt, 1); + DenseVector deltaLambda(faceCnt + 2 * vIntCnt, 1); + DenseVector deltaAlpha(edgeCnt, 1); + while (gradient > gradThreshold and gradDelta > gradThreshold and iters < maxIters) { if (std::isnan(gradient) or std::isinf(gradient)) { throw MeshException("Mesh gradient cannot be computed"); } - // Typedefs - using Triplet = Eigen::Triplet; - using SparseMatrix = Eigen::SparseMatrix; - using DenseVector = Eigen::Matrix; - - // Helpful parameters - auto vIntCnt = mesh->num_vertices_interior(); - auto edgeCnt = mesh->num_edges(); - auto faceCnt = mesh->num_faces(); // b1 = -alpha gradient - std::vector triplets; + triplets.clear(); std::size_t idx{0}; for (const auto& e : mesh->edges()) { triplets.emplace_back(idx, 0, -AlphaGrad(e)); ++idx; } - SparseMatrix b1(edgeCnt, 1); - b1.reserve(triplets.size()); b1.setFromTriplets(triplets.begin(), triplets.end()); // b2 = -lambda gradient @@ -2930,8 +2942,6 @@ class ABFPlusPlus triplets.emplace_back(vIntCnt + idx, 0, -LenGrad(v)); idx++; } - SparseMatrix b2(faceCnt + 2 * vIntCnt, 1); - b2.reserve(triplets.size()); b2.setFromTriplets(triplets.begin(), triplets.end()); // Compute J1 + J2 @@ -2958,8 +2968,6 @@ class ABFPlusPlus } ++idx; } - SparseMatrix J(faceCnt + 2 * vIntCnt, 3 * faceCnt); - J.reserve(triplets.size()); J.setFromTriplets(triplets.begin(), triplets.end()); // Lambda = diag(2/w) @@ -2971,15 +2979,13 @@ class ABFPlusPlus triplets.emplace_back(idx, idx, T(1) / (2 * e->weight)); ++idx; } - SparseMatrix LambdaInv(edgeCnt, edgeCnt); - LambdaInv.reserve(edgeCnt); LambdaInv.setFromTriplets(triplets.begin(), triplets.end()); // solve Eq. 16 auto bstar = J * LambdaInv * b1 - b2; auto JLiJt = J * LambdaInv * J.transpose(); - SparseMatrix LambdaStarInv = JLiJt.block(0, 0, faceCnt, faceCnt); + LambdaStarInv = JLiJt.block(0, 0, faceCnt, faceCnt); for (int k = 0; k < LambdaStarInv.outerSize(); ++k) { for (typename SparseMatrix::InnerIterator it(LambdaStarInv, k); it; ++it) { it.valueRef() = T(1) / it.value(); @@ -2992,8 +2998,8 @@ class ABFPlusPlus auto bstar2 = bstar.block(faceCnt, 0, 2 * vIntCnt, 1); // (J* Lam*^-1 J*^t - J**) delta_lambda_2 = J* Lam*^-1 b*_1 - b*_2 - SparseMatrix A = Jstar * LambdaStarInv * JstarT - Jstar2; - SparseMatrix b = Jstar * LambdaStarInv * bstar1 - bstar2; + A = Jstar * LambdaStarInv * JstarT - Jstar2; + b = Jstar * LambdaStarInv * bstar1 - bstar2; A.makeCompressed(); Solver solver; solver.compute(A); @@ -3009,11 +3015,11 @@ class ABFPlusPlus auto deltaLambda1 = LambdaStarInv * (bstar1 - JstarT * deltaLambda2); // Construct deltaLambda - DenseVector deltaLambda(deltaLambda1.rows() + deltaLambda2.rows(), 1); - deltaLambda << DenseVector(deltaLambda1), DenseVector(deltaLambda2); + deltaLambda.topRows(faceCnt) = deltaLambda1; + deltaLambda.bottomRows(2 * vIntCnt) = deltaLambda2; // Compute Eq. 10 -> delta_alpha - DenseVector deltaAlpha = LambdaInv * (b1 - J.transpose() * deltaLambda); + deltaAlpha = LambdaInv * (b1 - J.transpose() * deltaLambda); // lambda += delta_lambda for (auto& f : mesh->faces()) { @@ -3280,7 +3286,6 @@ auto BuildSystem(const typename MeshType::Pointer& mesh, const PinMap& pins) tripletsB.emplace_back(2 * s + 1, 0, uv[1]); } SparseMatrix bFixed(2 * numFixed, 1); - bFixed.reserve(tripletsB.size()); bFixed.setFromTriplets(tripletsB.begin(), tripletsB.end()); // Permutation for free vertices: maps mesh vertex idx → row-pair slot in A. @@ -3362,11 +3367,9 @@ auto BuildSystem(const typename MeshType::Pointer& mesh, const PinMap& pins) } SparseMatrix A(2 * numFaces, 2 * numFree); - A.reserve(tripletsA.size()); A.setFromTriplets(tripletsA.begin(), tripletsA.end()); SparseMatrix bFree(2 * numFaces, 2 * numFixed); - bFree.reserve(tripletsB.size()); bFree.setFromTriplets(tripletsB.begin(), tripletsB.end()); SparseMatrix b = bFree * bFixed * T(-1);