961 lines
20 KiB
C++
961 lines
20 KiB
C++
// SPDX-License-Identifier: Apache-2.0
|
|
//
|
|
// Copyright 2011-2017 Ryan Curtin (http://www.ratml.org/)
|
|
// Copyright 2017 National ICT Australia (NICTA)
|
|
//
|
|
// Licensed under the Apache License, Version 2.0 (the "License");
|
|
// you may not use this file except in compliance with the License.
|
|
// You may obtain a copy of the License at
|
|
// http://www.apache.org/licenses/LICENSE-2.0
|
|
//
|
|
// Unless required by applicable law or agreed to in writing, software
|
|
// distributed under the License is distributed on an "AS IS" BASIS,
|
|
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
|
// See the License for the specific language governing permissions and
|
|
// limitations under the License.
|
|
// ------------------------------------------------------------------------
|
|
|
|
#include <armadillo>
|
|
|
|
#include "catch.hpp"
|
|
|
|
using namespace arma;
|
|
|
|
#if defined(ARMA_USE_SUPERLU)
|
|
|
|
TEST_CASE("fn_spsolve_sparse_test", "[spsolve]")
|
|
{
|
|
// We want to spsolve a system of equations, AX = B, where we want to recover
|
|
// X and we have A and B, and A is sparse.
|
|
for (size_t t = 0; t < 10; ++t)
|
|
{
|
|
const uword size = 5 * (t + 1);
|
|
|
|
mat rX;
|
|
rX.randu(size, size);
|
|
|
|
sp_mat A;
|
|
A.sprandu(size, size, 0.25);
|
|
A.diag().randu();
|
|
A.diag() += 1;
|
|
|
|
mat B = A * rX;
|
|
|
|
mat X;
|
|
bool result = spsolve(X, A, B);
|
|
REQUIRE( result );
|
|
|
|
// Dense solver.
|
|
mat dA(A);
|
|
mat dX = solve(dA, B);
|
|
|
|
REQUIRE( X.n_cols == dX.n_cols );
|
|
REQUIRE( X.n_rows == dX.n_rows );
|
|
|
|
for (uword i = 0; i < dX.n_cols; ++i)
|
|
{
|
|
for (uword j = 0; j < dX.n_rows; ++j)
|
|
{
|
|
REQUIRE( (double) X(j, i) == Approx((double) dX(j, i)).margin(0.01) );
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
TEST_CASE("fn_spsolve_sparse_nonsymmetric_test", "[spsolve]")
|
|
{
|
|
for (size_t t = 0; t < 10; ++t)
|
|
{
|
|
const uword r_size = 5 * (t + 1);
|
|
const uword c_size = 3 * (t + 4);
|
|
|
|
mat rX;
|
|
rX.randu(r_size, c_size);
|
|
|
|
sp_mat A;
|
|
A.sprandu(r_size, r_size, 0.25);
|
|
A.diag().randu();
|
|
A.diag() += 1;
|
|
|
|
mat B = A * rX;
|
|
|
|
mat X;
|
|
bool result = spsolve(X, A, B);
|
|
REQUIRE( result );
|
|
|
|
// Dense solver.
|
|
mat dA(A);
|
|
mat dX = solve(dA, B);
|
|
|
|
REQUIRE( X.n_cols == dX.n_cols );
|
|
REQUIRE( X.n_rows == dX.n_rows );
|
|
|
|
for (uword i = 0; i < dX.n_cols; ++i)
|
|
{
|
|
for (uword j = 0; j < dX.n_rows; ++j)
|
|
{
|
|
REQUIRE( (double) X(j, i) == Approx((double) dX(j, i)).margin(0.01) );
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
TEST_CASE("fn_spsolve_sparse_float_test", "[spsolve]")
|
|
{
|
|
// We want to spsolve a system of equations, AX = B, where we want to recover
|
|
// X and we have A and B, and A is sparse.
|
|
size_t successes = 0;
|
|
for (size_t t = 0; t < 10; ++t)
|
|
{
|
|
const uword size = 5 * (t + 1);
|
|
|
|
fmat rX;
|
|
rX.randu(size, size);
|
|
|
|
SpMat<float> A;
|
|
A.sprandu(size, size, 0.25);
|
|
A.diag().randu();
|
|
A.diag() += 1;
|
|
|
|
fmat B = A * rX;
|
|
|
|
fmat X;
|
|
bool result = spsolve(X, A, B);
|
|
REQUIRE( result );
|
|
|
|
// Dense solver.
|
|
fmat dA(A);
|
|
fmat dX = solve(dA, B);
|
|
|
|
REQUIRE( X.n_cols == dX.n_cols );
|
|
REQUIRE( X.n_rows == dX.n_rows );
|
|
|
|
if (any(abs(vectorise(X - dX)) > 0.04))
|
|
{
|
|
// Skip this trial, it did not turn out correctly.
|
|
continue;
|
|
}
|
|
|
|
for (size_t i = 0; i < dX.n_cols; ++i)
|
|
{
|
|
for (size_t j = 0; j < dX.n_rows; ++j)
|
|
{
|
|
REQUIRE( (float) X(j, i) == Approx((float) dX(j, i)).margin(0.04) );
|
|
}
|
|
}
|
|
|
|
++successes;
|
|
}
|
|
|
|
REQUIRE( successes > 8 );
|
|
}
|
|
|
|
|
|
|
|
TEST_CASE("fn_spsolve_sparse_nonsymmetric_float_test", "[spsolve]")
|
|
{
|
|
for (size_t t = 0; t < 10; ++t)
|
|
{
|
|
const uword r_size = 5 * (t + 1);
|
|
const uword c_size = 3 * (t + 4);
|
|
|
|
fmat rX;
|
|
rX.randu(r_size, c_size);
|
|
|
|
SpMat<float> A;
|
|
A.sprandu(r_size, r_size, 0.25);
|
|
A.diag().randu();
|
|
A.diag() += 1;
|
|
|
|
fmat B = A * rX;
|
|
|
|
fmat X;
|
|
bool result = spsolve(X, A, B);
|
|
REQUIRE( result );
|
|
|
|
// Dense solver.
|
|
fmat dA(A);
|
|
fmat dX = solve(dA, B);
|
|
|
|
REQUIRE( X.n_cols == dX.n_cols );
|
|
REQUIRE( X.n_rows == dX.n_rows );
|
|
|
|
for (uword i = 0; i < dX.n_cols; ++i)
|
|
{
|
|
for (uword j = 0; j < dX.n_rows; ++j)
|
|
{
|
|
REQUIRE( (float) X(j, i) == Approx((float) dX(j, i)).margin(0.04) );
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
TEST_CASE("fn_spsolve_sparse_complex_float_test", "[spsolve]")
|
|
{
|
|
// We want to spsolve a system of equations, AX = B, where we want to recover
|
|
// X and we have A and B, and A is sparse.
|
|
for (size_t t = 0; t < 10; ++t)
|
|
{
|
|
const uword size = 5 * (t + 1);
|
|
|
|
Mat<cx_float> rX;
|
|
rX.randu(size, size);
|
|
|
|
SpMat<cx_float> A;
|
|
A.sprandu(size, size, 0.25);
|
|
A.diag().randu();
|
|
A.diag() += 1;
|
|
|
|
Mat<cx_float> B = A * rX;
|
|
|
|
Mat<cx_float> X;
|
|
bool result = spsolve(X, A, B);
|
|
REQUIRE( result );
|
|
|
|
// Dense solver.
|
|
Mat<cx_float> dA(A);
|
|
Mat<cx_float> dX = solve(dA, B);
|
|
|
|
REQUIRE( X.n_cols == dX.n_cols );
|
|
REQUIRE( X.n_rows == dX.n_rows );
|
|
|
|
for (uword i = 0; i < dX.n_cols; ++i)
|
|
{
|
|
for (uword j = 0; j < dX.n_rows; ++j)
|
|
{
|
|
REQUIRE( (float) std::abs((cx_float) X(j, i)) ==
|
|
Approx((float) std::abs((cx_float) dX(j, i))).margin(0.04) );
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
TEST_CASE("fn_spsolve_sparse_nonsymmetric_complex_float_test", "[spsolve]")
|
|
{
|
|
for (size_t t = 0; t < 10; ++t)
|
|
{
|
|
const uword r_size = 5 * (t + 1);
|
|
const uword c_size = 3 * (t + 4);
|
|
|
|
Mat<cx_float> rX;
|
|
rX.randu(r_size, c_size);
|
|
|
|
SpMat<cx_float> A;
|
|
A.sprandu(r_size, r_size, 0.25);
|
|
A.diag().randu();
|
|
A.diag() += 1;
|
|
|
|
Mat<cx_float> B = A * rX;
|
|
|
|
Mat<cx_float> X;
|
|
bool result = spsolve(X, A, B);
|
|
REQUIRE( result );
|
|
|
|
// Dense solver.
|
|
Mat<cx_float> dA(A);
|
|
Mat<cx_float> dX = solve(dA, B);
|
|
|
|
REQUIRE( X.n_cols == dX.n_cols );
|
|
REQUIRE( X.n_rows == dX.n_rows );
|
|
|
|
for (uword i = 0; i < dX.n_cols; ++i)
|
|
{
|
|
for (uword j = 0; j < dX.n_rows; ++j)
|
|
{
|
|
REQUIRE( (float) std::abs((cx_float) X(j, i)) ==
|
|
Approx((float) std::abs((cx_float) dX(j, i))).margin(0.04) );
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
TEST_CASE("fn_spsolve_sparse_complex_test", "[spsolve]")
|
|
{
|
|
// We want to spsolve a system of equations, AX = B, where we want to recover
|
|
// X and we have A and B, and A is sparse.
|
|
for (size_t t = 0; t < 10; ++t)
|
|
{
|
|
const uword size = 5 * (t + 1);
|
|
|
|
Mat<cx_double> rX;
|
|
rX.randu(size, size);
|
|
|
|
SpMat<cx_double> A;
|
|
A.sprandu(size, size, 0.25);
|
|
A.diag().randu();
|
|
A.diag() += 1;
|
|
|
|
Mat<cx_double> B = A * rX;
|
|
|
|
Mat<cx_double> X;
|
|
bool result = spsolve(X, A, B);
|
|
REQUIRE( result );
|
|
|
|
// Dense solver.
|
|
Mat<cx_double> dA(A);
|
|
Mat<cx_double> dX = solve(dA, B);
|
|
|
|
REQUIRE( X.n_cols == dX.n_cols );
|
|
REQUIRE( X.n_rows == dX.n_rows );
|
|
|
|
for (uword i = 0; i < dX.n_cols; ++i)
|
|
{
|
|
for (uword j = 0; j < dX.n_rows; ++j)
|
|
{
|
|
REQUIRE( (double) std::abs((cx_double) X(j, i)) ==
|
|
Approx((double) std::abs((cx_double) dX(j, i))).margin(0.01) );
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
TEST_CASE("fn_spsolve_sparse_nonsymmetric_complex_test", "[spsolve]")
|
|
{
|
|
for (size_t t = 0; t < 10; ++t)
|
|
{
|
|
const uword r_size = 5 * (t + 1);
|
|
const uword c_size = 3 * (t + 4);
|
|
|
|
Mat<cx_double> rX;
|
|
rX.randu(r_size, c_size);
|
|
|
|
SpMat<cx_double> A;
|
|
A.sprandu(r_size, r_size, 0.25);
|
|
A.diag().randu();
|
|
A.diag() += 1;
|
|
|
|
Mat<cx_double> B = A * rX;
|
|
|
|
Mat<cx_double> X;
|
|
bool result = spsolve(X, A, B);
|
|
REQUIRE( result );
|
|
|
|
// Dense solver.
|
|
Mat<cx_double> dA(A);
|
|
Mat<cx_double> dX = solve(dA, B);
|
|
|
|
REQUIRE( X.n_cols == dX.n_cols );
|
|
REQUIRE( X.n_rows == dX.n_rows );
|
|
|
|
for (uword i = 0; i < dX.n_cols; ++i)
|
|
{
|
|
for (uword j = 0; j < dX.n_rows; ++j)
|
|
{
|
|
REQUIRE( (double) std::abs((cx_double) X(j, i)) ==
|
|
Approx((double) std::abs((cx_double) dX(j, i))).margin(0.01) );
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
TEST_CASE("fn_spsolve_delayed_sparse_test", "[spsolve]")
|
|
{
|
|
const uword size = 10;
|
|
|
|
mat rX;
|
|
rX.randu(size, size);
|
|
|
|
sp_mat A;
|
|
A.sprandu(size, size, 0.25);
|
|
A.diag().randu();
|
|
A.diag() += 1;
|
|
|
|
mat B = A * rX;
|
|
|
|
mat X;
|
|
bool result = spsolve(X, A, B);
|
|
REQUIRE( result );
|
|
|
|
mat dX = spsolve(A, B);
|
|
|
|
REQUIRE( X.n_cols == dX.n_cols );
|
|
REQUIRE( X.n_rows == dX.n_rows );
|
|
|
|
for (uword i = 0; i < dX.n_cols; ++i)
|
|
{
|
|
for (uword j = 0; j < dX.n_rows; ++j)
|
|
{
|
|
REQUIRE( (double) X(j, i) == Approx((double) dX(j, i)).margin(0.01) );
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
TEST_CASE("fn_spsolve_superlu_solve_test", "[spsolve]")
|
|
{
|
|
// Solve this matrix, as in the examples:
|
|
// [[19 0 21 21 0]
|
|
// [12 21 0 0 0]
|
|
// [ 0 12 16 0 0]
|
|
// [ 0 0 0 5 21]
|
|
// [12 12 0 0 18]]
|
|
sp_mat b(5, 5);
|
|
b(0, 0) = 19;
|
|
b(0, 2) = 21;
|
|
b(0, 3) = 21;
|
|
b(1, 0) = 12;
|
|
b(1, 1) = 21;
|
|
b(2, 1) = 12;
|
|
b(2, 2) = 16;
|
|
b(3, 3) = 5;
|
|
b(3, 4) = 21;
|
|
b(4, 0) = 12;
|
|
b(4, 1) = 12;
|
|
b(4, 4) = 18;
|
|
|
|
mat db(b);
|
|
|
|
sp_mat a;
|
|
a.eye(5, 5);
|
|
mat da(a);
|
|
|
|
mat x;
|
|
spsolve(x, a, db);
|
|
|
|
mat dx = solve(da, db);
|
|
|
|
for (uword i = 0; i < x.n_cols; ++i)
|
|
{
|
|
for (uword j = 0; j < x.n_rows; ++j)
|
|
{
|
|
REQUIRE( (double) x(j, i) == Approx(dx(j, i)).margin(0.01) );
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
TEST_CASE("fn_spsolve_random_superlu_solve_test", "[spsolve]")
|
|
{
|
|
// Try to solve some random systems.
|
|
const size_t iterations = 10;
|
|
for (size_t it = 0; it < iterations; ++it)
|
|
{
|
|
sp_mat a;
|
|
a.sprandu(50, 50, 0.3);
|
|
sp_mat trueX;
|
|
trueX.sprandu(50, 50, 0.3);
|
|
|
|
sp_mat b = a * trueX;
|
|
|
|
// Get things into the right format.
|
|
mat db(b);
|
|
|
|
mat x;
|
|
|
|
spsolve(x, a, db);
|
|
|
|
for (uword i = 0; i < x.n_cols; ++i)
|
|
{
|
|
for (uword j = 0; j < x.n_rows; ++j)
|
|
{
|
|
REQUIRE( x(j, i) == Approx((double) trueX(j, i)).margin(0.01) );
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
TEST_CASE("fn_spsolve_float_superlu_solve_test", "[spsolve]")
|
|
{
|
|
// Solve this matrix, as in the examples:
|
|
// [[19 0 21 21 0]
|
|
// [12 21 0 0 0]
|
|
// [ 0 12 16 0 0]
|
|
// [ 0 0 0 5 21]
|
|
// [12 12 0 0 18]]
|
|
sp_fmat b(5, 5);
|
|
b(0, 0) = 19;
|
|
b(0, 2) = 21;
|
|
b(0, 3) = 21;
|
|
b(1, 0) = 12;
|
|
b(1, 1) = 21;
|
|
b(2, 1) = 12;
|
|
b(2, 2) = 16;
|
|
b(3, 3) = 5;
|
|
b(3, 4) = 21;
|
|
b(4, 0) = 12;
|
|
b(4, 1) = 12;
|
|
b(4, 4) = 18;
|
|
|
|
fmat db(b);
|
|
|
|
sp_fmat a;
|
|
a.eye(5, 5);
|
|
fmat da(a);
|
|
|
|
fmat x;
|
|
spsolve(x, a, db);
|
|
|
|
fmat dx = solve(da, db);
|
|
|
|
for (uword i = 0; i < x.n_cols; ++i)
|
|
{
|
|
for (uword j = 0; j < x.n_rows; ++j)
|
|
{
|
|
REQUIRE( (float) x(j, i) == Approx(dx(j, i)).margin(0.01) );
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
TEST_CASE("fn_spsolve_float_random_superlu_solve_test", "[spsolve]")
|
|
{
|
|
// Try to solve some random systems.
|
|
const size_t iterations = 10;
|
|
size_t successes = 0;
|
|
for (size_t it = 0; it < iterations; ++it)
|
|
{
|
|
sp_fmat a;
|
|
a.sprandu(50, 50, 0.3);
|
|
sp_fmat trueX;
|
|
trueX.sprandu(50, 50, 0.3);
|
|
|
|
sp_fmat b = a * trueX;
|
|
|
|
// Get things into the right format.
|
|
fmat db(b);
|
|
|
|
fmat x;
|
|
|
|
spsolve(x, a, db);
|
|
|
|
if (any(abs(vectorise(trueX - x)) > 0.04))
|
|
{
|
|
// Skip this trial.
|
|
continue;
|
|
}
|
|
|
|
for (uword i = 0; i < x.n_cols; ++i)
|
|
{
|
|
for (uword j = 0; j < x.n_rows; ++j)
|
|
{
|
|
if (std::abs(trueX(j, i)) < 0.002)
|
|
REQUIRE( std::abs(x(j, i)) < 0.04 );
|
|
else
|
|
REQUIRE( trueX(j, i) == Approx((float) x(j, i)).margin(0.04) );
|
|
}
|
|
}
|
|
|
|
++successes;
|
|
}
|
|
|
|
REQUIRE( successes > 8 );
|
|
}
|
|
|
|
|
|
|
|
TEST_CASE("fn_spsolve_cx_float_superlu_solve_test", "[spsolve]")
|
|
{
|
|
// Solve this matrix, as in the examples:
|
|
// [[19 0 21 21 0]
|
|
// [12 21 0 0 0]
|
|
// [ 0 12 16 0 0]
|
|
// [ 0 0 0 5 21]
|
|
// [12 12 0 0 18]] (imaginary part is the same)
|
|
SpMat<cx_float> b(5, 5);
|
|
b(0, 0) = cx_float(19, 19);
|
|
b(0, 2) = cx_float(21, 21);
|
|
b(0, 3) = cx_float(21, 21);
|
|
b(1, 0) = cx_float(12, 12);
|
|
b(1, 1) = cx_float(21, 21);
|
|
b(2, 1) = cx_float(12, 12);
|
|
b(2, 2) = cx_float(16, 16);
|
|
b(3, 3) = cx_float(5, 5);
|
|
b(3, 4) = cx_float(21, 21);
|
|
b(4, 0) = cx_float(12, 12);
|
|
b(4, 1) = cx_float(12, 12);
|
|
b(4, 4) = cx_float(18, 18);
|
|
|
|
Mat<cx_float> db(b);
|
|
|
|
SpMat<cx_float> a;
|
|
a.eye(5, 5);
|
|
Mat<cx_float> da(a);
|
|
|
|
Mat<cx_float> x;
|
|
spsolve(x, a, db);
|
|
|
|
Mat<cx_float> dx = solve(da, db);
|
|
|
|
for (uword i = 0; i < x.n_cols; ++i)
|
|
{
|
|
for (uword j = 0; j < x.n_rows; ++j)
|
|
{
|
|
if (std::abs(x(j, i)) < 0.002 )
|
|
{
|
|
REQUIRE( std::abs(dx(j, i)) < 0.02 );
|
|
}
|
|
else
|
|
{
|
|
REQUIRE( ((cx_float) x(j, i)).real() ==
|
|
Approx(dx(j, i).real()).margin(0.02) );
|
|
REQUIRE( ((cx_float) x(j, i)).imag() ==
|
|
Approx(dx(j, i).imag()).margin(0.02) );
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
TEST_CASE("fn_spsolve_cx_float_random_superlu_solve_test", "[spsolve]")
|
|
{
|
|
// Try to solve some random systems.
|
|
const size_t iterations = 10;
|
|
for (size_t it = 0; it < iterations; ++it)
|
|
{
|
|
SpMat<cx_float> a;
|
|
a.sprandu(50, 50, 0.3);
|
|
SpMat<cx_float> trueX;
|
|
trueX.sprandu(50, 50, 0.3);
|
|
|
|
SpMat<cx_float> b = a * trueX;
|
|
|
|
// Get things into the right format.
|
|
Mat<cx_float> db(b);
|
|
|
|
Mat<cx_float> x;
|
|
|
|
spsolve(x, a, db);
|
|
|
|
for (uword i = 0; i < x.n_cols; ++i)
|
|
{
|
|
for (uword j = 0; j < x.n_rows; ++j)
|
|
{
|
|
if (std::abs((cx_float) trueX(j, i)) < 0.002 )
|
|
{
|
|
REQUIRE( std::abs(x(j, i)) < 0.02 );
|
|
}
|
|
else
|
|
{
|
|
REQUIRE( ((cx_float) trueX(j, i)).real() ==
|
|
Approx(x(j, i).real()).margin(0.02) );
|
|
REQUIRE( ((cx_float) trueX(j, i)).imag() ==
|
|
Approx(x(j, i).imag()).margin(0.02) );
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
TEST_CASE("fn_spsolve_cx_superlu_solve_test", "[spsolve]")
|
|
{
|
|
// Solve this matrix, as in the examples:
|
|
// [[19 0 21 21 0]
|
|
// [12 21 0 0 0]
|
|
// [ 0 12 16 0 0]
|
|
// [ 0 0 0 5 21]
|
|
// [12 12 0 0 18]] (imaginary part is the same)
|
|
SpMat<cx_double> b(5, 5);
|
|
b(0, 0) = cx_double(19, 19);
|
|
b(0, 2) = cx_double(21, 21);
|
|
b(0, 3) = cx_double(21, 21);
|
|
b(1, 0) = cx_double(12, 12);
|
|
b(1, 1) = cx_double(21, 21);
|
|
b(2, 1) = cx_double(12, 12);
|
|
b(2, 2) = cx_double(16, 16);
|
|
b(3, 3) = cx_double(5, 5);
|
|
b(3, 4) = cx_double(21, 21);
|
|
b(4, 0) = cx_double(12, 12);
|
|
b(4, 1) = cx_double(12, 12);
|
|
b(4, 4) = cx_double(18, 18);
|
|
|
|
cx_mat db(b);
|
|
|
|
sp_cx_mat a;
|
|
a.eye(5, 5);
|
|
cx_mat da(a);
|
|
|
|
cx_mat x;
|
|
spsolve(x, a, db);
|
|
|
|
cx_mat dx = solve(da, db);
|
|
|
|
for (uword i = 0; i < x.n_cols; ++i)
|
|
{
|
|
for (uword j = 0; j < x.n_rows; ++j)
|
|
{
|
|
if (std::abs(x(j, i)) < 0.001)
|
|
{
|
|
REQUIRE( std::abs(dx(j, i)) < 0.005 );
|
|
}
|
|
else
|
|
{
|
|
REQUIRE( ((cx_double) x(j, i)).real() ==
|
|
Approx(dx(j, i).real()).margin(0.01) );
|
|
REQUIRE( ((cx_double) x(j, i)).imag() ==
|
|
Approx(dx(j, i).imag()).margin(0.01) );
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
TEST_CASE("fn_spsolve_cx_random_superlu_solve_test", "[spsolve]")
|
|
{
|
|
// Try to solve some random systems.
|
|
const size_t iterations = 10;
|
|
for (size_t it = 0; it < iterations; ++it)
|
|
{
|
|
sp_cx_mat a;
|
|
a.sprandu(50, 50, 0.3);
|
|
sp_cx_mat trueX;
|
|
trueX.sprandu(50, 50, 0.3);
|
|
|
|
sp_cx_mat b = a * trueX;
|
|
|
|
// Get things into the right format.
|
|
cx_mat db(b);
|
|
|
|
cx_mat x;
|
|
|
|
spsolve(x, a, db);
|
|
|
|
for (uword i = 0; i < x.n_cols; ++i)
|
|
{
|
|
for (uword j = 0; j < x.n_rows; ++j)
|
|
{
|
|
if (std::abs((cx_double) trueX(j, i)) < 0.001)
|
|
{
|
|
REQUIRE( std::abs(x(j, i)) < 0.005 );
|
|
}
|
|
else
|
|
{
|
|
REQUIRE( ((cx_double) trueX(j, i)).real() ==
|
|
Approx(x(j, i).real()).margin(0.01) );
|
|
REQUIRE( ((cx_double) trueX(j, i)).imag() ==
|
|
Approx(x(j, i).imag()).margin(0.01) );
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
TEST_CASE("fn_spsolve_function_test", "[spsolve]")
|
|
{
|
|
sp_mat a;
|
|
a.sprandu(50, 50, 0.3);
|
|
sp_mat trueX;
|
|
trueX.sprandu(50, 50, 0.3);
|
|
|
|
sp_mat b = a * trueX;
|
|
|
|
// Get things into the right format.
|
|
mat db(b);
|
|
|
|
mat x;
|
|
|
|
// Mostly these are compilation tests.
|
|
spsolve(x, a, db);
|
|
x = spsolve(a, db); // Test another overload.
|
|
x = spsolve(a, db + 0.0);
|
|
spsolve(x, a, db + 0.0);
|
|
|
|
for (uword i = 0; i < x.n_cols; ++i)
|
|
{
|
|
for (uword j = 0; j < x.n_rows; ++j)
|
|
{
|
|
REQUIRE( (double) trueX(j, i) == Approx(x(j, i)).margin(0.01) );
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
TEST_CASE("fn_spsolve_float_function_test", "[spsolve]")
|
|
{
|
|
sp_fmat a;
|
|
a.sprandu(50, 50, 0.3);
|
|
sp_fmat trueX;
|
|
trueX.sprandu(50, 50, 0.3);
|
|
|
|
sp_fmat b = a * trueX;
|
|
|
|
// Get things into the right format.
|
|
fmat db(b);
|
|
|
|
fmat x;
|
|
|
|
// Mostly these are compilation tests.
|
|
spsolve(x, a, db);
|
|
x = spsolve(a, db); // Test another overload.
|
|
x = spsolve(a, db + 0.0);
|
|
spsolve(x, a, db + 0.0);
|
|
|
|
for (uword i = 0; i < x.n_cols; ++i)
|
|
{
|
|
for (uword j = 0; j < x.n_rows; ++j)
|
|
{
|
|
if (std::abs(trueX(j, i)) < 0.001)
|
|
{
|
|
REQUIRE( std::abs(x(j, i)) < 0.02 );
|
|
}
|
|
else
|
|
{
|
|
REQUIRE( (float) trueX(j, i) == Approx(x(j, i)).margin(0.02) );
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
TEST_CASE("fn_spsolve_cx_function_test", "[spsolve]")
|
|
{
|
|
sp_cx_mat a;
|
|
a.sprandu(50, 50, 0.3);
|
|
sp_cx_mat trueX;
|
|
trueX.sprandu(50, 50, 0.3);
|
|
|
|
sp_cx_mat b = a * trueX;
|
|
|
|
// Get things into the right format.
|
|
cx_mat db(b);
|
|
|
|
cx_mat x;
|
|
|
|
// Mostly these are compilation tests.
|
|
spsolve(x, a, db);
|
|
x = spsolve(a, db); // Test another overload.
|
|
x = spsolve(a, db + cx_double(0.0));
|
|
spsolve(x, a, db + cx_double(0.0));
|
|
|
|
for (uword i = 0; i < x.n_cols; ++i)
|
|
{
|
|
for (uword j = 0; j < x.n_rows; ++j)
|
|
{
|
|
if (std::abs((cx_double) trueX(j, i)) < 0.001)
|
|
{
|
|
REQUIRE( std::abs(x(j, i)) < 0.005 );
|
|
}
|
|
else
|
|
{
|
|
REQUIRE( ((cx_double) trueX(j, i)).real() ==
|
|
Approx(x(j, i).real()).margin(0.01) );
|
|
REQUIRE( ((cx_double) trueX(j, i)).imag() ==
|
|
Approx(x(j, i).imag()).margin(0.01) );
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
TEST_CASE("fn_spsolve_cx_float_function_test", "[spsolve]")
|
|
{
|
|
sp_cx_fmat a;
|
|
a.sprandu(50, 50, 0.3);
|
|
sp_cx_fmat trueX;
|
|
trueX.sprandu(50, 50, 0.3);
|
|
|
|
sp_cx_fmat b = a * trueX;
|
|
|
|
// Get things into the right format.
|
|
cx_fmat db(b);
|
|
|
|
cx_fmat x;
|
|
|
|
// Mostly these are compilation tests.
|
|
spsolve(x, a, db);
|
|
x = spsolve(a, db); // Test another overload.
|
|
x = spsolve(a, db + cx_float(0.0));
|
|
spsolve(x, a, db + cx_float(0.0));
|
|
|
|
for (uword i = 0; i < x.n_cols; ++i)
|
|
{
|
|
for (uword j = 0; j < x.n_rows; ++j)
|
|
{
|
|
if (std::abs((cx_float) trueX(j, i)) < 0.002 )
|
|
{
|
|
REQUIRE( std::abs(x(j, i)) < 0.02 );
|
|
}
|
|
else
|
|
{
|
|
REQUIRE( ((cx_float) trueX(j, i)).real() ==
|
|
Approx(x(j, i).real()).margin(0.02) );
|
|
REQUIRE( ((cx_float) trueX(j, i)).imag() ==
|
|
Approx(x(j, i).imag()).margin(0.02) );
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
TEST_CASE("spsolve_factoriser_test", "[spsolve]")
|
|
{
|
|
sp_mat A;
|
|
A.sprandu(100, 100, 0.2);
|
|
A.diag().randu();
|
|
A.diag() += 1;
|
|
|
|
vec B(100, fill::randu);
|
|
|
|
vec X1;
|
|
bool X1_status = spsolve(X1, A, B);
|
|
REQUIRE( X1_status );
|
|
|
|
|
|
spsolve_factoriser SF;
|
|
bool SF_status = SF.factorise(A);
|
|
REQUIRE( SF_status );
|
|
|
|
double rcond_value = SF.rcond();
|
|
REQUIRE( rcond_value > 0.0 );
|
|
|
|
vec X2;
|
|
bool X2_status = SF.solve(X2, B);
|
|
REQUIRE( X2_status );
|
|
|
|
vec X3;
|
|
bool X3_status = SF.solve(X3, B);
|
|
REQUIRE( X3_status );
|
|
|
|
REQUIRE( X1.n_rows == X2.n_rows );
|
|
REQUIRE( X1.n_rows == X3.n_rows );
|
|
|
|
REQUIRE( approx_equal(X1, X2, "absdiff", 100.0*datum::eps) );
|
|
REQUIRE( approx_equal(X1, X3, "absdiff", 100.0*datum::eps) );
|
|
|
|
vec C(100, fill::randu);
|
|
vec Y;
|
|
|
|
bool Y_status = SF.solve(Y, C);
|
|
REQUIRE( Y_status );
|
|
|
|
REQUIRE( approx_equal(X1, Y, "absdiff", 100.0*datum::eps) == false );
|
|
|
|
vec D(101, fill::randu);
|
|
vec Z(100, fill::randu);
|
|
|
|
bool Z_status = SF.solve(Z, D);
|
|
REQUIRE( Z_status == false );
|
|
|
|
REQUIRE( Z.n_elem == 0 );
|
|
}
|
|
|
|
#endif
|