GetFEM  5.4.3
gmm_solver_Schwarz_additive.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2002-2020 Yves Renard
5 
6  This file is a part of GetFEM
7 
8  GetFEM is free software; you can redistribute it and/or modify it
9  under the terms of the GNU Lesser General Public License as published
10  by the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version along with the GCC Runtime Library
12  Exception either version 3.1 or (at your option) any later version.
13  This program is distributed in the hope that it will be useful, but
14  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  License and GCC Runtime Library Exception for more details.
17  You should have received a copy of the GNU Lesser General Public License
18  along with this program; if not, write to the Free Software Foundation,
19  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
20 
21  As a special exception, you may use this file as it is a part of a free
22  software library without restriction. Specifically, if other files
23  instantiate templates or use macros or inline functions from this file,
24  or you compile this file and link it with other files to produce an
25  executable, this file does not by itself cause the resulting executable
26  to be covered by the GNU Lesser General Public License. This exception
27  does not however invalidate any other reasons why the executable file
28  might be covered by the GNU Lesser General Public License.
29 
30 ===========================================================================*/
31 
32 /**@file gmm_solver_Schwarz_additive.h
33  @author Yves Renard <[email protected]>
34  @author Michel Fournie <[email protected]>
35  @date October 13, 2002.
36 */
37 
38 #ifndef GMM_SOLVERS_SCHWARZ_ADDITIVE_H__
39 #define GMM_SOLVERS_SCHWARZ_ADDITIVE_H__
40 
41 #include "gmm_kernel.h"
42 #if defined(GMM_USES_SUPERLU)
43 #include "gmm_superlu_interface.h"
44 #else
45 #include "gmm_MUMPS_interface.h"
46 #endif
47 #include "gmm_solver_cg.h"
48 #include "gmm_solver_gmres.h"
49 #include "gmm_solver_bicgstab.h"
50 #include "gmm_solver_qmr.h"
51 
52 namespace gmm {
53 
54  /* ******************************************************************** */
55  /* Additive Schwarz interfaced local solvers */
56  /* ******************************************************************** */
57 
58  struct using_cg {};
59  struct using_gmres {};
60  struct using_bicgstab {};
61  struct using_qmr {};
62 
63  template <typename P, typename local_solver, typename Matrix>
64  struct actual_precond {
65  typedef P APrecond;
66  static const APrecond &transform(const P &PP) { return PP; }
67  };
68 
69  template <typename Matrix1, typename Precond, typename Vector>
70  void AS_local_solve(using_cg, const Matrix1 &A, Vector &x, const Vector &b,
71  const Precond &P, iteration &iter)
72  { cg(A, x, b, P, iter); }
73 
74  template <typename Matrix1, typename Precond, typename Vector>
75  void AS_local_solve(using_gmres, const Matrix1 &A, Vector &x,
76  const Vector &b, const Precond &P, iteration &iter)
77  { gmres(A, x, b, P, 100, iter); }
78 
79  template <typename Matrix1, typename Precond, typename Vector>
80  void AS_local_solve(using_bicgstab, const Matrix1 &A, Vector &x,
81  const Vector &b, const Precond &P, iteration &iter)
82  { bicgstab(A, x, b, P, iter); }
83 
84  template <typename Matrix1, typename Precond, typename Vector>
85  void AS_local_solve(using_qmr, const Matrix1 &A, Vector &x,
86  const Vector &b, const Precond &P, iteration &iter)
87  { qmr(A, x, b, P, iter); }
88 
89 #if defined(GMM_USES_SUPERLU)
90  struct using_superlu {};
91 
92  template <typename P, typename Matrix>
93  struct actual_precond<P, using_superlu, Matrix> {
94  typedef typename linalg_traits<Matrix>::value_type value_type;
95  typedef SuperLU_factor<value_type> APrecond;
96  template <typename PR>
97  static APrecond transform(const PR &) { return APrecond(); }
98  static const APrecond &transform(const APrecond &PP) { return PP; }
99  };
100 
101  template <typename Matrix1, typename Precond, typename Vector>
102  void AS_local_solve(using_superlu, const Matrix1 &, Vector &x,
103  const Vector &b, const Precond &P, iteration &iter)
104  { P.solve(x, b); iter.set_iteration(1); }
105 #endif
106 
107  /* ******************************************************************** */
108  /* Additive Schwarz Linear system */
109  /* ******************************************************************** */
110 
111  template <typename Matrix1, typename Matrix2, typename Precond,
112  typename local_solver>
113  struct add_schwarz_mat{
114  typedef typename linalg_traits<Matrix1>::value_type value_type;
115 
116  const Matrix1 *A;
117  const std::vector<Matrix2> *vB;
118  std::vector<Matrix2> vAloc;
119  mutable iteration iter;
120  double residual;
121  mutable size_type itebilan;
122  mutable std::vector<std::vector<value_type> > gi, fi;
123  std::vector<typename actual_precond<Precond, local_solver,
124  Matrix1>::APrecond> precond1;
125 
126  void init(const Matrix1 &A_, const std::vector<Matrix2> &vB_,
127  iteration iter_, const Precond &P, double residual_);
128 
129  add_schwarz_mat() {}
130  add_schwarz_mat(const Matrix1 &A_, const std::vector<Matrix2> &vB_,
131  iteration iter_, const Precond &P, double residual_)
132  { init(A_, vB_, iter_, P, residual_); }
133  };
134 
135  template <typename Matrix1, typename Matrix2, typename Precond,
136  typename local_solver>
137  void add_schwarz_mat<Matrix1, Matrix2, Precond, local_solver>::init(
138  const Matrix1 &A_, const std::vector<Matrix2> &vB_,
139  iteration iter_, const Precond &P, double residual_) {
140 
141  vB = &vB_; A = &A_; iter = iter_;
142  residual = residual_;
143 
144  size_type nb_sub = vB->size();
145  vAloc.resize(nb_sub);
146  gi.resize(nb_sub); fi.resize(nb_sub);
147  precond1.resize(nb_sub);
148  std::fill(precond1.begin(), precond1.end(),
149  actual_precond<Precond, local_solver, Matrix1>::transform(P));
150  itebilan = 0;
151 
152  if (iter.get_noisy()) cout << "Init pour sub dom ";
153 #ifdef GMM_USES_MPI
154  int size,tranche,borne_sup,borne_inf,rank,tag1=11,tag2=12,tag3=13,sizepr = 0;
155  // int tab[4];
156  double t_ref,t_final;
157  MPI_Status status;
158  t_ref=MPI_Wtime();
159  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
160  MPI_Comm_size(MPI_COMM_WORLD, &size);
161  tranche=nb_sub/size;
162  borne_inf=rank*tranche;
163  borne_sup=(rank+1)*tranche;
164  // if (rank==size-1) borne_sup = nb_sub;
165 
166  cout << "Nombre de sous domaines " << borne_sup - borne_inf << endl;
167 
168  int sizeA = mat_nrows(*A);
169  gmm::csr_matrix<value_type> Acsr(sizeA, sizeA), Acsrtemp(sizeA, sizeA);
170  gmm::copy(gmm::eff_matrix(*A), Acsr);
171  int next = (rank + 1) % size;
172  int previous = (rank + size - 1) % size;
173  //communication of local information on ring pattern
174  //Each process receive Nproc-1 contributions
175 
176  for (int nproc = 0; nproc < size; ++nproc) {
177  for (size_type i = size_type(borne_inf); i < size_type(borne_sup); ++i) {
178 // for (size_type i = 0; i < nb_sub/size; ++i) {
179 // for (size_type i = 0; i < nb_sub; ++i) {
180  // size_type i=(rank+size*(j-1)+nb_sub)%nb_sub;
181 
182  cout << "Sous domaines " << i << " : " << mat_ncols((*vB)[i]) << endl;
183 #else
184  for (size_type i = 0; i < nb_sub; ++i) {
185 #endif
186  if (iter.get_noisy())
187  cout << i << " " << std::flush;
188  Matrix2 Maux(mat_ncols((*vB)[i]), mat_nrows((*vB)[i]));
189 
190 #ifdef GMM_USES_MPI
191  Matrix2 Maux2(mat_ncols((*vB)[i]), mat_ncols((*vB)[i]));
192  if (nproc == 0) {
193  gmm::resize(vAloc[i], mat_ncols((*vB)[i]), mat_ncols((*vB)[i]));
194  gmm::clear(vAloc[i]);
195  }
196  gmm::mult(gmm::transposed((*vB)[i]), Acsr, Maux);
197  gmm::mult(Maux, (*vB)[i], Maux2);
198  gmm::add(Maux2, vAloc[i]);
199 #else
200  gmm::resize(vAloc[i], mat_ncols((*vB)[i]), mat_ncols((*vB)[i]));
201  gmm::mult(gmm::transposed((*vB)[i]), *A, Maux);
202  gmm::mult(Maux, (*vB)[i], vAloc[i]);
203 #endif
204 
205 #ifdef GMM_USES_MPI
206  if (nproc == size - 1 ) {
207 #endif
208  precond1[i].build_with(vAloc[i]);
209  gmm::resize(fi[i], mat_ncols((*vB)[i]));
210  gmm::resize(gi[i], mat_ncols((*vB)[i]));
211 #ifdef GMM_USES_MPI
212  }
213 #else
214  }
215 #endif
216 #ifdef GMM_USES_MPI
217  }
218  if (nproc != size - 1) {
219  MPI_Sendrecv(&(Acsr.jc[0]), sizeA+1, MPI_INT, next, tag2,
220  &(Acsrtemp.jc[0]), sizeA+1, MPI_INT, previous, tag2,
221  MPI_COMM_WORLD, &status);
222  if (Acsrtemp.jc[sizeA] > size_type(sizepr)) {
223  sizepr = Acsrtemp.jc[sizeA];
224  gmm::resize(Acsrtemp.pr, sizepr);
225  gmm::resize(Acsrtemp.ir, sizepr);
226  }
227  MPI_Sendrecv(&(Acsr.ir[0]), Acsr.jc[sizeA], MPI_INT, next, tag1,
228  &(Acsrtemp.ir[0]), Acsrtemp.jc[sizeA], MPI_INT, previous, tag1,
229  MPI_COMM_WORLD, &status);
230 
231  MPI_Sendrecv(&(Acsr.pr[0]), Acsr.jc[sizeA], mpi_type(value_type()), next, tag3,
232  &(Acsrtemp.pr[0]), Acsrtemp.jc[sizeA], mpi_type(value_type()), previous, tag3,
233  MPI_COMM_WORLD, &status);
234  gmm::copy(Acsrtemp, Acsr);
235  }
236  }
237  t_final=MPI_Wtime();
238  cout<<"temps boucle precond "<< t_final-t_ref<<endl;
239 #endif
240  if (iter.get_noisy()) cout << "\n";
241  }
242 
243  template <typename Matrix1, typename Matrix2, typename Precond,
244  typename Vector2, typename Vector3, typename local_solver>
245  void mult(const add_schwarz_mat<Matrix1, Matrix2, Precond, local_solver> &M,
246  const Vector2 &p, Vector3 &q) {
247  size_type itebilan = 0;
248 #ifdef GMM_USES_MPI
249  static double tmult_tot = 0.0;
250  double t_ref = MPI_Wtime();
251 #endif
252  // cout << "tmult AS begin " << endl;
253  mult(*(M.A), p, q);
254 #ifdef GMM_USES_MPI
255  tmult_tot += MPI_Wtime()-t_ref;
256  cout << "tmult_tot = " << tmult_tot << endl;
257 #endif
258  std::vector<double> qbis(gmm::vect_size(q));
259  std::vector<double> qter(gmm::vect_size(q));
260 #ifdef GMM_USES_MPI
261  // MPI_Status status;
262  // MPI_Request request,request1;
263  // int tag=111;
264  int size,tranche,borne_sup,borne_inf,rank;
265  size_type nb_sub=M.fi.size();
266  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
267  MPI_Comm_size(MPI_COMM_WORLD, &size);
268  tranche=nb_sub/size;
269  borne_inf=rank*tranche;
270  borne_sup=(rank+1)*tranche;
271  // if (rank==size-1) borne_sup=nb_sub;
272  // int next = (rank + 1) % size;
273  // int previous = (rank + size - 1) % size;
274  t_ref = MPI_Wtime();
275  for (size_type i = size_type(borne_inf); i < size_type(borne_sup); ++i)
276 // for (size_type i = 0; i < nb_sub/size; ++i)
277  // for (size_type j = 0; j < nb_sub; ++j)
278 #else
279  for (size_type i = 0; i < M.fi.size(); ++i)
280 #endif
281  {
282 #ifdef GMM_USES_MPI
283  // size_type i=j; // (rank+size*(j-1)+nb_sub)%nb_sub;
284 #endif
285  gmm::mult(gmm::transposed((*(M.vB))[i]), q, M.fi[i]);
286  M.iter.init();
287  AS_local_solve(local_solver(), (M.vAloc)[i], (M.gi)[i],
288  (M.fi)[i],(M.precond1)[i],M.iter);
289  itebilan = std::max(itebilan, M.iter.get_iteration());
290  }
291 
292 #ifdef GMM_USES_MPI
293  cout << "First AS loop time " << MPI_Wtime() - t_ref << endl;
294 #endif
295 
296  gmm::clear(q);
297 #ifdef GMM_USES_MPI
298  t_ref = MPI_Wtime();
299  // for (size_type j = 0; j < nb_sub; ++j)
300  for (size_type i = size_type(borne_inf); i < size_type(borne_sup); ++i)
301 
302 #else
303  for (size_type i = 0; i < M.gi.size(); ++i)
304 #endif
305  {
306 
307 #ifdef GMM_USES_MPI
308  // size_type i=j; // (rank+size*(j-1)+nb_sub)%nb_sub;
309 // gmm::mult((*(M.vB))[i], M.gi[i], qbis,qbis);
310  gmm::mult((*(M.vB))[i], M.gi[i], qter);
311  add(qter,qbis,qbis);
312 #else
313  gmm::mult((*(M.vB))[i], M.gi[i], q, q);
314 #endif
315  }
316 #ifdef GMM_USES_MPI
317  //WARNING this add only if you use the ring pattern below
318  // need to do this below if using a n explicit ring pattern communication
319 
320 // add(qbis,q,q);
321  cout << "Second AS loop time " << MPI_Wtime() - t_ref << endl;
322 #endif
323 
324 
325 #ifdef GMM_USES_MPI
326  // int tag1=11;
327  static double t_tot = 0.0;
328  double t_final;
329  t_ref=MPI_Wtime();
330 // int next = (rank + 1) % size;
331 // int previous = (rank + size - 1) % size;
332  //communication of local information on ring pattern
333  //Each process receive Nproc-1 contributions
334 
335 // if (size > 1) {
336 // for (int nproc = 0; nproc < size-1; ++nproc)
337 // {
338 
339 // MPI_Sendrecv(&(qbis[0]), gmm::vect_size(q), MPI_DOUBLE, next, tag1,
340 // &(qter[0]), gmm::vect_size(q),MPI_DOUBLE,previous,tag1,
341 // MPI_COMM_WORLD,&status);
342 // gmm::copy(qter, qbis);
343 // add(qbis,q,q);
344 // }
345 // }
346  MPI_Allreduce(&(qbis[0]), &(q[0]),gmm::vect_size(q), MPI_DOUBLE,
347  MPI_SUM,MPI_COMM_WORLD);
348  t_final=MPI_Wtime();
349  t_tot += t_final-t_ref;
350  cout<<"["<< rank<<"] temps reduce Resol "<< t_final-t_ref << " t_tot = " << t_tot << endl;
351 #endif
352 
353  if (M.iter.get_noisy() > 0) cout << "itebloc = " << itebilan << endl;
354  M.itebilan += itebilan;
355  M.iter.set_resmax((M.iter.get_resmax() + M.residual) * 0.5);
356  }
357 
358  template <typename Matrix1, typename Matrix2, typename Precond,
359  typename Vector2, typename Vector3, typename local_solver>
360  void mult(const add_schwarz_mat<Matrix1, Matrix2, Precond, local_solver> &M,
361  const Vector2 &p, const Vector3 &q) {
362  mult(M, p, const_cast<Vector3 &>(q));
363  }
364 
365  template <typename Matrix1, typename Matrix2, typename Precond,
366  typename Vector2, typename Vector3, typename Vector4,
367  typename local_solver>
368  void mult(const add_schwarz_mat<Matrix1, Matrix2, Precond, local_solver> &M,
369  const Vector2 &p, const Vector3 &p2, Vector4 &q)
370  { mult(M, p, q); add(p2, q); }
371 
372  template <typename Matrix1, typename Matrix2, typename Precond,
373  typename Vector2, typename Vector3, typename Vector4,
374  typename local_solver>
375  void mult(const add_schwarz_mat<Matrix1, Matrix2, Precond, local_solver> &M,
376  const Vector2 &p, const Vector3 &p2, const Vector4 &q)
377  { mult(M, p, const_cast<Vector4 &>(q)); add(p2, q); }
378 
379  /* ******************************************************************** */
380  /* Additive Schwarz interfaced global solvers */
381  /* ******************************************************************** */
382 
383  template <typename ASM_type, typename Vect>
384  void AS_global_solve(using_cg, const ASM_type &ASM, Vect &x,
385  const Vect &b, iteration &iter)
386  { cg(ASM, x, b, *(ASM.A), identity_matrix(), iter); }
387 
388  template <typename ASM_type, typename Vect>
389  void AS_global_solve(using_gmres, const ASM_type &ASM, Vect &x,
390  const Vect &b, iteration &iter)
391  { gmres(ASM, x, b, identity_matrix(), 100, iter); }
392 
393  template <typename ASM_type, typename Vect>
394  void AS_global_solve(using_bicgstab, const ASM_type &ASM, Vect &x,
395  const Vect &b, iteration &iter)
396  { bicgstab(ASM, x, b, identity_matrix(), iter); }
397 
398  template <typename ASM_type, typename Vect>
399  void AS_global_solve(using_qmr,const ASM_type &ASM, Vect &x,
400  const Vect &b, iteration &iter)
401  { qmr(ASM, x, b, identity_matrix(), iter); }
402 
403 #if defined(GMM_USES_SUPERLU)
404  template <typename ASM_type, typename Vect>
405  void AS_global_solve(using_superlu, const ASM_type &, Vect &,
406  const Vect &, iteration &) {
407  GMM_ASSERT1(false, "You cannot use SuperLU as "
408  "global solver in additive Schwarz meethod");
409  }
410 #endif
411 
412  /* ******************************************************************** */
413  /* Linear Additive Schwarz method */
414  /* ******************************************************************** */
415  /* ref : Domain decomposition algorithms for the p-version finite */
416  /* element method for elliptic problems, Luca F. Pavarino, */
417  /* PhD thesis, Courant Institute of Mathematical Sciences, 1992. */
418  /* ******************************************************************** */
419 
420  /** Function to call if the ASM matrix is precomputed for successive solve
421  * with the same system.
422  */
423  template <typename Matrix1, typename Matrix2,
424  typename Vector2, typename Vector3, typename Precond,
425  typename local_solver, typename global_solver>
426  void additive_schwarz(
427  add_schwarz_mat<Matrix1, Matrix2, Precond, local_solver> &ASM, Vector3 &u,
428  const Vector2 &f, iteration &iter, const global_solver&) {
429 
430  typedef typename linalg_traits<Matrix1>::value_type value_type;
431 
432  size_type nb_sub = ASM.vB->size(), nb_dof = gmm::vect_size(f);
433  ASM.itebilan = 0;
434  std::vector<value_type> g(nb_dof);
435  std::vector<value_type> gbis(nb_dof);
436 #ifdef GMM_USES_MPI
437  double t_init=MPI_Wtime();
438  int size,tranche,borne_sup,borne_inf,rank;
439  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
440  MPI_Comm_size(MPI_COMM_WORLD, &size);
441  tranche=nb_sub/size;
442  borne_inf=rank*tranche;
443  borne_sup=(rank+1)*tranche;
444  // if (rank==size-1) borne_sup=nb_sub*size;
445  for (size_type i = size_type(borne_inf); i < size_type(borne_sup); ++i)
446 // for (size_type i = 0; i < nb_sub/size; ++i)
447  // for (size_type j = 0; j < nb_sub; ++j)
448  // for (size_type i = rank; i < nb_sub; i+=size)
449 #else
450  for (size_type i = 0; i < nb_sub; ++i)
451 #endif
452  {
453 
454 #ifdef GMM_USES_MPI
455  // size_type i=j; // (rank+size*(j-1)+nb_sub)%nb_sub;
456 #endif
457  gmm::mult(gmm::transposed((*(ASM.vB))[i]), f, ASM.fi[i]);
458  ASM.iter.init();
459  AS_local_solve(local_solver(), ASM.vAloc[i], ASM.gi[i], ASM.fi[i],
460  ASM.precond1[i], ASM.iter);
461  ASM.itebilan = std::max(ASM.itebilan, ASM.iter.get_iteration());
462 #ifdef GMM_USES_MPI
463  gmm::mult((*(ASM.vB))[i], ASM.gi[i], gbis,gbis);
464 #else
465  gmm::mult((*(ASM.vB))[i], ASM.gi[i], g, g);
466 #endif
467  }
468 #ifdef GMM_USES_MPI
469  cout<<"temps boucle init "<< MPI_Wtime()-t_init<<endl;
470  double t_ref,t_final;
471  t_ref=MPI_Wtime();
472  MPI_Allreduce(&(gbis[0]), &(g[0]),gmm::vect_size(g), MPI_DOUBLE,
473  MPI_SUM,MPI_COMM_WORLD);
474  t_final=MPI_Wtime();
475  cout<<"temps reduce init "<< t_final-t_ref<<endl;
476 #endif
477 #ifdef GMM_USES_MPI
478  t_ref=MPI_Wtime();
479  cout<<"begin global AS"<<endl;
480 #endif
481  AS_global_solve(global_solver(), ASM, u, g, iter);
482 #ifdef GMM_USES_MPI
483  t_final=MPI_Wtime();
484  cout<<"temps AS Global Solve "<< t_final-t_ref<<endl;
485 #endif
486  if (iter.get_noisy())
487  cout << "Total number of internal iterations : " << ASM.itebilan << endl;
488  }
489 
490  /** Global function. Compute the ASM matrix and call the previous function.
491  * The ASM matrix represent the preconditionned linear system.
492  */
493  template <typename Matrix1, typename Matrix2,
494  typename Vector2, typename Vector3, typename Precond,
495  typename local_solver, typename global_solver>
496  void additive_schwarz(const Matrix1 &A, Vector3 &u,
497  const Vector2 &f, const Precond &P,
498  const std::vector<Matrix2> &vB,
499  iteration &iter,
500  local_solver, global_solver) {
501  iter.set_rhsnorm(vect_norm2(f));
502  if (iter.get_rhsnorm() == 0.0) { gmm::clear(u); return; }
503  iteration iter2 = iter; iter2.reduce_noisy();
504  iter2.set_maxiter(size_type(-1));
505  add_schwarz_mat<Matrix1, Matrix2, Precond, local_solver>
506  ASM(A, vB, iter2, P, iter.get_resmax());
507  additive_schwarz(ASM, u, f, iter, global_solver());
508  }
509 
510  /* ******************************************************************** */
511  /* Sequential Non-Linear Additive Schwarz method */
512  /* ******************************************************************** */
513  /* ref : Nonlinearly Preconditionned Inexact Newton Algorithms, */
514  /* Xiao-Chuan Cai, David E. Keyes, */
515  /* SIAM J. Sci. Comp. 24: p183-200. l */
516  /* ******************************************************************** */
517 
518  template <typename Matrixt, typename MatrixBi>
519  class NewtonAS_struct {
520 
521  public :
522  typedef Matrixt tangent_matrix_type;
523  typedef MatrixBi B_matrix_type;
524  typedef typename linalg_traits<Matrixt>::value_type value_type;
525  typedef std::vector<value_type> Vector;
526 
527  virtual size_type size() = 0;
528  virtual const std::vector<MatrixBi> &get_vB() = 0;
529 
530  virtual void compute_F(Vector &f, Vector &x) = 0;
531  virtual void compute_tangent_matrix(Matrixt &M, Vector &x) = 0;
532  // compute Bi^T grad(F(X)) Bi
533  virtual void compute_sub_tangent_matrix(Matrixt &Mloc, Vector &x,
534  size_type i) = 0;
535  // compute Bi^T F(X)
536  virtual void compute_sub_F(Vector &fi, Vector &x, size_type i) = 0;
537 
538  virtual ~NewtonAS_struct() {}
539  };
540 
541  template <typename Matrixt, typename MatrixBi>
542  struct AS_exact_gradient {
543  const std::vector<MatrixBi> &vB;
544  std::vector<Matrixt> vM;
545  std::vector<Matrixt> vMloc;
546 
547  void init() {
548  for (size_type i = 0; i < vB.size(); ++i) {
549  Matrixt aux(gmm::mat_ncols(vB[i]), gmm::mat_ncols(vM[i]));
550  gmm::resize(vMloc[i], gmm::mat_ncols(vB[i]), gmm::mat_ncols(vB[i]));
551  gmm::mult(gmm::transposed(vB[i]), vM[i], aux);
552  gmm::mult(aux, vB[i], vMloc[i]);
553  }
554  }
555  AS_exact_gradient(const std::vector<MatrixBi> &vB_) : vB(vB_) {
556  vM.resize(vB.size()); vMloc.resize(vB.size());
557  for (size_type i = 0; i < vB.size(); ++i) {
558  gmm::resize(vM[i], gmm::mat_nrows(vB[i]), gmm::mat_nrows(vB[i]));
559  }
560  }
561  };
562 
563  template <typename Matrixt, typename MatrixBi,
564  typename Vector2, typename Vector3>
565  void mult(const AS_exact_gradient<Matrixt, MatrixBi> &M,
566  const Vector2 &p, Vector3 &q) {
567  gmm::clear(q);
568  typedef typename gmm::linalg_traits<Vector3>::value_type T;
569  std::vector<T> v(gmm::vect_size(p)), w, x;
570  for (size_type i = 0; i < M.vB.size(); ++i) {
571  w.resize(gmm::mat_ncols(M.vB[i]));
572  x.resize(gmm::mat_ncols(M.vB[i]));
573  gmm::mult(M.vM[i], p, v);
574  gmm::mult(gmm::transposed(M.vB[i]), v, w);
575 #if defined(GMM_USES_SUPERLU)
576  double rcond;
577  gmm::SuperLU_solve(M.vMloc[i], x, w, rcond);
578 #else
579  gmm::MUMPS_solve(M.vMloc[i], x, w);
580 #endif
581  // gmm::iteration iter(1E-10, 0, 100000);
582  //gmm::gmres(M.vMloc[i], x, w, gmm::identity_matrix(), 50, iter);
583  gmm::mult_add(M.vB[i], x, q);
584  }
585  }
586 
587  template <typename Matrixt, typename MatrixBi,
588  typename Vector2, typename Vector3>
589  void mult(const AS_exact_gradient<Matrixt, MatrixBi> &M,
590  const Vector2 &p, const Vector3 &q) {
591  mult(M, p, const_cast<Vector3 &>(q));
592  }
593 
594  template <typename Matrixt, typename MatrixBi,
595  typename Vector2, typename Vector3, typename Vector4>
596  void mult(const AS_exact_gradient<Matrixt, MatrixBi> &M,
597  const Vector2 &p, const Vector3 &p2, Vector4 &q)
598  { mult(M, p, q); add(p2, q); }
599 
600  template <typename Matrixt, typename MatrixBi,
601  typename Vector2, typename Vector3, typename Vector4>
602  void mult(const AS_exact_gradient<Matrixt, MatrixBi> &M,
603  const Vector2 &p, const Vector3 &p2, const Vector4 &q)
604  { mult(M, p, const_cast<Vector4 &>(q)); add(p2, q); }
605 
606  struct S_default_newton_line_search {
607 
608  double conv_alpha, conv_r;
609  size_t it, itmax, glob_it;
610 
611  double alpha, alpha_old, alpha_mult, first_res, alpha_max_ratio;
612  double alpha_min_ratio, alpha_min;
613  size_type count, count_pat;
614  bool max_ratio_reached;
615  double alpha_max_ratio_reached, r_max_ratio_reached;
616  size_type it_max_ratio_reached;
617 
618 
619  double converged_value() { return conv_alpha; };
620  double converged_residual() { return conv_r; };
621 
622  virtual void init_search(double r, size_t git, double = 0.0) {
623  alpha_min_ratio = 0.9;
624  alpha_min = 1e-10;
625  alpha_max_ratio = 10.0;
626  alpha_mult = 0.25;
627  itmax = size_type(-1);
628  glob_it = git; if (git <= 1) count_pat = 0;
629  conv_alpha = alpha = alpha_old = 1.;
630  conv_r = first_res = r; it = 0;
631  count = 0;
632  max_ratio_reached = false;
633  }
634  virtual double next_try() {
635  alpha_old = alpha;
636  if (alpha >= 0.4) alpha *= 0.5; else alpha *= alpha_mult; ++it;
637  return alpha_old;
638  }
639  virtual bool is_converged(double r, double = 0.0) {
640  // cout << "r = " << r << " alpha = " << alpha / alpha_mult << " count_pat = " << count_pat << endl;
641  if (!max_ratio_reached && r < first_res * alpha_max_ratio) {
642  alpha_max_ratio_reached = alpha_old; r_max_ratio_reached = r;
643  it_max_ratio_reached = it; max_ratio_reached = true;
644  }
645  if (max_ratio_reached && r < r_max_ratio_reached * 0.5
646  && r > first_res * 1.1 && it <= it_max_ratio_reached+1) {
647  alpha_max_ratio_reached = alpha_old; r_max_ratio_reached = r;
648  it_max_ratio_reached = it;
649  }
650  if (count == 0 || r < conv_r)
651  { conv_r = r; conv_alpha = alpha_old; count = 1; }
652  if (conv_r < first_res) ++count;
653 
654  if (r < first_res * alpha_min_ratio)
655  { count_pat = 0; return true; }
656  if (count >= 5 || (alpha < alpha_min && max_ratio_reached)) {
657  if (conv_r < first_res * 0.99) count_pat = 0;
658  if (/*gmm::random() * 50. < -log(conv_alpha)-4.0 ||*/ count_pat >= 3)
659  { conv_r=r_max_ratio_reached; conv_alpha=alpha_max_ratio_reached; }
660  if (conv_r >= first_res * 0.9999) count_pat++;
661  return true;
662  }
663  return false;
664  }
665  S_default_newton_line_search() { count_pat = 0; }
666  };
667 
668 
669 
670  template <typename Matrixt, typename MatrixBi, typename Vector,
671  typename Precond, typename local_solver, typename global_solver>
672  void Newton_additive_Schwarz(NewtonAS_struct<Matrixt, MatrixBi> &NS,
673  const Vector &u_,
674  iteration &iter, const Precond &P,
675  local_solver, global_solver) {
676  Vector &u = const_cast<Vector &>(u_);
677  typedef typename linalg_traits<Vector>::value_type value_type;
678  typedef typename number_traits<value_type>::magnitude_type mtype;
679  typedef actual_precond<Precond, local_solver, Matrixt> chgt_precond;
680 
681  double residual = iter.get_resmax();
682 
683  S_default_newton_line_search internal_ls;
684  S_default_newton_line_search external_ls;
685 
686  typename chgt_precond::APrecond PP = chgt_precond::transform(P);
687  iter.set_rhsnorm(mtype(1));
688  iteration iternc(iter);
689  iternc.reduce_noisy(); iternc.set_maxiter(size_type(-1));
690  iteration iter2(iternc);
691  iteration iter3(iter2); iter3.reduce_noisy();
692  iteration iter4(iter3);
693  iternc.set_name("Local Newton");
694  iter2.set_name("Linear System for Global Newton");
695  iternc.set_resmax(residual/100.0);
696  iter3.set_resmax(residual/10000.0);
697  iter2.set_resmax(residual/1000.0);
698  iter4.set_resmax(residual/1000.0);
699  std::vector<value_type> rhs(NS.size()), x(NS.size()), d(NS.size());
700  std::vector<value_type> xi, xii, fi, di;
701 
702  std::vector< std::vector<value_type> > vx(NS.get_vB().size());
703  for (size_type i = 0; i < NS.get_vB().size(); ++i) // for exact gradient
704  vx[i].resize(NS.size()); // for exact gradient
705 
706  Matrixt Mloc, M(NS.size(), NS.size());
707  NS.compute_F(rhs, u);
708  mtype act_res=gmm::vect_norm2(rhs), act_res_new(0), precond_res = act_res;
709  mtype alpha;
710 
711  while(!iter.finished(std::min(act_res, precond_res))) {
712  for (int SOR_step = 0; SOR_step >= 0; --SOR_step) {
713  gmm::clear(rhs);
714  for (size_type isd = 0; isd < NS.get_vB().size(); ++isd) {
715  const MatrixBi &Bi = (NS.get_vB())[isd];
716  size_type si = mat_ncols(Bi);
717  gmm::resize(Mloc, si, si);
718  xi.resize(si); xii.resize(si); fi.resize(si); di.resize(si);
719 
720  iternc.init();
721  iternc.set_maxiter(30); // ?
722  if (iternc.get_noisy())
723  cout << "Non-linear local problem " << isd << endl;
724  gmm::clear(xi);
725  gmm::copy(u, x);
726  NS.compute_sub_F(fi, x, isd); gmm::scale(fi, value_type(-1));
727  mtype r = gmm::vect_norm2(fi), r_t(r);
728  if (r > value_type(0)) {
729  iternc.set_rhsnorm(std::max(r, mtype(1)));
730  while(!iternc.finished(r)) {
731  NS.compute_sub_tangent_matrix(Mloc, x, isd);
732 
733  PP.build_with(Mloc);
734  iter3.init();
735  AS_local_solve(local_solver(), Mloc, di, fi, PP, iter3);
736 
737  internal_ls.init_search(r, iternc.get_iteration());
738  do {
739  alpha = internal_ls.next_try();
740  gmm::add(xi, gmm::scaled(di, -alpha), xii);
741  gmm::mult(Bi, gmm::scaled(xii, -1.0), u, x);
742  NS.compute_sub_F(fi, x, isd); gmm::scale(fi, value_type(-1));
743  r_t = gmm::vect_norm2(fi);
744  } while (!internal_ls.is_converged(r_t));
745 
746  if (alpha != internal_ls.converged_value()) {
747  alpha = internal_ls.converged_value();
748  gmm::add(xi, gmm::scaled(di, -alpha), xii);
749  gmm::mult(Bi, gmm::scaled(xii, -1.0), u, x);
750  NS.compute_sub_F(fi, x, isd); gmm::scale(fi, value_type(-1));
751  r_t = gmm::vect_norm2(fi);
752  }
753  gmm::copy(x, vx[isd]); // for exact gradient
754 
755  if (iternc.get_noisy()) cout << "(step=" << alpha << ")\t";
756  ++iternc; r = r_t; gmm::copy(xii, xi);
757  }
758  if (SOR_step) gmm::mult(Bi, gmm::scaled(xii, -1.0), u, u);
759  gmm::mult(Bi, gmm::scaled(xii, -1.0), rhs, rhs);
760  }
761  }
762  precond_res = gmm::vect_norm2(rhs);
763  if (SOR_step) cout << "SOR step residual = " << precond_res << endl;
764  if (precond_res < residual) break;
765  cout << "Precond residual = " << precond_res << endl;
766  }
767 
768  iter2.init();
769  // solving linear system for the global Newton method
770  if (0) {
771  NS.compute_tangent_matrix(M, u);
772  add_schwarz_mat<Matrixt, MatrixBi, Precond, local_solver>
773  ASM(M, NS.get_vB(), iter4, P, iter.get_resmax());
774  AS_global_solve(global_solver(), ASM, d, rhs, iter2);
775  }
776  else { // for exact gradient
777  AS_exact_gradient<Matrixt, MatrixBi> eg(NS.get_vB());
778  for (size_type i = 0; i < NS.get_vB().size(); ++i) {
779  NS.compute_tangent_matrix(eg.vM[i], vx[i]);
780  }
781  eg.init();
782  gmres(eg, d, rhs, gmm::identity_matrix(), 50, iter2);
783  }
784 
785  // gmm::add(gmm::scaled(rhs, 0.1), u); ++iter;
786  external_ls.init_search(act_res, iter.get_iteration());
787  do {
788  alpha = external_ls.next_try();
789  gmm::add(gmm::scaled(d, alpha), u, x);
790  NS.compute_F(rhs, x);
791  act_res_new = gmm::vect_norm2(rhs);
792  } while (!external_ls.is_converged(act_res_new));
793 
794  if (alpha != external_ls.converged_value()) {
795  alpha = external_ls.converged_value();
796  gmm::add(gmm::scaled(d, alpha), u, x);
797  NS.compute_F(rhs, x);
798  act_res_new = gmm::vect_norm2(rhs);
799  }
800 
801  if (iter.get_noisy() > 1) cout << endl;
802  act_res = act_res_new;
803  if (iter.get_noisy())
804  cout << "(step=" << alpha
805  << ")\t unprecond res = " << act_res << " ";
806 
807  ++iter; gmm::copy(x, u);
808  }
809  }
810 
811 }
812 
813 
814 #endif // GMM_SOLVERS_SCHWARZ_ADDITIVE_H__
gmm::iteration
The Iteration object calculates whether the solution has reached the desired accuracy,...
Definition: gmm_iter.h:53
gmm_MUMPS_interface.h
Interface with MUMPS (LU direct solver for sparse matrices).
gmm::mult_add
void mult_add(const L1 &l1, const L2 &l2, L3 &l3)
*‍/
Definition: gmm_blas.h:1791
gmm::copy
void copy(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:978
gmm::vect_norm2
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norm2(const V &v)
Euclidean norm of a vector.
Definition: gmm_blas.h:558
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:59
gmm::resize
void resize(V &v, size_type n)
*‍/
Definition: gmm_blas.h:210
gmm::mult
void mult(const L1 &l1, const L2 &l2, L3 &l3)
*‍/
Definition: gmm_blas.h:1664
gmm::add
void add(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:1277
gmm_kernel.h
Include the base gmm files.
gmm::additive_schwarz
void additive_schwarz(add_schwarz_mat< Matrix1, Matrix2, Precond, local_solver > &ASM, Vector3 &u, const Vector2 &f, iteration &iter, const global_solver &)
Function to call if the ASM matrix is precomputed for successive solve with the same system.
Definition: gmm_solver_Schwarz_additive.h:426
gmm_solver_bicgstab.h
BiCGStab iterative solver.
gmm_solver_cg.h
Conjugate gradient iterative solver.
gmm_solver_gmres.h
GMRES (Generalized Minimum Residual) iterative solver.
gmm::gmres
void gmres(const Mat &A, Vec &x, const VecB &b, const Precond &M, int restart, iteration &outer, Basis &KS)
Generalized Minimum Residual.
Definition: gmm_solver_gmres.h:90
gmm_solver_qmr.h
Quasi-Minimal Residual iterative solver.
gmm::qmr
void qmr(const Matrix &A, Vector &x, const VectorB &b, const Precond1 &M1, iteration &iter)
Quasi-Minimal Residual.
Definition: gmm_solver_qmr.h:94
gmm_superlu_interface.h
Interface with SuperLU (LU direct solver for sparse matrices).
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
bgeot::alpha
size_type alpha(short_type n, short_type d)
Return the value of which is the number of monomials of a polynomial of variables and degree .
Definition: bgeot_poly.cc:47