-
Notifications
You must be signed in to change notification settings - Fork 1
/
conjugateGradient.hpp
executable file
·167 lines (135 loc) · 5.5 KB
/
conjugateGradient.hpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
//
// conjugateGradient.hpp
// conjugateGradient
//
// Created by Adithya Vijaykumar on 08/02/2019.
// Copyright © 2019 Adithya Vijaykumar. All rights reserved.
//
#ifndef conjugateGradient_hpp
#define conjugateGradient_hpp
namespace Eigen {
namespace internal {
template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
bool conjugategradients(const MatrixType& mat, const Rhs& rhs, Dest& x,
const Preconditioner& precond, Index& iters,
typename Dest::RealScalar& tol_error)
{
using std::sqrt;
using std::abs;
typedef typename Dest::RealScalar RealScalar;
typedef typename Dest::Scalar Scalar;
typedef Matrix<Scalar,Dynamic,1> VectorType;
RealScalar tol = tol_error;
Index maxIters = iters;
VectorType r = rhs - mat * x;
VectorType d = r;
RealScalar deltaNew = r.squaredNorm();
RealScalar delta0 = deltaNew;
RealScalar rhs_sqnorm = rhs.squaredNorm();
if(rhs_sqnorm == 0)
{
x.setZero();
return true;
}
RealScalar tol2 = tol*tol*delta0;
Index i = 0;
while ( r.squaredNorm() > tol2 && i<iters )
{
std::cout << x(0) << " " << x(1) << std::endl;
VectorType q = mat * d;
RealScalar alpha = deltaNew/d.dot(q);
x = x + alpha*d;
if (i%50 == 0)
r = rhs - mat * x;
else
r = r - alpha * q;
RealScalar deltaOld = deltaNew;
deltaNew = r.squaredNorm();
RealScalar beta = deltaNew/deltaOld;
d = r+beta*d;
++i;
}
tol_error = sqrt(r.squaredNorm()/rhs_sqnorm);
iters = i;
return true;
}
}
template< typename _MatrixType,
typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> >
class ConjugateGradients;
namespace internal {
template< typename _MatrixType, typename _Preconditioner>
struct traits<ConjugateGradients<_MatrixType,_Preconditioner> >
{
typedef _MatrixType MatrixType;
typedef _Preconditioner Preconditioner;
};
}
template< typename _MatrixType, typename _Preconditioner>
class ConjugateGradients : public IterativeSolverBase<ConjugateGradients<_MatrixType,_Preconditioner> >
{
typedef IterativeSolverBase<ConjugateGradients> Base;
using Base::matrix;
using Base::m_error;
using Base::m_iterations;
using Base::m_info;
using Base::m_isInitialized;
public:
typedef _MatrixType MatrixType;
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
typedef _Preconditioner Preconditioner;
public:
/** Default constructor. */
ConjugateGradients() : Base() {}
/** Initialize the solver with matrix \a A for further \c Ax=b solving.
*
* This constructor is a shortcut for the default constructor followed
* by a call to compute().
*
* \warning this class stores a reference to the matrix A as well as some
* precomputed values that depend on it. Therefore, if \a A is changed
* this class becomes invalid. Call compute() to update it with the new
* matrix A, or modify a copy of A.
*/
template<typename MatrixDerived>
explicit ConjugateGradients(const EigenBase<MatrixDerived>& A) : Base(A.derived()) {}
~ConjugateGradients() {}
/** \internal */
/** Loops over the number of columns of b and does the following:
1. sets the tolerence and maxIterations
2. Calls the function that has the core solver routine
*/
template<typename Rhs,typename Dest>
void _solve_with_guess_impl(const Rhs& b, Dest& x) const
{
bool failed = false;
for(Index j=0; j<b.cols(); ++j)
{
m_iterations = Base::maxIterations();
//******************MANUALLY SET NUM ITERATIONS
//m_iterations = 30;
m_error = Base::m_tolerance;
typename Dest::ColXpr xj(x,j);
if(!internal::conjugategradients(matrix(), b.col(j), xj, Base::m_preconditioner, m_iterations, m_error))
failed = true;
}
m_info = failed ? NumericalIssue
: m_error <= Base::m_tolerance ? Success
: NoConvergence;
m_isInitialized = true;
}
/** \internal */
/** Resizes the x vector to match the dimenstion of b and sets the elements to zero*/
using Base::_solve_impl;
template<typename Rhs,typename Dest>
void _solve_impl(const MatrixBase<Rhs>& b, Dest& x) const
{
x.resize(this->rows(),b.cols());
x.setZero();
_solve_with_guess_impl(b,x);
}
protected:
};
}
#endif /* conjugateGradient_hpp */