Free Electron
BlockPCG.h
Go to the documentation of this file.
1 /* Copyright (C) 2003-2021 Free Electron Organization
2  Any use of this software requires a license. If a valid license
3  was not distributed with this file, visit freeelectron.org. */
4 
5 /** @file */
6 
7 #ifndef __BlockPCG_h__
8 #define __BlockPCG_h__
9 
10 #if 0
11 #include <boost/numeric/ublas/vector.hpp>
12 #include <boost/numeric/ublas/vector_proxy.hpp>
13 #include <boost/numeric/ublas/matrix.hpp>
14 #include <boost/numeric/ublas/triangular.hpp>
15 #include <boost/numeric/ublas/lu.hpp>
16 #include <boost/numeric/ublas/io.hpp>
17 
18 namespace ublas = boost::numeric::ublas;
19 #endif
20 
21 namespace fe
22 {
23 namespace ext
24 {
25 
26 template <typename T>
27 class FE_DL_EXPORT Preconditioner
28 {
29  public:
30 typedef T t_real;
31 typedef Vector<3, t_real> t_v3;
32 typedef Matrix<3,3,t_real> t_matrix;
33 virtual ~Preconditioner(void){}
34 virtual void apply( std::vector<t_v3> &a_r,
35  const std::vector<t_v3> &a_x) const = 0;
36 };
37 
38 template <typename T>
39 class FE_DL_EXPORT BlockDiagonalPreconditioner : public Preconditioner<T>
40 {
41  public:
42 typedef T t_real;
43 typedef Vector<3, t_real> t_v3;
44 typedef Matrix<3,3,t_real> t_matrix;
45  BlockDiagonalPreconditioner(void);
46 virtual ~BlockDiagonalPreconditioner(void);
47 
48 virtual void apply( std::vector<t_v3> &a_r,
49  const std::vector<t_v3> &a_x) const;
50 
51  std::vector<t_matrix> &diagonal(void);
52 
53  private:
54  std::vector<t_matrix> m_diagonal;
55 };
56 
57 #if 0
58 class FE_DL_EXPORT OverlappedInversionPreconditioner : public Preconditioner
59 {
60  public:
61  OverlappedInversionPreconditioner(void){}
62 virtual ~OverlappedInversionPreconditioner(void){}
63 
64 virtual void apply( std::vector<t_v3> &a_r,
65  const std::vector<t_v3> &a_x) const;
66 
67  class Block
68  {
69  public:
70  ublas::matrix<t_real> m_inverse;
71  unsigned int m_L;
72  unsigned int m_H;
73  unsigned int m_m0;
74  unsigned int m_m1;
75  };
76 
77  std::vector<Block> &block(void) { return m_block; }
78 
79  private:
80  std::vector<Block> m_block;
81 };
82 #endif
83 
84 template <typename T>
85 class FE_DL_EXPORT FilterConstraint
86 {
87  public:
88 typedef T t_real;
89 typedef Vector<3, t_real> t_v3;
90 typedef Matrix<3,3,t_real> t_matrix;
91  FilterConstraint(void);
92 virtual ~FilterConstraint(void);
93 
94  const t_matrix &filter(void) const;
95  void setIndex(unsigned long a_index);
96  unsigned long index(void) const;
97 
98  void makePointConstraint(void);
99  void makeLineConstraint(const t_v3 &a_direction);
100  void makePlaneConstraint(const t_v3 &a_normal);
101  void makeInequalityConstraint(const t_v3 &a_axis, t_real a_delta);
102 
103  const t_v3 &axis(void) const {return m_axis;}
104  t_real delta(void) const;
105  bool isInequality(void) const;
106 
107  bool m_inequality;
108  private:
109  t_matrix m_filter;
110  unsigned long m_index;
111  t_real m_delta;
112  t_v3 m_axis;
113 };
114 
115 
116 template <typename T>
117 class FE_DL_EXPORT BlockPCG
118 {
119  public:
120 typedef T t_real;
121 typedef Vector<3, t_real> t_v3;
122 typedef Matrix<3,3,t_real> t_matrix;
123  BlockPCG(void);
124 virtual ~BlockPCG(void);
125 
126  bool solve( std::vector<t_v3> &a_x,
127  const sp< SparseMatrix3x3<t_real> > a_A,
128  const std::vector<t_v3> &a_b,
129  const Preconditioner<T> &a_P,
130  std::vector< FilterConstraint<T> > &a_filters);
131 
132  const std::vector<t_v3> &residual(void) const;
133 
134  unsigned int iterations(void) const { return m_iter; }
135 
136  private:
137  t_real inner( const std::vector<t_v3> &a_left,
138  const std::vector<t_v3> &a_right);
139 
140  void residual( std::vector<t_v3> &a_r,
141  const sp< SparseMatrix3x3<t_real> > a_A,
142  const std::vector<t_v3> &a_x,
143  const std::vector<t_v3> &a_b);
144 
145  void setSize( unsigned long a_size);
146 
147  void applyFilters(std::vector<t_v3> &a_r,
148  std::vector< FilterConstraint<T> > &a_filters);
149 
150  bool m_restart;
151  t_real m_tol;
152  unsigned int m_maxIterations;
153  unsigned int m_iter;
154  std::vector<t_v3> m_r;
155  std::vector<t_v3> m_q;
156  std::vector<t_v3> m_d;
157  std::vector<t_v3> m_s;
158  std::vector<t_v3> m_tmp0;
159  std::vector<t_v3> m_tmp1;
160  unsigned int m_n;
161  std::vector<t_v3> *m_x;
162 
163 };
164 
165 template <typename T>
166 inline
167 BlockPCG<T>::BlockPCG(void)
168 {
169  m_tol = fe::tol;
170 // m_tol = 1.0e-10;
171  m_maxIterations = 100000;
172  m_iter = 0;
173  m_n = 0;
174  m_restart = false;
175 }
176 
177 template <typename T>
178 inline
179 BlockPCG<T>::~BlockPCG(void)
180 {
181 }
182 
183 template <typename T>
184 inline
185 T BlockPCG<T>::inner(const std::vector<t_v3> &a_left, const std::vector<t_v3> &a_right)
186 {
187  t_real s = 0.0;
188  unsigned int n = a_left.size();
189  for(unsigned int i = 0; i < n; i++)
190  {
191  s += dot(a_left[i], a_right[i]);
192  }
193  return s;
194 }
195 
196 template <typename T>
197 inline
198 void BlockPCG<T>::residual(std::vector<t_v3> &a_r, sp< SparseMatrix3x3<t_real> > a_A, const std::vector<t_v3> &a_x, const std::vector<t_v3> &a_b)
199 {
200  unsigned int n = a_b.size();
201 
202  a_A->multiply(a_r, a_x);
203 
204  for(unsigned int i = 0; i < n; i++)
205  {
206  a_r[i] = a_b[i] - a_r[i];
207  }
208 }
209 
210 template <typename T>
211 inline
212 void BlockPCG<T>::setSize(unsigned long a_size)
213 {
214  if(a_size != m_n)
215  {
216  m_n = a_size;
217  m_r.resize(m_n);
218  m_q.resize(m_n);
219  m_d.resize(m_n);
220  m_s.resize(m_n);
221  m_tmp0.resize(m_n);
222  m_tmp1.resize(m_n);
223  }
224 }
225 
226 template <typename T>
227 inline
228 void BlockPCG<T>::applyFilters(std::vector<t_v3> &a_r, std::vector< FilterConstraint<T> > &a_filters)
229 {
230  for(unsigned int i = 0; i < a_filters.size(); i++)
231  {
232  unsigned long index = a_filters[i].index();
233 
234 // TODO: probably remove this inequality stuff
235  if(a_filters[i].isInequality())
236  {
237  t_v3 &v = (*m_x)[index];
238  t_real d = dot(a_filters[i].axis(), v);
239 
240  d = (*m_x)[index][2];
241 
242 //fe_fprintf(stderr, "index %d delta %f d %f\n", index, a_filters[i].delta(), d);
243  if(d < a_filters[i].delta())
244  {
245 //fe_fprintf(stderr, "<***********************\n");
246  //(*m_x)[index] = (*m_x)[index] - a_filters[i].axis() * d;
247  //(*m_x)[index] = (*m_x)[index] + a_filters[i].axis() * a_filters[i].delta();
248  (*m_x)[index][2] = a_filters[i].delta();
249  a_r[index] = fe::multiply(a_filters[i].filter(), a_r[index]);
250  a_filters[i].m_inequality = false;
251  m_restart = true;
252  }
253  }
254  else
255  {
256  a_r[index] = fe::multiply(a_filters[i].filter(), a_r[index]);
257  }
258  }
259 //fe_fprintf(stderr, "-----\n");
260 }
261 
262 template <typename T>
263 inline
264 const std::vector< Vector<3,T> > &BlockPCG<T>::residual(void) const
265 {
266  return m_r;
267 }
268 
269 //#define FE_BLOCKPCG_RESTART
270 template <typename T>
271 inline
272 bool BlockPCG<T>::solve(std::vector<t_v3> &a_x, const sp< SparseMatrix3x3<t_real> > a_A, const std::vector<t_v3> &a_b, const Preconditioner<T> &a_P, std::vector< FilterConstraint<T> > &a_filters)
273 {
274  bool converged = false;
275 
276  setSize(a_b.size());
277 
278  t_real delta0, deltaNew, deltaOld, deltaStop;
279 
280 #ifdef FE_BLOCKPCG_RESTART
281  unsigned int modn = 50;
282 #endif
283 
284  m_x = &(a_x);
285 
286  m_tmp0 = a_b;
287  applyFilters(m_tmp0, a_filters);
288  a_P.apply(m_tmp1, m_tmp0);
289  delta0 = inner(m_tmp0, m_tmp1);
290 
291  residual(m_r, a_A, a_x, a_b);
292  applyFilters(m_r, a_filters);
293 
294  a_P.apply(m_d, m_r);
295  applyFilters(m_d, a_filters);
296 
297  deltaNew = inner(m_r,m_d);
298  deltaStop = m_tol * m_tol * delta0;
299 
300 
301  m_iter = 0;
302  while((m_iter < m_maxIterations) && (deltaNew > deltaStop))
303  {
304  a_A->multiply(m_q, m_d);
305  applyFilters(m_q, a_filters);
306 
307  t_real alpha = deltaNew / inner(m_d, m_q);
308 
309  for(unsigned int i = 0; i < m_n; i++)
310  {
311  a_x[i] += alpha * m_d[i];
312  }
313 
314 #if 0
315  if(m_restart)
316  {
317  m_restart = false;
318  residual(m_r, a_A, a_x, a_b);
319  applyFilters(m_r, a_filters);
320  }
321  else
322  {
323  for(unsigned int i = 0; i < m_n; i++)
324  {
325  m_r[i] -= alpha * m_q[i];
326  }
327  }
328 #endif
329  for(unsigned int i = 0; i < m_n; i++)
330  {
331  m_r[i] -= alpha * m_q[i];
332  }
333 
334 #if 0
335  //if(!(m_iter % 1))
336  if(m_restart)
337  {
338 fe_fprintf(stderr, "CONV HARD RESTART %d\n", m_iter);
339  m_restart = false;
340  m_tmp0 = a_b;
341  applyFilters(m_tmp0, a_filters);
342  a_P.apply(m_tmp1, m_tmp0);
343  delta0 = inner(m_tmp0, m_tmp1);
344 
345  residual(m_r, a_A, a_x, a_b);
346  applyFilters(m_r, a_filters);
347 
348  a_P.apply(m_d, m_r);
349  applyFilters(m_d, a_filters);
350 
351  deltaNew = inner(m_r,m_d);
352  deltaStop = m_tol * m_tol * delta0;
353 
354  a_A->multiply(m_q, m_d);
355  applyFilters(m_q, a_filters);
356 
357  t_real alpha = deltaNew / inner(m_d, m_q);
358 
359  for(unsigned int i = 0; i < m_n; i++)
360  {
361  a_x[i] += alpha * m_d[i];
362  }
363  }
364  else
365  {
366  for(unsigned int i = 0; i < m_n; i++)
367  {
368  m_r[i] -= alpha * m_q[i];
369  }
370  }
371 #endif
372 
373 #if 0
374 #ifdef FE_BLOCKPCG_RESTART
375  if(!(m_iter % modn))
376  {
377  residual(m_r, a_A, a_x, a_b);
378  applyFilters(m_r, a_filters);
379  }
380  else
381  {
382 #endif
383  for(unsigned int i = 0; i < m_n; i++)
384  {
385  m_r[i] -= alpha * m_q[i];
386  }
387 #ifdef FE_BLOCKPCG_RESTART
388  }
389 #endif
390 #endif
391 
392  a_P.apply(m_s, m_r);
393 
394  deltaOld = deltaNew;
395 
396  deltaNew = inner(m_r, m_s);
397 
398  t_real beta = deltaNew/deltaOld;
399 
400  for(unsigned int i = 0; i < m_n; i++)
401  {
402  m_d[i] = m_s[i] + beta*m_d[i];
403  }
404 
405  applyFilters(m_d, a_filters);
406 
407  m_iter++;
408  }
409 
410  if(m_iter < m_maxIterations)
411  {
412  converged = true;
413  }
414  residual(m_r, a_A, a_x, a_b);
415 
416  //feLogGroup("BlockPCG", "CONVERGED IN %d ITERATIONS\n", m_iter);
417 
418  return converged;
419 }
420 
421 template <typename T>
422 inline
423 BlockDiagonalPreconditioner<T>::BlockDiagonalPreconditioner(void)
424 {
425 }
426 
427 template <typename T>
428 inline
429 BlockDiagonalPreconditioner<T>::~BlockDiagonalPreconditioner(void)
430 {
431 }
432 
433 template <typename T>
434 inline
435 std::vector< Matrix<3,3,T> > &BlockDiagonalPreconditioner<T>::diagonal(void)
436 {
437  return m_diagonal;
438 }
439 
440 template <typename T>
441 inline
442 void BlockDiagonalPreconditioner<T>::apply(std::vector<t_v3> &a_r, const std::vector<t_v3> &a_x) const
443 {
444  unsigned int n = a_r.size();
445  for(unsigned int i = 0; i < n; i++)
446  {
447  multiply(a_r[i], m_diagonal[i], a_x[i]);
448 //fe_fprintf(stderr, "apply P %f %f %f\n", a_x[i][0], a_x[i][1], a_x[i][2]);
449 //fe_fprintf(stderr, "apply P %f %f %f\n", a_r[i][0], a_r[i][1], a_r[i][2]);
450  }
451 }
452 
453 #if 0
454 void OverlappedInversionPreconditioner::apply(std::vector<t_v3> &a_r, const std::vector<t_v3> &a_x) const
455 {
456 fe_fprintf(stderr, ".");
457  for(unsigned int i = 0; i < a_r.size(); i++)
458  {
459  a_r[i] = a_x[i];
460  }
461  unsigned int n = m_block.size();
462  for(unsigned int i = 0; i < n; i++)
463  {
464  ublas::vector<t_real> v, vr;
465  unsigned int nn = 3*(m_block[i].m_H - m_block[i].m_L);
466  v.resize(nn);
467  for(unsigned int j = m_block[i].m_L; j < m_block[i].m_H; j++)
468  {
469 //fe_fprintf(stderr, "%d <- %d\n",(j-m_block[i].m_L)*3+0, j);
470  v[(j-m_block[i].m_L)*3+0] = a_x[j][0];
471  v[(j-m_block[i].m_L)*3+1] = a_x[j][1];
472  v[(j-m_block[i].m_L)*3+2] = a_x[j][2];
473  }
474 #if 0
475 fe_fprintf(stderr, ">>>>>>> ");
476 for(unsigned int k = 0; k < v.size(); k++)
477 {
478 fe_fprintf(stderr, "%f ", v[k]);
479 }
480 fe_fprintf(stderr, "\n");
481 std::cerr << m_block[i].m_inverse << "\n";
482 #endif
483  v = prod(m_block[i].m_inverse, v);
484 #if 0
485 fe_fprintf(stderr, "<<<<<<< ");
486 for(unsigned int k = 0; k < v.size(); k++)
487 {
488 fe_fprintf(stderr, "%f ", v[k]);
489 }
490 fe_fprintf(stderr, "\n");
491 #endif
492 
493  //for(unsigned int j = m_block[i].m_m0; j < m_block[i].m_m1; j++)
494  for(unsigned int j = m_block[i].m_L; j < m_block[i].m_H; j++)
495  {
496  if(j < m_block[i].m_m0) continue;
497  if(j >= m_block[i].m_m1) continue;
498 //fe_fprintf(stderr, "%d -> %d\n",j, (j-m_block[i].m_L)*3+0);
499  a_r[j][0] = v[(j-m_block[i].m_L)*3+0];
500  a_r[j][1] = v[(j-m_block[i].m_L)*3+1];
501  a_r[j][2] = v[(j-m_block[i].m_L)*3+2];
502  }
503  }
504 }
505 #endif
506 
507 template <typename T>
508 inline
509 FilterConstraint<T>::FilterConstraint(void)
510 {
511  m_inequality = false;
512 }
513 
514 template <typename T>
515 inline
516 FilterConstraint<T>::~FilterConstraint(void)
517 {
518 }
519 
520 template <typename T>
521 inline
522 void FilterConstraint<T>::setIndex(unsigned long a_index)
523 {
524  m_index = a_index;
525 }
526 
527 template <typename T>
528 inline
529 const Matrix<3,3,T> &FilterConstraint<T>::filter(void) const
530 {
531  return m_filter;
532 }
533 
534 template <typename T>
535 inline
536 unsigned long FilterConstraint<T>::index(void) const
537 {
538  return m_index;
539 }
540 
541 template <typename T>
542 inline
543 void FilterConstraint<T>::makePointConstraint(void)
544 {
545  setAll(m_filter, 0.0);
546  m_inequality = false;
547 }
548 
549 template <typename T>
550 inline
551 void FilterConstraint<T>::makeLineConstraint(const t_v3 &a_direction)
552 {
553  setAll(m_filter, 0.0);
554  outerProduct(m_filter, a_direction, a_direction);
555  m_inequality = false;
556 }
557 
558 template <typename T>
559 inline
560 void FilterConstraint<T>::makePlaneConstraint(const t_v3 &a_normal)
561 {
562 
563  setAll(m_filter, 0.0);
564  m_filter(0,0) = 1.0 - a_normal[0] * a_normal[0];
565  m_filter(0,1) = - a_normal[0] * a_normal[1];
566  m_filter(0,2) = - a_normal[0] * a_normal[2];
567  m_filter(1,0) = m_filter(0,1);
568  m_filter(1,1) = 1.0 - a_normal[1] * a_normal[1];
569  m_filter(1,2) = - a_normal[1] * a_normal[2];
570  m_filter(2,0) = m_filter(0,2);
571  m_filter(2,1) = m_filter(1,2);
572  m_filter(2,2) = 1.0 - a_normal[2] * a_normal[2];
573 
574  m_inequality = false;
575 }
576 
577 template <typename T>
578 inline
579 void FilterConstraint<T>::makeInequalityConstraint(const t_v3 &a_axis, t_real a_delta)
580 {
581  m_axis = a_axis;
582  normalize(m_axis);
583  makePlaneConstraint(m_axis);
584 
585  m_inequality = true;
586  m_delta = a_delta;
587 }
588 
589 template <typename T>
590 inline
591 bool FilterConstraint<T>::isInequality(void) const
592 {
593  return m_inequality;
594 }
595 
596 template <typename T>
597 inline
598 T FilterConstraint<T>::delta(void) const
599 {
600  return m_delta;
601 }
602 
603 
604 } /* namespace ext */
605 } /* namespace fe */
606 
607 #endif /* __BlockPCG_h__ */
608 
fe
kernel
Definition: namespace.dox:3
Generated by   doxygen 1.8.13