ProteoWizard
LinearSolverTest.cpp
Go to the documentation of this file.
1//
2// $Id$
3//
4//
5// Original author: Darren Kessner <darren@proteowizard.org>
6//
7// Copyright 2006 Louis Warschaw Prostate Cancer Center
8// Cedars Sinai Medical Center, Los Angeles, California 90048
9//
10// Licensed under the Apache License, Version 2.0 (the "License");
11// you may not use this file except in compliance with the License.
12// You may obtain a copy of the License at
13//
14// http://www.apache.org/licenses/LICENSE-2.0
15//
16// Unless required by applicable law or agreed to in writing, software
17// distributed under the License is distributed on an "AS IS" BASIS,
18// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
19// See the License for the specific language governing permissions and
20// limitations under the License.
21//
22
23
24#include "LinearSolver.hpp"
26#include <boost/numeric/ublas/matrix_sparse.hpp>
27#include <boost/numeric/ublas/banded.hpp>
29#include <cstring>
30
31
32using namespace pwiz::util;
33using namespace pwiz::math;
34
35
36namespace ublas = boost::numeric::ublas;
37
38
39ostream* os_ = 0;
40
41
43{
44 if (os_) *os_ << "testDouble()\n";
45
46 LinearSolver<> solver;
47
48 ublas::matrix<double> A(2,2);
49 A(0,0) = 1; A(0,1) = 2;
50 A(1,0) = 3; A(1,1) = 4;
51
52 ublas::vector<double> y(2);
53 y(0) = 5;
54 y(1) = 11;
55
56 ublas::vector<double> x = solver.solve(A, y);
57
58 if (os_) *os_ << "A: " << A << endl;
59 if (os_) *os_ << "y: " << y << endl;
60 if (os_) *os_ << "x: " << x << endl;
61
62 unit_assert(x(0) == 1.);
63 unit_assert(x(1) == 2.);
64}
65
66
68{
69 if (os_) *os_ << "testComplex()\n";
70
71 LinearSolver<> solver;
72
73 ublas::matrix< complex<double> > A(2,2);
74 A(0,0) = 1; A(0,1) = 2;
75 A(1,0) = 3; A(1,1) = 4;
76
77 ublas::vector< complex<double> > y(2);
78 y(0) = 5;
79 y(1) = 11;
80
81 ublas::vector< complex<double> > x = solver.solve(A, y);
82
83 if (os_) *os_ << "A: " << A << endl;
84 if (os_) *os_ << "y: " << y << endl;
85 if (os_) *os_ << "x: " << x << endl;
86
87 unit_assert(x(0) == 1.);
88 unit_assert(x(1) == 2.);
89}
90
92{
93 if (os_) *os_ << "testDoubleQR()\n";
94
96
97 ublas::matrix<double> A(2,2);
98 A(0,0) = 1.; A(0,1) = 2.;
99 A(1,0) = 3.; A(1,1) = 4.;
100
101 ublas::vector<double> y(2);
102 y(0) = 5.;
103 y(1) = 11.;
104
105 ublas::vector<double> x = solver.solve(A, y);
106
107 if (os_) *os_ << "A: " << A << endl;
108 if (os_) *os_ << "y: " << y << endl;
109 if (os_) *os_ << "x: " << x << endl;
110
111 if (os_) *os_ << x(0) << " - 1. = " << x(0) - 1. << endl;
112
113 unit_assert_equal(x(0), 1., 1e-14);
114 unit_assert_equal(x(1), 2., 1e-14);
115}
116
117/*
118void testComplexQR()
119{
120 if (os_) *os_ << "testComplex()\n";
121
122 LinearSolver<LinearSolverType_QR> solver;
123
124 ublas::matrix< complex<double> > A(2,2);
125 A(0,0) = 1; A(0,1) = 2;
126 A(1,0) = 3; A(1,1) = 4;
127
128 ublas::vector< complex<double> > y(2);
129 y(0) = 5;
130 y(1) = 11;
131
132 ublas::vector< complex<double> > x = solver.solve(A, y);
133
134 if (os_) *os_ << "A: " << A << endl;
135 if (os_) *os_ << "y: " << y << endl;
136 if (os_) *os_ << "x: " << x << endl;
137
138 unit_assert(x(0) == 1.);
139 unit_assert(x(1) == 2.);
140}
141*/
142
143
145{
146 if (os_) *os_ << "testSparse()\n";
147
148 LinearSolver<> solver;
149
150 ublas::mapped_matrix<double> A(2,2,4);
151 A(0,0) = 1.; A(0,1) = 2.;
152 A(1,0) = 3.; A(1,1) = 4.;
153
154 ublas::vector<double> y(2);
155 y(0) = 5.;
156 y(1) = 11.;
157
158 ublas::vector<double> x = solver.solve(A, y);
159
160 if (os_) *os_ << "A: " << A << endl;
161 if (os_) *os_ << "y: " << y << endl;
162 if (os_) *os_ << "x: " << x << endl;
163
164 unit_assert_equal(x(0), 1., 1e-14);
165 unit_assert_equal(x(1), 2., 1e-14);
166}
167
168
169/*
170void testSparseComplex()
171{
172 if (os_) *os_ << "testSparseComplex()\n";
173
174 LinearSolver<> solver;
175
176 ublas::mapped_matrix< complex<double> > A(2,2,4);
177 A(0,0) = 1.; A(0,1) = 2.;
178 A(1,0) = 3.; A(1,1) = 4.;
179
180 ublas::vector< complex<double> > y(2);
181 y(0) = 5.;
182 y(1) = 11.;
183
184 ublas::vector< complex<double> > x = solver.solve(A, y);
185
186 if (os_) *os_ << "A: " << A << endl;
187 if (os_) *os_ << "y: " << y << endl;
188 if (os_) *os_ << "x: " << x << endl;
189
190 unit_assert(norm(x(0)-1.) < 1e-14);
191 unit_assert(norm(x(1)-2.) < 1e-14);
192}
193*/
194
195
197{
198 if (os_) *os_ << "testBanded()\n";
199
200 LinearSolver<> solver;
201
202 ublas::banded_matrix<double> A(2,2,1,1);
203 A(0,0) = 1.; A(0,1) = 2.;
204 A(1,0) = 3.; A(1,1) = 4.;
205
206 ublas::vector<double> y(2);
207 y(0) = 5.;
208 y(1) = 11.;
209
210 ublas::vector<double> x = solver.solve(A, y);
211
212 if (os_) *os_ << "A: " << A << endl;
213 if (os_) *os_ << "y: " << y << endl;
214 if (os_) *os_ << "x: " << x << endl;
215
216 unit_assert_equal(x(0), 1., 1e-14);
217 unit_assert_equal(x(1), 2., 1e-14);
218}
219
220
222{
223 if (os_) *os_ << "testBandedComplex()\n";
224
225 LinearSolver<> solver;
226
227 ublas::banded_matrix< complex<double> > A(2,2,1,1);
228 A(0,0) = 1.; A(0,1) = 2.;
229 A(1,0) = 3.; A(1,1) = 4.;
230
231 ublas::vector< complex<double> > y(2);
232 y(0) = 5.;
233 y(1) = 11.;
234
235 ublas::vector< complex<double> > x = solver.solve(A, y);
236
237 if (os_) *os_ << "A: " << A << endl;
238 if (os_) *os_ << "y: " << y << endl;
239 if (os_) *os_ << "x: " << x << endl;
240
241 unit_assert(norm(x(0)-1.) < 1e-14);
242 unit_assert(norm(x(1)-2.) < 1e-14);
243}
244
245
246int main(int argc, char* argv[])
247{
248 TEST_PROLOG(argc, argv)
249
250 try
251 {
252 if (argc>1 && !strcmp(argv[1],"-v")) os_ = &cout;
253 if (os_) *os_ << "LinearSolverTest\n";
254
255 testDouble();
256 testComplex();
257 testDoubleQR();
258 //testComplexQR();
259 testSparse();
260 //testSparseComplex(); // lu_factorize doesn't like mapped_matrix<complex>
261 testBanded();
262 //testBandedComplex(); // FIXME: GCC 4.2 doesn't like this test with link=shared
263 }
264 catch (exception& e)
265 {
266 TEST_FAILED(e.what())
267 }
268 catch (...)
269 {
270 TEST_FAILED("Caught unknown exception.")
271 }
272
274}
275
#define A
int main(int argc, char *argv[])
void testBandedComplex()
void testSparse()
void testBanded()
void testComplex()
void testDoubleQR()
void testDouble()
ostream * os_
KernelTraitsBase< Kernel >::space_type::abscissa_type x
KernelTraitsBase< Kernel >::space_type::ordinate_type y
#define unit_assert(x)
Definition unit.hpp:85
#define TEST_EPILOG
Definition unit.hpp:183
#define TEST_FAILED(x)
Definition unit.hpp:177
#define unit_assert_equal(x, y, epsilon)
Definition unit.hpp:99
#define TEST_PROLOG(argc, argv)
Definition unit.hpp:175