Line data Source code
1 : /*
2 : * File: TestCholesky.cpp
3 : * Author: Pantelis Sopasakis
4 : *
5 : * Created on Aug 5, 2015, 12:28:55 AM
6 : *
7 : * ForBES is free software: you can redistribute it and/or modify
8 : * it under the terms of the GNU Lesser General Public License as published by
9 : * the Free Software Foundation, either version 3 of the License, or
10 : * (at your option) any later version.
11 : *
12 : * ForBES is distributed in the hope that it will be useful,
13 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 : * GNU Lesser General Public License for more details.
16 : *
17 : * You should have received a copy of the GNU Lesser General Public License
18 : * along with ForBES. If not, see <http://www.gnu.org/licenses/>.
19 : */
20 :
21 : #include "TestCholesky.h"
22 : #include <cmath>
23 :
24 :
25 1 : CPPUNIT_TEST_SUITE_REGISTRATION(TestCholesky);
26 :
27 4 : TestCholesky::TestCholesky() {
28 4 : }
29 :
30 8 : TestCholesky::~TestCholesky() {
31 8 : }
32 :
33 4 : void TestCholesky::setUp() {
34 4 : }
35 :
36 4 : void TestCholesky::tearDown() {
37 4 : }
38 :
39 1 : void TestCholesky::testCholeskyDense() {
40 1 : size_t n = 30;
41 1 : size_t repetitions = 80;
42 1 : const double tol = 1e-7;
43 1 : Matrix A = MatrixFactory::MakeRandomMatrix(n, n, 0.0, 1.0, Matrix::MATRIX_DENSE);
44 2 : Matrix At = A;
45 1 : At.transpose();
46 1 : A += At;
47 1 : A *= 0.5;
48 31 : for (size_t i = 0; i < n; i++) {
49 30 : A.set(i, i, A.get(i, i) + 1.2 * n);
50 : }
51 2 : Matrix A_copy(A);
52 2 : Matrix b;
53 :
54 1 : CholeskyFactorization * cholFactorization = new CholeskyFactorization(A);
55 1 : _ASSERT_EQ(ForBESUtils::STATUS_OK, cholFactorization ->factorize());
56 :
57 : // Goodbye A...
58 1 : A = Matrix();
59 :
60 2 : Matrix x;
61 2 : Matrix err;
62 81 : for (size_t k = 0; k < repetitions; k++) {
63 80 : b = MatrixFactory::MakeRandomMatrix(n, 1, 0.0, 1.0, Matrix::MATRIX_DENSE);
64 80 : cholFactorization -> solve(b, x);
65 80 : err = A_copy * x - b;
66 2480 : for (size_t i = 0; i < n; i++) {
67 2400 : _ASSERT(std::abs(err.get(i, 0)) < tol);
68 : }
69 : }
70 2 : delete cholFactorization;
71 1 : }
72 :
73 1 : void TestCholesky::testCholeskySparse() {
74 1 : const double tol = 1e-7;
75 1 : size_t n = 25;
76 1 : size_t nnz = 2 * n - 1;
77 1 : Matrix A = MatrixFactory::MakeSparseSymmetric(n, nnz);
78 1 : _ASSERT_EQ(Matrix::MATRIX_SPARSE, A.getType());
79 2 : Matrix b(n, 1);
80 :
81 26 : for (size_t i = 0; i < n; i++) {
82 25 : A.set(i, i, n + 2.5);
83 25 : b.set(i, 0, i + 1);
84 : }
85 25 : for (size_t i = 1; i < n; i++) { /* Set the LT part only */
86 24 : A.set(i, i - 1, 0.5);
87 : }
88 :
89 :
90 2 : Matrix A_copy(A);
91 2 : Matrix x;
92 1 : FactoredSolver * solver = new CholeskyFactorization(A);
93 1 : _ASSERT_EQ(ForBESUtils::STATUS_OK, solver -> factorize());
94 1 : A = Matrix(); // Goodbye A...
95 1 : _ASSERT_EQ(ForBESUtils::STATUS_OK, solver -> solve(b, x));
96 :
97 2 : Matrix ax = A_copy * x;
98 26 : for (size_t i = 0; i < n; i++) {
99 25 : _ASSERT(std::abs(ax.get(i, 0) - i - 1.0) < tol);
100 : }
101 2 : delete solver;
102 1 : }
103 :
104 1 : void TestCholesky::testCholeskySymmetric2() {
105 1 : const size_t n = 4;
106 1 : Matrix *A = new Matrix(n, n, Matrix::MATRIX_SYMMETRIC);
107 1 : _ASSERT(A->isSymmetric());
108 1 : _ASSERT_EQ(Matrix::MATRIX_SYMMETRIC, A -> getType());
109 5 : for (size_t i = 0; i < n; i++) {
110 14 : for (size_t j = 0; j <= i; j++) {
111 10 : A -> set(i, j, 3.2 * i + 0.2 * j + 0.45);
112 10 : if (i == j) A -> set(i, i, A->get(i, i) + 20.0);
113 : }
114 : }
115 :
116 1 : Matrix L;
117 2 : Matrix b(n, 1);
118 5 : for (size_t i = 0; i < n; i++) {
119 4 : b.set(i, 0, i + 1.0);
120 : }
121 :
122 1 : FactoredSolver * solver = new CholeskyFactorization(*A);
123 1 : _ASSERT_EQ(ForBESUtils::STATUS_OK, solver -> factorize());
124 :
125 : // Bye bye Matrix A
126 1 : delete A;
127 :
128 : double sol_expected_data[4] = {-0.032619563872026,
129 : 0.020584586456240,
130 : 0.070701048816070,
131 1 : 0.110212353853490};
132 2 : Matrix sol_expected(4, 1, sol_expected_data);
133 :
134 2 : Matrix sol;
135 1 : _ASSERT_EQ(ForBESUtils::STATUS_OK, solver -> solve(b, sol));
136 :
137 1 : _ASSERT_EQ(sol_expected, sol);
138 :
139 2 : delete solver;
140 1 : }
141 :
142 1 : void TestCholesky::testCholeskySymmetric() {
143 1 : const double tol = 1e-7;
144 1 : size_t n = 30;
145 1 : size_t repetitions = 10;
146 1 : Matrix A = MatrixFactory::MakeRandomMatrix(n, n, 0.0, 1.0, Matrix::MATRIX_SYMMETRIC);
147 31 : for (size_t i = 0; i < n; i++) {
148 30 : A.set(i, i, A.get(i, i) + 1.2 * n);
149 : }
150 2 : Matrix A_copy(A);
151 2 : Matrix b;
152 2 : Matrix x;
153 2 : Matrix err;
154 : FactoredSolver * cholFactorization;
155 1 : cholFactorization = new CholeskyFactorization(A);
156 1 : _ASSERT_EQ(ForBESUtils::STATUS_OK, cholFactorization ->factorize());
157 :
158 1 : A = Matrix();
159 11 : for (size_t k = 0; k < repetitions; k++) {
160 10 : b = MatrixFactory::MakeRandomMatrix(n, 1, 0.0, 1.0, Matrix::MATRIX_DENSE);
161 10 : cholFactorization -> solve(b, x);
162 10 : err = A_copy * x - b;
163 310 : for (size_t i = 0; i < n; i++) {
164 300 : _ASSERT(std::abs(err.get(i, 0)) < tol);
165 : }
166 : }
167 :
168 2 : delete cholFactorization;
169 4 : }
170 :
171 :
172 :
|
