ProteoWizard
qrTest.cpp
Go to the documentation of this file.
1//
2// $Id$
3//
4// Licensed under the Apache License, Version 2.0 (the "License");
5// you may not use this file except in compliance with the License.
6// You may obtain a copy of the License at
7//
8// http://www.apache.org/licenses/LICENSE-2.0
9//
10// Unless required by applicable law or agreed to in writing, software
11// distributed under the License is distributed on an "AS IS" BASIS,
12// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13// See the License for the specific language governing permissions and
14// limitations under the License.
15//
16
17#include <boost/numeric/ublas/matrix.hpp>
18#include <boost/numeric/ublas/vector.hpp>
19#include <boost/numeric/ublas/io.hpp>
20
22//#include <execinfo.h>
23
25#include "Types.hpp"
26#include "qr.hpp"
27
28using namespace boost::numeric::ublas;
29using namespace pwiz::math;
30using namespace pwiz::math::types;
31using namespace pwiz::util;
32
33
34ostream* os_ = 0;
35const double epsilon = 1e-12;
36
37
38template<class matrix_type>
39bool isUpperTriangular(const matrix_type& A, double eps)
40{
41 typedef typename matrix_type::size_type size_type;
42 bool upper = true;
43
44 for (size_type i=1; i<A.size2() && upper; i++)
45 {
46 for (size_type j=0; j<i && upper; j++)
47 {
48 upper = fabs(A(i,j)) < eps;
49 }
50 }
51
52 return upper;
53}
54
55
57{
58 if (os_) *os_ << "testReflector() begin" << endl;
59
60 dmatrix F(3,3);
61 dvector x(3);
62
63 x(0) = 1;
64 x(1) = 1;
65 x(2) = 0;
66
67 Reflector(x, F);
68
69 dvector v = prod(F, x);
70
71 unit_assert_equal(fabs(v(0)), norm_2(x), epsilon);
72 unit_assert_equal(v(1), 0, epsilon);
73 unit_assert_equal(v(2), 0, epsilon);
74
75 x(0) = -1;
76 x(1) = 1;
77 x(2) = 0;
78
79 if (os_) *os_ << "testReflector() end" << endl;
80}
81
82void testQR()
83{
84 if (os_) *os_ << "testQR() begin" << endl;
85
86 dmatrix A(3,3);
87 dmatrix Q(3,3);
88 dmatrix R(3,3);
89
90 for (dmatrix::size_type i=0; i< A.size1(); i++)
91 {
92 for (dmatrix::size_type j=0; j<A.size2(); j++)
93 {
94 A(i,j) = ((i==j) || (j == 0));
95 }
96 }
97
98 try {
99 qr<dmatrix>(A, Q, R);
100
102
103 dmatrix diff = (A - prod(Q, R));
104
105 unit_assert_equal(norm_1(diff), 0, epsilon);
106
107 identity_matrix<dmatrix::value_type> eye(Q.size1());
108 diff = prod(Q, herm(Q)) - eye;
109
110 unit_assert_equal(norm_1(diff), 0, epsilon);
111 }
112 catch (boost::numeric::ublas::bad_argument ba)
113 {
114 if (os_) *os_ << "exception: " << ba.what() << endl;
115 }
116
117 if (os_) *os_ << "testQR() end" << endl;
118}
119
121{
122 if (os_) *os_ << "testRectangularQR() begin" << endl;
123
124 dmatrix A(5,3);
125 dmatrix Q;
126 dmatrix R;
127
128 for (dmatrix::size_type i=0; i< A.size1(); i++)
129 {
130 for (dmatrix::size_type j=0; j<A.size2(); j++)
131 {
132 A(i,j) = ((i==j) || (j == 0));
133 }
134 }
135
136 try {
137 qr<dmatrix>(A, Q, R);
138
140
141 dmatrix diff = (A - prod(Q, R));
142
143 unit_assert_equal(norm_1(diff), 0, epsilon);
144
145 identity_matrix<dmatrix::value_type> eye(Q.size1());
146 diff = prod(Q, herm(Q)) - eye;
147
148 unit_assert_equal(norm_1(diff), 0, epsilon);
149 }
150 catch (boost::numeric::ublas::bad_argument ba)
151 {
152 if (os_) *os_ << "exception: " << ba.what() << endl;
153 }
154
155 if (os_) *os_ << "testRectangularQR() end" << endl;
156}
157
158int main(int argc, char** argv)
159{
160 if (argc>1 && !strcmp(argv[1],"-v")) os_ = &cout;
161 if (os_) *os_ << "qrTest\n";
163 testQR();
165
166 return 0;
167}
F
Definition Chemistry.hpp:81
#define A
void diff(const string &filename1, const string &filename2)
KernelTraitsBase< Kernel >::space_type::abscissa_type x
boost::numeric::ublas::vector< double > dvector
Definition Types.hpp:34
boost::numeric::ublas::matrix< double > dmatrix
Definition Types.hpp:33
void Reflector(const vector_type &x, matrix_type &F)
Definition qr.hpp:39
void testQR()
Definition qrTest.cpp:82
int main(int argc, char **argv)
Definition qrTest.cpp:158
bool isUpperTriangular(const matrix_type &A, double eps)
Definition qrTest.cpp:39
void testReflector()
Definition qrTest.cpp:56
const double epsilon
Definition qrTest.cpp:35
void testRectangularQR()
Definition qrTest.cpp:120
ostream * os_
Definition qrTest.cpp:34
#define unit_assert(x)
Definition unit.hpp:85
#define unit_assert_equal(x, y, epsilon)
Definition unit.hpp:99