Rasmus Munk LarsenandRasmus Munk Larsen
d6dfad0329
libeigen/eigen!2568 Closes #649 Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>
104 lines
3.1 KiB
C++
104 lines
3.1 KiB
C++
// This file is part of Eigen, a lightweight C++ template library
|
|
// for linear algebra.
|
|
//
|
|
// Copyright (C) 2026 Rasmus Munk Larsen <rmlarsen@gmail.com>
|
|
//
|
|
// This Source Code Form is subject to the terms of the Mozilla
|
|
// Public License v. 2.0. If a copy of the MPL was not distributed
|
|
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
|
|
// SPDX-License-Identifier: MPL-2.0
|
|
|
|
#include "sparse.h"
|
|
#include <Eigen/OrderingMethods>
|
|
|
|
template <typename Idx>
|
|
void verify_is_permutation(const PermutationMatrix<Dynamic, Dynamic, Idx>& perm) {
|
|
typedef Matrix<Idx, Dynamic, 1> Vec;
|
|
Vec sorted = perm.indices();
|
|
std::sort(sorted.data(), sorted.data() + sorted.size());
|
|
Vec expected = Vec::LinSpaced(perm.size(), 0, Idx(perm.size() - 1));
|
|
VERIFY_IS_EQUAL(sorted, expected);
|
|
}
|
|
|
|
template <typename Scalar>
|
|
void test_colamd_uncompressed() {
|
|
typedef int Idx;
|
|
typedef SparseMatrix<Scalar, ColMajor, Idx> SpMat;
|
|
const int n = 32;
|
|
SpMat A(n, n);
|
|
A.reserve(VectorXi::Constant(n, 5));
|
|
for (int j = 0; j < n; ++j) {
|
|
A.insert(j, j) = Scalar(2) + Scalar(j);
|
|
if (j > 0) A.insert(j - 1, j) = Scalar(-1);
|
|
if (j + 1 < n) A.insert(j + 1, j) = Scalar(-1);
|
|
}
|
|
A.insert(0, n - 1) = Scalar(0.5);
|
|
A.insert(n - 1, 0) = Scalar(0.25);
|
|
VERIFY(!A.isCompressed());
|
|
|
|
// Regression for #649: COLAMD on uncompressed input used to silently corrupt
|
|
// the permutation in release builds.
|
|
PermutationMatrix<Dynamic, Dynamic, Idx> perm_u, perm_c;
|
|
COLAMDOrdering<Idx> ord;
|
|
ord(A, perm_u);
|
|
verify_is_permutation(perm_u);
|
|
|
|
SpMat B = A;
|
|
B.makeCompressed();
|
|
ord(B, perm_c);
|
|
verify_is_permutation(perm_c);
|
|
|
|
VERIFY_IS_EQUAL(perm_u.indices(), perm_c.indices());
|
|
}
|
|
|
|
void test_colamd_mixed_index_types() {
|
|
typedef SparseMatrix<double, ColMajor, int> SpMat;
|
|
const int n = 8;
|
|
SpMat A(n, n);
|
|
A.reserve(VectorXi::Constant(n, 3));
|
|
for (int j = 0; j < n; ++j) {
|
|
A.insert(j, j) = 2.0 + j;
|
|
if (j > 0) A.insert(j - 1, j) = -1.0;
|
|
if (j + 1 < n) A.insert(j + 1, j) = -1.0;
|
|
}
|
|
A.makeCompressed();
|
|
|
|
// COLAMDOrdering<StorageIndex> is a public API choice independent of the
|
|
// source matrix's StorageIndex. Keep accepting mixed index types.
|
|
PermutationMatrix<Dynamic, Dynamic, long> perm;
|
|
COLAMDOrdering<long> ord;
|
|
ord(A, perm);
|
|
verify_is_permutation(perm);
|
|
}
|
|
|
|
template <typename Scalar>
|
|
void test_amd_uncompressed() {
|
|
typedef int Idx;
|
|
typedef SparseMatrix<Scalar, ColMajor, Idx> SpMat;
|
|
const int n = 24;
|
|
SpMat A(n, n);
|
|
A.reserve(VectorXi::Constant(n, 4));
|
|
for (int j = 0; j < n; ++j) {
|
|
A.insert(j, j) = Scalar(2) + Scalar(j);
|
|
if (j > 0) {
|
|
A.insert(j - 1, j) = Scalar(-1);
|
|
A.insert(j, j - 1) = Scalar(-1);
|
|
}
|
|
}
|
|
VERIFY(!A.isCompressed());
|
|
|
|
PermutationMatrix<Dynamic, Dynamic, Idx> perm;
|
|
AMDOrdering<Idx> ord;
|
|
ord(A, perm);
|
|
verify_is_permutation(perm);
|
|
}
|
|
|
|
EIGEN_DECLARE_TEST(sparse_ordering) {
|
|
for (int i = 0; i < g_repeat; ++i) {
|
|
CALL_SUBTEST_1(test_colamd_uncompressed<double>());
|
|
CALL_SUBTEST_2(test_colamd_uncompressed<std::complex<double> >());
|
|
CALL_SUBTEST_3(test_amd_uncompressed<double>());
|
|
CALL_SUBTEST_4(test_colamd_mixed_index_types());
|
|
}
|
|
}
|