Skip to content

Reuse matrices during iteration #102

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
maekclena wants to merge 1 commit into
base: develop
Choose a base branch
from
+92 −86
Open

Reuse matrices during iteration #102

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
34 changes: 18 additions & 16 deletions include/OpenABF/ABF.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -307,23 +307,30 @@ class ABF
}
}

// Typedefs
using Triplet = Eigen::Triplet<T>;
using SparseMatrix = Eigen::SparseMatrix<T>;
using DenseVector = Eigen::Matrix<T, Eigen::Dynamic, 1>;

// Helpful parameters
auto vIntCnt = mesh->num_vertices_interior();
auto edgeCnt = mesh->num_edges();
auto faceCnt = mesh->num_faces();

std::vector<Triplet> 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<T>;
using SparseMatrix = Eigen::SparseMatrix<T>;
using DenseVector = Eigen::Matrix<T, Eigen::Dynamic, 1>;

// Helpful parameters
auto vIntCnt = mesh->num_vertices_interior();
auto edgeCnt = mesh->num_edges();
auto faceCnt = mesh->num_faces();

//// RHS ////
// b1 = -alpha gradient
std::vector<Triplet> triplets;
triplets.clear();
std::size_t idx{0};
for (const auto& e : mesh->edges()) {
triplets.emplace_back(idx, 0, -AlphaGrad<T>(e));
Expand All @@ -342,8 +349,6 @@ class ABF
triplets.emplace_back(vIntCnt + idx, 0, -LenGrad<T>(v));
++idx;
}
SparseMatrix b(edgeCnt + faceCnt + 2 * vIntCnt, 1);
b.reserve(triplets.size());
b.setFromTriplets(triplets.begin(), triplets.end());

///// LHS /////
Expand Down Expand Up @@ -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();
Expand All @@ -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");
}
Expand Down
52 changes: 28 additions & 24 deletions include/OpenABF/ABFPlusPlus.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -120,29 +120,39 @@ class ABFPlusPlus
}
}

// Typedefs
using Triplet = Eigen::Triplet<T>;
using SparseMatrix = Eigen::SparseMatrix<T>;
using DenseVector = Eigen::Matrix<T, Eigen::Dynamic, 1>;

// Helpful parameters
auto vIntCnt = mesh->num_vertices_interior();
auto edgeCnt = mesh->num_edges();
auto faceCnt = mesh->num_faces();

std::vector<Triplet> 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<T>;
using SparseMatrix = Eigen::SparseMatrix<T>;
using DenseVector = Eigen::Matrix<T, Eigen::Dynamic, 1>;

// Helpful parameters
auto vIntCnt = mesh->num_vertices_interior();
auto edgeCnt = mesh->num_edges();
auto faceCnt = mesh->num_faces();

// b1 = -alpha gradient
std::vector<Triplet> triplets;
triplets.clear();
std::size_t idx{0};
for (const auto& e : mesh->edges()) {
triplets.emplace_back(idx, 0, -AlphaGrad<T>(e));
++idx;
}
SparseMatrix b1(edgeCnt, 1);
b1.reserve(triplets.size());
b1.setFromTriplets(triplets.begin(), triplets.end());

// b2 = -lambda gradient
Expand All @@ -159,8 +169,6 @@ class ABFPlusPlus
triplets.emplace_back(vIntCnt + idx, 0, -LenGrad<T>(v));
idx++;
}
SparseMatrix b2(faceCnt + 2 * vIntCnt, 1);
b2.reserve(triplets.size());
b2.setFromTriplets(triplets.begin(), triplets.end());

// Compute J1 + J2
Expand All @@ -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)
Expand All @@ -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();
Expand All @@ -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);
Expand All @@ -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()) {
Expand Down
3 changes: 0 additions & 3 deletions include/OpenABF/detail/LSCMSystem.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -187,7 +187,6 @@ auto BuildSystem(const typename MeshType::Pointer& mesh, const PinMap<T>& 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.
Expand Down Expand Up @@ -269,11 +268,9 @@ auto BuildSystem(const typename MeshType::Pointer& mesh, const PinMap<T>& 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);
Expand Down
Loading