Blender  V3.3
ConstrainedConjugateGradient.h
Go to the documentation of this file.
1 /* SPDX-License-Identifier: GPL-2.0-or-later
2  * Copyright Blender Foundation. All rights reserved. */
3 
4 #pragma once
5 
10 #include <Eigen/Core>
11 
12 namespace Eigen {
13 
14 namespace internal {
15 
27 template<typename MatrixType,
28  typename Rhs,
29  typename Dest,
30  typename FilterMatrixType,
31  typename Preconditioner>
32 EIGEN_DONT_INLINE void constrained_conjugate_gradient(const MatrixType &mat,
33  const Rhs &rhs,
34  Dest &x,
35  const FilterMatrixType &filter,
36  const Preconditioner &precond,
37  int &iters,
38  typename Dest::RealScalar &tol_error)
39 {
40  using std::abs;
41  using std::sqrt;
42  typedef typename Dest::RealScalar RealScalar;
43  typedef typename Dest::Scalar Scalar;
44  typedef Matrix<Scalar, Dynamic, 1> VectorType;
45 
46  RealScalar tol = tol_error;
47  int maxIters = iters;
48 
49  int n = mat.cols();
50 
51  VectorType residual = filter * (rhs - mat * x); /* initial residual */
52 
53  RealScalar rhsNorm2 = (filter * rhs).squaredNorm();
54  if (rhsNorm2 == 0) {
55  /* XXX TODO: set constrained result here. */
56  x.setZero();
57  iters = 0;
58  tol_error = 0;
59  return;
60  }
61  RealScalar threshold = tol * tol * rhsNorm2;
62  RealScalar residualNorm2 = residual.squaredNorm();
63  if (residualNorm2 < threshold) {
64  iters = 0;
65  tol_error = sqrt(residualNorm2 / rhsNorm2);
66  return;
67  }
68 
69  VectorType p(n);
70  p = filter * precond.solve(residual); /* initial search direction */
71 
72  VectorType z(n), tmp(n);
73  RealScalar absNew = numext::real(
74  residual.dot(p)); /* the square of the absolute value of r scaled by invM */
75  int i = 0;
76  while (i < maxIters) {
77  tmp.noalias() = filter * (mat * p); /* the bottleneck of the algorithm */
78 
79  Scalar alpha = absNew / p.dot(tmp); /* the amount we travel on dir */
80  x += alpha * p; /* update solution */
81  residual -= alpha * tmp; /* update residue */
82 
83  residualNorm2 = residual.squaredNorm();
84  if (residualNorm2 < threshold) {
85  break;
86  }
87 
88  z = precond.solve(residual); /* approximately solve for "A z = residual" */
89 
90  RealScalar absOld = absNew;
91  absNew = numext::real(residual.dot(z)); /* update the absolute value of r */
92  RealScalar beta =
93  absNew /
94  absOld; /* calculate the Gram-Schmidt value used to create the new search direction */
95  p = filter * (z + beta * p); /* update search direction */
96  i++;
97  }
98  tol_error = sqrt(residualNorm2 / rhsNorm2);
99  iters = i;
100 }
101 
102 } // namespace internal
103 
104 #if 0 /* unused */
105 template<typename MatrixType> struct MatrixFilter {
106  MatrixFilter() : m_cmat(NULL)
107  {
108  }
109 
110  MatrixFilter(const MatrixType &cmat) : m_cmat(&cmat)
111  {
112  }
113 
114  void setMatrix(const MatrixType &cmat)
115  {
116  m_cmat = &cmat;
117  }
118 
119  template<typename VectorType> void apply(VectorType v) const
120  {
121  v = (*m_cmat) * v;
122  }
123 
124  protected:
125  const MatrixType *m_cmat;
126 };
127 #endif
128 
129 template<typename _MatrixType,
130  int _UpLo = Lower,
131  typename _FilterMatrixType = _MatrixType,
132  typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar>>
133 class ConstrainedConjugateGradient;
134 
135 namespace internal {
136 
137 template<typename _MatrixType, int _UpLo, typename _FilterMatrixType, typename _Preconditioner>
138 struct traits<
139  ConstrainedConjugateGradient<_MatrixType, _UpLo, _FilterMatrixType, _Preconditioner>> {
140  typedef _MatrixType MatrixType;
141  typedef _FilterMatrixType FilterMatrixType;
142  typedef _Preconditioner Preconditioner;
143 };
144 
145 } // namespace internal
146 
196 template<typename _MatrixType, int _UpLo, typename _FilterMatrixType, typename _Preconditioner>
198  : public IterativeSolverBase<
199  ConstrainedConjugateGradient<_MatrixType, _UpLo, _FilterMatrixType, _Preconditioner>> {
200  typedef IterativeSolverBase<ConstrainedConjugateGradient> Base;
201  using Base::m_error;
202  using Base::m_info;
203  using Base::m_isInitialized;
204  using Base::m_iterations;
205  using Base::mp_matrix;
206 
207  public:
208  typedef _MatrixType MatrixType;
209  typedef typename MatrixType::Scalar Scalar;
210  typedef typename MatrixType::Index Index;
211  typedef typename MatrixType::RealScalar RealScalar;
212  typedef _FilterMatrixType FilterMatrixType;
213  typedef _Preconditioner Preconditioner;
214 
215  enum { UpLo = _UpLo };
216 
217  public:
220  {
221  }
222 
234  {
235  }
236 
238  {
239  }
240 
242  {
243  return m_filter;
244  }
245  const FilterMatrixType &filter() const
246  {
247  return m_filter;
248  }
249 
255  template<typename Rhs, typename Guess>
256  inline const internal::solve_retval_with_guess<ConstrainedConjugateGradient, Rhs, Guess>
257  solveWithGuess(const MatrixBase<Rhs> &b, const Guess &x0) const
258  {
259  eigen_assert(m_isInitialized && "ConjugateGradient is not initialized.");
260  eigen_assert(
261  Base::rows() == b.rows() &&
262  "ConjugateGradient::solve(): invalid number of rows of the right hand side matrix b");
263  return internal::solve_retval_with_guess<ConstrainedConjugateGradient, Rhs, Guess>(
264  *this, b.derived(), x0);
265  }
266 
268  template<typename Rhs, typename Dest> void _solveWithGuess(const Rhs &b, Dest &x) const
269  {
270  m_iterations = Base::maxIterations();
271  m_error = Base::m_tolerance;
272 
273  for (int j = 0; j < b.cols(); j++) {
274  m_iterations = Base::maxIterations();
275  m_error = Base::m_tolerance;
276 
277  typename Dest::ColXpr xj(x, j);
278  internal::constrained_conjugate_gradient(mp_matrix->template selfadjointView<UpLo>(),
279  b.col(j),
280  xj,
281  m_filter,
282  Base::m_preconditioner,
283  m_iterations,
284  m_error);
285  }
286 
287  m_isInitialized = true;
288  m_info = m_error <= Base::m_tolerance ? Success : NoConvergence;
289  }
290 
292  template<typename Rhs, typename Dest> void _solve(const Rhs &b, Dest &x) const
293  {
294  x.setOnes();
295  _solveWithGuess(b, x);
296  }
297 
298  protected:
300 };
301 
302 namespace internal {
303 
304 template<typename _MatrixType, int _UpLo, typename _Filter, typename _Preconditioner, typename Rhs>
305 struct solve_retval<ConstrainedConjugateGradient<_MatrixType, _UpLo, _Filter, _Preconditioner>,
306  Rhs>
307  : solve_retval_base<ConstrainedConjugateGradient<_MatrixType, _UpLo, _Filter, _Preconditioner>,
308  Rhs> {
310  EIGEN_MAKE_SOLVE_HELPERS(Dec, Rhs)
311 
312  template<typename Dest> void evalTo(Dest &dst) const
313  {
314  dec()._solve(rhs(), dst);
315  }
316 };
317 
318 } // end namespace internal
319 
320 } // end namespace Eigen
sqrt(x)+1/max(0
_GL_VOID GLfloat value _GL_VOID_RET _GL_VOID const GLuint GLboolean *residences _GL_BOOL_RET _GL_VOID GLsizei GLfloat GLfloat GLfloat GLfloat const GLubyte *bitmap _GL_VOID_RET _GL_VOID GLenum const void *lists _GL_VOID_RET _GL_VOID const GLdouble *equation _GL_VOID_RET _GL_VOID GLdouble GLdouble blue _GL_VOID_RET _GL_VOID GLfloat GLfloat blue _GL_VOID_RET _GL_VOID GLint GLint blue _GL_VOID_RET _GL_VOID GLshort GLshort blue _GL_VOID_RET _GL_VOID GLubyte GLubyte blue _GL_VOID_RET _GL_VOID GLuint GLuint blue _GL_VOID_RET _GL_VOID GLushort GLushort blue _GL_VOID_RET _GL_VOID GLbyte GLbyte GLbyte alpha _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble alpha _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat alpha _GL_VOID_RET _GL_VOID GLint GLint GLint alpha _GL_VOID_RET _GL_VOID GLshort GLshort GLshort alpha _GL_VOID_RET _GL_VOID GLubyte GLubyte GLubyte alpha _GL_VOID_RET _GL_VOID GLuint GLuint GLuint alpha _GL_VOID_RET _GL_VOID GLushort GLushort GLushort alpha _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLint GLsizei GLsizei GLenum type _GL_VOID_RET _GL_VOID GLsizei GLenum GLenum const void *pixels _GL_VOID_RET _GL_VOID const void *pointer _GL_VOID_RET _GL_VOID GLdouble v _GL_VOID_RET _GL_VOID GLfloat v _GL_VOID_RET _GL_VOID GLint GLint i2 _GL_VOID_RET _GL_VOID GLint j _GL_VOID_RET _GL_VOID GLfloat param _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble GLdouble GLdouble zFar _GL_VOID_RET _GL_UINT GLdouble *equation _GL_VOID_RET _GL_VOID GLenum GLint *params _GL_VOID_RET _GL_VOID GLenum GLfloat *v _GL_VOID_RET _GL_VOID GLenum GLfloat *params _GL_VOID_RET _GL_VOID GLfloat *values _GL_VOID_RET _GL_VOID GLushort *values _GL_VOID_RET _GL_VOID GLenum GLfloat *params _GL_VOID_RET _GL_VOID GLenum GLdouble *params _GL_VOID_RET _GL_VOID GLenum GLint *params _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_BOOL GLfloat param _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID GLenum GLfloat param _GL_VOID_RET _GL_VOID GLenum GLint param _GL_VOID_RET _GL_VOID GLushort pattern _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLint const GLdouble *points _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLint GLdouble GLdouble GLint GLint const GLdouble *points _GL_VOID_RET _GL_VOID GLdouble GLdouble u2 _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLdouble GLdouble v2 _GL_VOID_RET _GL_VOID GLenum GLfloat param _GL_VOID_RET _GL_VOID GLenum GLint param _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLdouble GLdouble nz _GL_VOID_RET _GL_VOID GLfloat GLfloat nz _GL_VOID_RET _GL_VOID GLint GLint nz _GL_VOID_RET _GL_VOID GLshort GLshort nz _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_VOID GLsizei const GLfloat *values _GL_VOID_RET _GL_VOID GLsizei const GLushort *values _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID const GLuint const GLclampf *priorities _GL_VOID_RET _GL_VOID GLdouble y _GL_VOID_RET _GL_VOID GLfloat y _GL_VOID_RET _GL_VOID GLint y _GL_VOID_RET _GL_VOID GLshort y _GL_VOID_RET _GL_VOID GLdouble GLdouble z _GL_VOID_RET _GL_VOID GLfloat GLfloat z _GL_VOID_RET _GL_VOID GLint GLint z _GL_VOID_RET _GL_VOID GLshort GLshort z _GL_VOID_RET _GL_VOID GLdouble GLdouble z
ATTR_WARN_UNUSED_RESULT const BMVert * v
#define A
A conjugate gradient solver for sparse self-adjoint problems with additional constraints.
void _solve(const Rhs &b, Dest &x) const
void _solveWithGuess(const Rhs &b, Dest &x) const
const internal::solve_retval_with_guess< ConstrainedConjugateGradient, Rhs, Guess > solveWithGuess(const MatrixBase< Rhs > &b, const Guess &x0) const
float Scalar
Definition: eigen_utils.h:26
DO_INLINE void filter(lfVector *V, fmatrix3x3 *S)
ccl_gpu_kernel_postfix ccl_global float int int int int float threshold
EIGEN_DONT_INLINE void constrained_conjugate_gradient(const MatrixType &mat, const Rhs &rhs, Dest &x, const FilterMatrixType &filter, const Preconditioner &precond, int &iters, typename Dest::RealScalar &tol_error)
double real
Definition: Precision.h:12
T abs(const T &a)
static const pxr::TfToken b("b", pxr::TfToken::Immortal)
ccl_device_inline float beta(float x, float y)
Definition: util/math.h:775