123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900 |
- // 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")
- {
- // 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);
- for (uword i = 0; i < size; ++i)
- {
- A(i, i) += rand();
- }
- 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)) );
- }
- }
- }
- }
- TEST_CASE("fn_spsolve_sparse_nonsymmetric_test")
- {
- 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);
- for (uword i = 0; i < r_size; ++i)
- {
- A(i, i) += rand();
- }
- 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)) );
- }
- }
- }
- }
- TEST_CASE("fn_spsolve_sparse_float_test")
- {
- // 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);
- fmat rX;
- rX.randu(size, size);
- SpMat<float> A;
- A.sprandu(size, size, 0.25);
- for (uword i = 0; i < size; ++i)
- {
- A(i, i) += rand();
- }
- 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 (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)) );
- }
- }
- }
- }
- TEST_CASE("fn_spsolve_sparse_nonsymmetric_float_test")
- {
- 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);
- for (uword i = 0; i < r_size; ++i)
- {
- A(i, i) += rand();
- }
- 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)) );
- }
- }
- }
- }
- TEST_CASE("fn_spsolve_sparse_complex_float_test")
- {
- // 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<std::complex<float> > rX;
- rX.randu(size, size);
- SpMat<std::complex<float> > A;
- A.sprandu(size, size, 0.25);
- for(uword i = 0; i < size; ++i)
- {
- A(i, i) += rand();
- }
- Mat<std::complex<float> > B = A * rX;
- Mat<std::complex<float> > X;
- bool result = spsolve(X, A, B);
- REQUIRE( result );
- // Dense solver.
- Mat<std::complex<float> > dA(A);
- Mat<std::complex<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((std::complex<float>) X(j, i)) ==
- Approx((float) std::abs((std::complex<float>) dX(j, i))) );
- }
- }
- }
- }
- TEST_CASE("fn_spsolve_sparse_nonsymmetric_complex_float_test")
- {
- for (size_t t = 0; t < 10; ++t)
- {
- const uword r_size = 5 * (t + 1);
- const uword c_size = 3 * (t + 4);
- Mat<std::complex<float> > rX;
- rX.randu(r_size, c_size);
- SpMat<std::complex<float> > A;
- A.sprandu(r_size, r_size, 0.25);
- for (uword i = 0; i < r_size; ++i)
- {
- A(i, i) += rand();
- }
- Mat<std::complex<float> > B = A * rX;
- Mat<std::complex<float> > X;
- bool result = spsolve(X, A, B);
- REQUIRE( result );
- // Dense solver.
- Mat<std::complex<float> > dA(A);
- Mat<std::complex<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((std::complex<float>) X(j, i)) ==
- Approx((float) std::abs((std::complex<float>) dX(j, i))) );
- }
- }
- }
- }
- TEST_CASE("fn_spsolve_sparse_complex_test")
- {
- // 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<std::complex<double> > rX;
- rX.randu(size, size);
- SpMat<std::complex<double> > A;
- A.sprandu(size, size, 0.25);
- for (uword i = 0; i < size; ++i)
- {
- A(i, i) += rand();
- }
- Mat<std::complex<double> > B = A * rX;
- Mat<std::complex<double> > X;
- bool result = spsolve(X, A, B);
- REQUIRE( result );
- // Dense solver.
- Mat<std::complex<double> > dA(A);
- Mat<std::complex<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((std::complex<double>) X(j, i)) ==
- Approx((double) std::abs((std::complex<double>) dX(j, i))) );
- }
- }
- }
- }
- TEST_CASE("fn_spsolve_sparse_nonsymmetric_complex_test")
- {
- for (size_t t = 0; t < 10; ++t)
- {
- const uword r_size = 5 * (t + 1);
- const uword c_size = 3 * (t + 4);
- Mat<std::complex<double> > rX;
- rX.randu(r_size, c_size);
- SpMat<std::complex<double> > A;
- A.sprandu(r_size, r_size, 0.25);
- for (uword i = 0; i < r_size; ++i)
- {
- A(i, i) += rand();
- }
- Mat<std::complex<double> > B = A * rX;
- Mat<std::complex<double> > X;
- bool result = spsolve(X, A, B);
- REQUIRE( result );
- // Dense solver.
- Mat<std::complex<double> > dA(A);
- Mat<std::complex<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((std::complex<double>) X(j, i)) ==
- Approx((double) std::abs((std::complex<double>) dX(j, i))) );
- }
- }
- }
- }
- TEST_CASE("fn_spsolve_delayed_sparse_test")
- {
- const uword size = 10;
- mat rX;
- rX.randu(size, size);
- sp_mat A;
- A.sprandu(size, size, 0.25);
- for (uword i = 0; i < size; ++i)
- {
- A(i, i) += rand();
- }
- 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)) );
- }
- }
- }
- TEST_CASE("fn_spsolve_superlu_solve_test")
- {
- // 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)) );
- }
- }
- }
- TEST_CASE("fn_spsolve_random_superlu_solve_test")
- {
- // 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)) );
- }
- }
- }
- }
- TEST_CASE("fn_spsolve_float_superlu_solve_test")
- {
- // 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)) );
- }
- }
- }
- TEST_CASE("fn_spsolve_float_random_superlu_solve_test")
- {
- // Try to solve some random systems.
- const size_t iterations = 10;
- 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);
- 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.005 );
- else
- REQUIRE( trueX(j, i) == Approx((float) x(j, i)).epsilon(0.01) );
- }
- }
- }
- }
- TEST_CASE("fn_spsolve_cx_float_superlu_solve_test")
- {
- // 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<std::complex<float> > b(5, 5);
- b(0, 0) = std::complex<float>(19, 19);
- b(0, 2) = std::complex<float>(21, 21);
- b(0, 3) = std::complex<float>(21, 21);
- b(1, 0) = std::complex<float>(12, 12);
- b(1, 1) = std::complex<float>(21, 21);
- b(2, 1) = std::complex<float>(12, 12);
- b(2, 2) = std::complex<float>(16, 16);
- b(3, 3) = std::complex<float>(5, 5);
- b(3, 4) = std::complex<float>(21, 21);
- b(4, 0) = std::complex<float>(12, 12);
- b(4, 1) = std::complex<float>(12, 12);
- b(4, 4) = std::complex<float>(18, 18);
- Mat<std::complex<float> > db(b);
- SpMat<std::complex<float> > a;
- a.eye(5, 5);
- Mat<std::complex<float> > da(a);
- Mat<std::complex<float> > x;
- spsolve(x, a, db);
- Mat<std::complex<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.001 )
- {
- REQUIRE( std::abs(dx(j, i)) < 0.005 );
- }
- else
- {
- REQUIRE( ((std::complex<float>) x(j, i)).real() ==
- Approx(dx(j, i).real()).epsilon(0.01) );
- REQUIRE( ((std::complex<float>) x(j, i)).imag() ==
- Approx(dx(j, i).imag()).epsilon(0.01) );
- }
- }
- }
- }
- TEST_CASE("fn_spsolve_cx_float_random_superlu_solve_test")
- {
- // Try to solve some random systems.
- const size_t iterations = 10;
- for (size_t it = 0; it < iterations; ++it)
- {
- SpMat<std::complex<float> > a;
- a.sprandu(50, 50, 0.3);
- SpMat<std::complex<float> > trueX;
- trueX.sprandu(50, 50, 0.3);
- SpMat<std::complex<float> > b = a * trueX;
- // Get things into the right format.
- Mat<std::complex<float> > db(b);
- Mat<std::complex<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((std::complex<float>) trueX(j, i)) < 0.001 )
- {
- REQUIRE( std::abs(x(j, i)) < 0.001 );
- }
- else
- {
- REQUIRE( ((std::complex<float>) trueX(j, i)).real() ==
- Approx(x(j, i).real()).epsilon(0.01) );
- REQUIRE( ((std::complex<float>) trueX(j, i)).imag() ==
- Approx(x(j, i).imag()).epsilon(0.01) );
- }
- }
- }
- }
- }
- TEST_CASE("fn_spsolve_cx_superlu_solve_test")
- {
- // 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<std::complex<double> > b(5, 5);
- b(0, 0) = std::complex<double>(19, 19);
- b(0, 2) = std::complex<double>(21, 21);
- b(0, 3) = std::complex<double>(21, 21);
- b(1, 0) = std::complex<double>(12, 12);
- b(1, 1) = std::complex<double>(21, 21);
- b(2, 1) = std::complex<double>(12, 12);
- b(2, 2) = std::complex<double>(16, 16);
- b(3, 3) = std::complex<double>(5, 5);
- b(3, 4) = std::complex<double>(21, 21);
- b(4, 0) = std::complex<double>(12, 12);
- b(4, 1) = std::complex<double>(12, 12);
- b(4, 4) = std::complex<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( ((std::complex<double>) x(j, i)).real() ==
- Approx(dx(j, i).real()).epsilon(0.01) );
- REQUIRE( ((std::complex<double>) x(j, i)).imag() ==
- Approx(dx(j, i).imag()).epsilon(0.01) );
- }
- }
- }
- }
- TEST_CASE("fn_spsolve_cx_random_superlu_solve_test")
- {
- // 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((std::complex<double>) trueX(j, i)) < 0.001)
- {
- REQUIRE( std::abs(x(j, i)) < 0.005 );
- }
- else
- {
- REQUIRE( ((std::complex<double>) trueX(j, i)).real() ==
- Approx(x(j, i).real()).epsilon(0.01) );
- REQUIRE( ((std::complex<double>) trueX(j, i)).imag() ==
- Approx(x(j, i).imag()).epsilon(0.01) );
- }
- }
- }
- }
- }
- TEST_CASE("fn_spsolve_function_test")
- {
- 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)) );
- }
- }
- }
- TEST_CASE("fn_spsolve_float_function_test")
- {
- 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.001 );
- }
- else
- {
- REQUIRE( (float) trueX(j, i) == Approx(x(j, i)).epsilon(0.01) );
- }
- }
- }
- }
- TEST_CASE("fn_spsolve_cx_function_test")
- {
- 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 + std::complex<double>(0.0));
- spsolve(x, a, db + std::complex<double>(0.0));
- for (uword i = 0; i < x.n_cols; ++i)
- {
- for (uword j = 0; j < x.n_rows; ++j)
- {
- if (std::abs((std::complex<double>) trueX(j, i)) < 0.001)
- {
- REQUIRE( std::abs(x(j, i)) < 0.005 );
- }
- else
- {
- REQUIRE( ((std::complex<double>) trueX(j, i)).real() ==
- Approx(x(j, i).real()).epsilon(0.01) );
- REQUIRE( ((std::complex<double>) trueX(j, i)).imag() ==
- Approx(x(j, i).imag()).epsilon(0.01) );
- }
- }
- }
- }
- TEST_CASE("fn_spsolve_cx_float_function_test")
- {
- 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 + std::complex<float>(0.0));
- spsolve(x, a, db + std::complex<float>(0.0));
- for (uword i = 0; i < x.n_cols; ++i)
- {
- for (uword j = 0; j < x.n_rows; ++j)
- {
- if (std::abs((std::complex<float>) trueX(j, i)) < 0.001 )
- {
- REQUIRE( std::abs(x(j, i)) < 0.005 );
- }
- else
- {
- REQUIRE( ((std::complex<float>) trueX(j, i)).real() ==
- Approx(x(j, i).real()).epsilon(0.01) );
- REQUIRE( ((std::complex<float>) trueX(j, i)).imag() ==
- Approx(x(j, i).imag()).epsilon(0.01) );
- }
- }
- }
- }
- #endif
|