ProteoWizard
LinearSolver.hpp
Go to the documentation of this file.
1//
2// $Id$
3//
4//
5// Original author: Darren Kessner <darren@proteowizard.org>
6//
7// Copyright 2007 Spielberg Family Center for Applied Proteomics
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#ifndef _LINEARSOLVER_HPP_
25#define _LINEARSOLVER_HPP_
26
27
29#include <boost/numeric/ublas/vector.hpp>
30#include <boost/numeric/ublas/matrix.hpp>
31#include <boost/numeric/ublas/lu.hpp>
32#include <boost/numeric/ublas/triangular.hpp>
33#include <boost/numeric/ublas/vector_proxy.hpp>
34#include <boost/numeric/ublas/io.hpp>
35#include <stdexcept>
36
37#include <iostream>
38#include <iomanip>
39#include <math.h>
40
41#include "qr.hpp"
42
43namespace pwiz {
44namespace math {
45
46
47enum PWIZ_API_DECL LinearSolverType {LinearSolverType_LU, LinearSolverType_QR};
48
49
50template <LinearSolverType solver_type = LinearSolverType_LU>
52
53
54template<>
56{
57public:
58
59 /// solve system of linear equations Ax = y using boost::ublas;
60 /// note: extra copying inefficiencies for ease of client use
61 template<typename matrix_type, typename vector_type>
62 vector_type solve(const matrix_type& A,
63 const vector_type& y)
64 {
65 namespace ublas = boost::numeric::ublas;
66
67 matrix_type A_factorized = A;
68 ublas::permutation_matrix<size_t> pm(y.size());
69
70 int singular = lu_factorize(A_factorized, pm);
71 if (singular) throw std::runtime_error("[LinearSolver<LU>::solve()] A is singular.");
72
73 vector_type result(y);
74 lu_substitute(A_factorized, pm, result);
75
76 return result;
77 }
78};
79
80
81template<>
82class LinearSolver<LinearSolverType_QR>
83{
84public:
85
86 /// solve system of linear equations Ax = y using boost::ublas;
87 /// note: extra copying inefficiencies for ease of client use
88 template<typename matrix_type, typename vector_type>
89 vector_type solve(const matrix_type& A, const vector_type& y)
90 {
91 typedef typename matrix_type::size_type size_type;
92 typedef typename matrix_type::value_type value_type;
93
94 namespace ublas = boost::numeric::ublas;
95
96 matrix_type Q(A.size1(), A.size2()), R(A.size1(), A.size2());
97
98 qr (A, Q, R);
99
100 vector_type b = prod(trans(Q), y);
101
102 vector_type result;
103 if (R.size1() > R.size2())
104 {
105 size_type min = (R.size1() < R.size2() ? R.size1() : R.size2());
106
107 result = ublas::solve(subrange(R, 0, min, 0, min),
108 subrange(b, 0, min),
109 ublas::upper_tag());
110 }
111 else
112 {
113 result = ublas::solve(R, b, ublas::upper_tag());
114 }
115 return result;
116 }
117};
118
119} // namespace math
120} // namespace pwiz
121
122
123#endif // _LINEARSOLVER_HPP_
124
#define PWIZ_API_DECL
Definition Export.hpp:32
#define A
LinearSolverType_LU
KernelTraitsBase< Kernel >::space_type::ordinate_type y
vector_type solve(const matrix_type &A, const vector_type &y)
solve system of linear equations Ax = y using boost::ublas; note: extra copying inefficiencies for ea...
vector_type solve(const matrix_type &A, const vector_type &y)
solve system of linear equations Ax = y using boost::ublas; note: extra copying inefficiencies for ea...
void qr(const matrix_type &A, matrix_type &Q, matrix_type &R)
Definition qr.hpp:81