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