MVE - Multi-View Environment mve-devel
Loading...
Searching...
No Matches
ba_conjugate_gradient.h
Go to the documentation of this file.
1/*
2 * Copyright (C) 2015, Fabian Langguth, Simon Fuhrmann
3 * TU Darmstadt - Graphics, Capture and Massively Parallel Computing
4 * All rights reserved.
5 *
6 * This software may be modified and distributed under the terms
7 * of the BSD 3-Clause license. See the LICENSE.txt file for details.
8 */
9
10#ifndef SFM_BA_CONJUGATE_GRADIENT_HEADER
11#define SFM_BA_CONJUGATE_GRADIENT_HEADER
12
13#include "sfm/defines.h"
14#include "sfm/ba_dense_vector.h"
16
19
20template <typename T>
22{
23public:
26
28 {
31 CG_INVALID_INPUT
32 };
33
34 struct Options
35 {
36 Options (void);
39 };
40
41 struct Status
42 {
43 Status (void);
46 };
47
48 class Functor
49 {
50 public:
51 virtual Vector multiply (Vector const& x) const = 0;
52 virtual std::size_t input_size (void) const = 0;
53 virtual std::size_t output_size (void) const = 0;
54 };
55
56public:
57 ConjugateGradient (Options const& opts);
58
59 Status solve (Matrix const& A, Vector const& b, Vector* x,
60 Matrix const* P = nullptr);
61
62 Status solve (Functor const& A, Vector const& b, Vector* x,
63 Functor const* P = nullptr);
64
65private:
66 Options opts;
67 Status status;
68};
69
70template <typename T>
71class CGBasicMatrixFunctor : public ConjugateGradient<T>::Functor
72{
73public:
75 DenseVector<T> multiply (DenseVector<T> const& x) const;
76 std::size_t input_size (void) const;
77 std::size_t output_size (void) const;
78
79private:
80 SparseMatrix<T> const* A;
81};
82
83/* ------------------------ Implementation ------------------------ */
84
85template <typename T>
86inline
88 : max_iterations(100)
89 , tolerance(1e-20)
90{
91}
92
93template <typename T>
94inline
96 : num_iterations(0)
97{
98}
99
100template <typename T>
101inline
103 (Options const& options)
104 : opts(options)
105{
106}
107
108template <typename T>
109inline typename ConjugateGradient<T>::Status
111 Matrix const* P)
112{
113 CGBasicMatrixFunctor<T> A_functor(A);
114 CGBasicMatrixFunctor<T> P_functor(*P);
115 return this->solve(A_functor, b, x, P == nullptr ? nullptr : &P_functor);
116}
117
118template <typename T>
121 Functor const* P)
122{
123 if (x == nullptr)
124 throw std::invalid_argument("RHS must not be null");
125
126 if (A.output_size() != b.size())
127 {
128 this->status.info = CG_INVALID_INPUT;
129 return this->status;
130 }
131
132 /* Set intial x = 0. */
133 if (x->size() != A.input_size())
134 {
135 x->clear();
136 x->resize(A.input_size(), T(0));
137 }
138 else
139 {
140 x->fill(T(0));
141 }
142
143 /* Initial residual is r = b - Ax with x = 0. */
144 Vector r = b;
145
146 /* Regular search direction. */
147 Vector d;
148 /* Preconditioned search direction. */
149 Vector z;
150 /* Norm of residual. */
151 T r_dot_r;
152
153 /* Compute initial search direction. */
154 if (P == nullptr)
155 {
156 r_dot_r = r.dot(r);
157 d = r;
158 }
159 else
160 {
161 z = (*P).multiply(r);
162 r_dot_r = z.dot(r);
163 d = z;
164 }
165
166 for (this->status.num_iterations = 0;
167 this->status.num_iterations < this->opts.max_iterations;
168 this->status.num_iterations += 1)
169 {
170 /* Compute step size in search direction. */
171 Vector Ad = A.multiply(d);
172 T alpha = r_dot_r / d.dot(Ad);
173
174 /* Update parameter vector. */
175 *x = (*x).add(d.multiply(alpha));
176
177 /* Compute new residual and its norm. */
178 r = r.subtract(Ad.multiply(alpha));
179 T new_r_dot_r = r.dot(r);
180
181 /* Check tolerance condition. */
182 if (new_r_dot_r < this->opts.tolerance)
183 {
184 this->status.info = CG_CONVERGENCE;
185 return this->status;
186 }
187
188 /* Precondition residual if necessary. */
189 if (P != nullptr)
190 {
191 z = P->multiply(r);
192 new_r_dot_r = z.dot(r);
193 }
194
195 /*
196 * Update search direction.
197 * The next residual will be orthogonal to new Krylov space.
198 */
199 T beta = new_r_dot_r / r_dot_r;
200 if (P != nullptr)
201 d = z.add(d.multiply(beta));
202 else
203 d = r.add(d.multiply(beta));
204
205 /* Update residual norm. */
206 r_dot_r = new_r_dot_r;
207 }
208
209 this->status.info = CG_MAX_ITERATIONS;
210 return this->status;
211}
212
213/* ---------------------------------------------------------------- */
214
215template <typename T>
216inline
221
222template <typename T>
223inline DenseVector<T>
225{
226 return this->A->multiply(x);
227}
228
229template <typename T>
230inline std::size_t
232{
233 return A->num_cols();
234}
235
236template <typename T>
237inline std::size_t
239{
240 return A->num_rows();
241}
242
245
246#endif /* SFM_BA_CONJUGATE_GRADIENT_HEADER */
CGBasicMatrixFunctor(SparseMatrix< T > const &A)
DenseVector< T > multiply(DenseVector< T > const &x) const
std::size_t output_size(void) const
virtual std::size_t output_size(void) const =0
virtual Vector multiply(Vector const &x) const =0
virtual std::size_t input_size(void) const =0
ConjugateGradient(Options const &opts)
Status solve(Matrix const &A, Vector const &b, Vector *x, Matrix const *P=nullptr)
DenseVector add(DenseVector const &rhs) const
DenseVector multiply(T const &factor) const
void fill(T const &value)
DenseVector subtract(DenseVector const &rhs) const
void resize(std::size_t size, T const &value=T(0))
std::size_t size(void) const
T dot(DenseVector const &rhs) const
Sparse matrix class in Yale format for column-major matrices.
#define SFM_BA_NAMESPACE_BEGIN
Definition defines.h:22
#define SFM_NAMESPACE_END
Definition defines.h:14
#define SFM_NAMESPACE_BEGIN
Definition defines.h:13
#define SFM_BA_NAMESPACE_END
Definition defines.h:23