dune-istl  2.8.0
gsetc.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_ISTL_GSETC_HH
4 #define DUNE_ISTL_GSETC_HH
5 
6 #include <cmath>
7 #include <complex>
8 #include <iostream>
9 #include <iomanip>
10 #include <string>
11 
12 #include <dune/common/hybridutilities.hh>
13 
14 #include "multitypeblockvector.hh"
15 #include "multitypeblockmatrix.hh"
16 
17 #include "istlexception.hh"
18 
19 
25 namespace Dune {
26 
37  //============================================================
38  // parameter types
39  //============================================================
40 
42  template<int l>
43  struct BL {
44  enum {recursion_level = l};
45  };
46 
47  enum WithDiagType {
49  nodiag=0
50  };
51 
54  norelax=0
55  };
56 
57  //============================================================
58  // generic triangular solves
59  // consider block decomposition A = L + D + U
60  // we can invert L, L+D, U, U+D
61  // we can apply relaxation or not
62  // we can recurse over a fixed number of levels
63  //============================================================
64 
65  // template meta program for triangular solves
66  template<int I, WithDiagType diag, WithRelaxType relax>
67  struct algmeta_btsolve {
68  template<class M, class X, class Y, class K>
69  static void bltsolve (const M& A, X& v, const Y& d, const K& w)
70  {
71  // iterator types
72  typedef typename M::ConstRowIterator rowiterator;
73  typedef typename M::ConstColIterator coliterator;
74  typedef typename Y::block_type bblock;
75 
76  // local solve at each block and immediate update
77  rowiterator endi=A.end();
78  for (rowiterator i=A.begin(); i!=endi; ++i)
79  {
80  bblock rhs(d[i.index()]);
81  coliterator j;
82  for (j=(*i).begin(); j.index()<i.index(); ++j)
83  (*j).mmv(v[j.index()],rhs);
84  algmeta_btsolve<I-1,diag,relax>::bltsolve(*j,v[i.index()],rhs,w);
85  }
86  }
87  template<class M, class X, class Y, class K>
88  static void butsolve (const M& A, X& v, const Y& d, const K& w)
89  {
90  // iterator types
91  typedef typename M::ConstRowIterator rowiterator;
92  typedef typename M::ConstColIterator coliterator;
93  typedef typename Y::block_type bblock;
94 
95  // local solve at each block and immediate update
96  rowiterator rendi=A.beforeBegin();
97  for (rowiterator i=A.beforeEnd(); i!=rendi; --i)
98  {
99  bblock rhs(d[i.index()]);
100  coliterator j;
101  for (j=(*i).beforeEnd(); j.index()>i.index(); --j)
102  (*j).mmv(v[j.index()],rhs);
103  algmeta_btsolve<I-1,diag,relax>::butsolve(*j,v[i.index()],rhs,w);
104  }
105  }
106  };
107 
108  // recursion end ...
109  template<>
111  template<class M, class X, class Y, class K>
112  static void bltsolve (const M& A, X& v, const Y& d, const K& w)
113  {
114  A.solve(v,d);
115  v *= w;
116  }
117  template<class M, class X, class Y, class K>
118  static void butsolve (const M& A, X& v, const Y& d, const K& w)
119  {
120  A.solve(v,d);
121  v *= w;
122  }
123  };
124  template<>
126  template<class M, class X, class Y, class K>
127  static void bltsolve (const M& A, X& v, const Y& d, const K& /*w*/)
128  {
129  A.solve(v,d);
130  }
131  template<class M, class X, class Y, class K>
132  static void butsolve (const M& A, X& v, const Y& d, const K& /*w*/)
133  {
134  A.solve(v,d);
135  }
136  };
137  template<>
139  template<class M, class X, class Y, class K>
140  static void bltsolve (const M& /*A*/, X& v, const Y& d, const K& w)
141  {
142  v = d;
143  v *= w;
144  }
145  template<class M, class X, class Y, class K>
146  static void butsolve (const M& /*A*/, X& v, const Y& d, const K& w)
147  {
148  v = d;
149  v *= w;
150  }
151  };
152  template<>
154  template<class M, class X, class Y, class K>
155  static void bltsolve (const M& /*A*/, X& v, const Y& d, const K& /*w*/)
156  {
157  v = d;
158  }
159  template<class M, class X, class Y, class K>
160  static void butsolve (const M& /*A*/, X& v, const Y& d, const K& /*w*/)
161  {
162  v = d;
163  }
164  };
165 
166 
167  // user calls
168 
169  // default block recursion level = 1
170 
172  template<class M, class X, class Y>
173  void bltsolve (const M& A, X& v, const Y& d)
174  {
175  typename X::field_type w=1;
177  }
179  template<class M, class X, class Y, class K>
180  void bltsolve (const M& A, X& v, const Y& d, const K& w)
181  {
183  }
185  template<class M, class X, class Y>
186  void ubltsolve (const M& A, X& v, const Y& d)
187  {
188  typename X::field_type w=1;
190  }
192  template<class M, class X, class Y, class K>
193  void ubltsolve (const M& A, X& v, const Y& d, const K& w)
194  {
196  }
197 
199  template<class M, class X, class Y>
200  void butsolve (const M& A, X& v, const Y& d)
201  {
202  typename X::field_type w=1;
204  }
206  template<class M, class X, class Y, class K>
207  void butsolve (const M& A, X& v, const Y& d, const K& w)
208  {
210  }
212  template<class M, class X, class Y>
213  void ubutsolve (const M& A, X& v, const Y& d)
214  {
215  typename X::field_type w=1;
217  }
219  template<class M, class X, class Y, class K>
220  void ubutsolve (const M& A, X& v, const Y& d, const K& w)
221  {
223  }
224 
225  // general block recursion level >= 0
226 
228  template<class M, class X, class Y, int l>
229  void bltsolve (const M& A, X& v, const Y& d, BL<l> /*bl*/)
230  {
231  typename X::field_type w=1;
233  }
235  template<class M, class X, class Y, class K, int l>
236  void bltsolve (const M& A, X& v, const Y& d, const K& w, BL<l> /*bl*/)
237  {
239  }
241  template<class M, class X, class Y, int l>
242  void ubltsolve (const M& A, X& v, const Y& d, BL<l> /*bl*/)
243  {
244  typename X::field_type w=1;
246  }
248  template<class M, class X, class Y, class K, int l>
249  void ubltsolve (const M& A, X& v, const Y& d, const K& w, BL<l> /*bl*/)
250  {
252  }
253 
255  template<class M, class X, class Y, int l>
256  void butsolve (const M& A, X& v, const Y& d, BL<l> bl)
257  {
258  typename X::field_type w=1;
260  }
262  template<class M, class X, class Y, class K, int l>
263  void butsolve (const M& A, X& v, const Y& d, const K& w, BL<l> bl)
264  {
266  }
268  template<class M, class X, class Y, int l>
269  void ubutsolve (const M& A, X& v, const Y& d, BL<l> bl)
270  {
271  typename X::field_type w=1;
273  }
275  template<class M, class X, class Y, class K, int l>
276  void ubutsolve (const M& A, X& v, const Y& d, const K& w, BL<l> bl)
277  {
279  }
280 
281 
282 
283  //============================================================
284  // generic block diagonal solves
285  // consider block decomposition A = L + D + U
286  // we can apply relaxation or not
287  // we can recurse over a fixed number of levels
288  //============================================================
289 
290  // template meta program for diagonal solves
291  template<int I, WithRelaxType relax>
293  template<class M, class X, class Y, class K>
294  static void bdsolve (const M& A, X& v, const Y& d, const K& w)
295  {
296  // iterator types
297  typedef typename M::ConstRowIterator rowiterator;
298  typedef typename M::ConstColIterator coliterator;
299 
300  // local solve at each block and immediate update
301  rowiterator rendi=A.beforeBegin();
302  for (rowiterator i=A.beforeEnd(); i!=rendi; --i)
303  {
304  coliterator ii=(*i).find(i.index());
305  algmeta_bdsolve<I-1,relax>::bdsolve(*ii,v[i.index()],d[i.index()],w);
306  }
307  }
308  };
309 
310  // recursion end ...
311  template<>
313  template<class M, class X, class Y, class K>
314  static void bdsolve (const M& A, X& v, const Y& d, const K& w)
315  {
316  A.solve(v,d);
317  v *= w;
318  }
319  };
320  template<>
322  template<class M, class X, class Y, class K>
323  static void bdsolve (const M& A, X& v, const Y& d, const K& /*w*/)
324  {
325  A.solve(v,d);
326  }
327  };
328 
329  // user calls
330 
331  // default block recursion level = 1
332 
334  template<class M, class X, class Y>
335  void bdsolve (const M& A, X& v, const Y& d)
336  {
337  typename X::field_type w=1;
339  }
341  template<class M, class X, class Y, class K>
342  void bdsolve (const M& A, X& v, const Y& d, const K& w)
343  {
345  }
346 
347  // general block recursion level >= 0
348 
350  template<class M, class X, class Y, int l>
351  void bdsolve (const M& A, X& v, const Y& d, BL<l> /*bl*/)
352  {
353  typename X::field_type w=1;
355  }
357  template<class M, class X, class Y, class K, int l>
358  void bdsolve (const M& A, X& v, const Y& d, const K& w, BL<l> /*bl*/)
359  {
361  }
362 
363 
364  //============================================================
365  // generic steps of iteration methods
366  // Jacobi, Gauss-Seidel, SOR, SSOR
367  // work directly on Ax=b, ie solve M(x^{i+1}-x^i) = w (b-Ax^i)
368  // we can recurse over a fixed number of levels
369  //============================================================
370 
371  // template meta program for iterative solver steps
372  template<int I, typename M>
374 
375  template<class X, class Y, class K>
376  static void dbgs (const M& A, X& x, const Y& b, const K& w)
377  {
378  typedef typename M::ConstRowIterator rowiterator;
379  typedef typename M::ConstColIterator coliterator;
380  typedef typename Y::block_type bblock;
381  bblock rhs;
382 
383  X xold(x); // remember old x
384 
385  rowiterator endi=A.end();
386  for (rowiterator i=A.begin(); i!=endi; ++i)
387  {
388  rhs = b[i.index()]; // rhs = b_i
389  coliterator endj=(*i).end();
390  coliterator j=(*i).begin();
391  if constexpr (IsNumber<typename M::block_type>())
392  {
393  for (; j.index()<i.index(); ++j)
394  rhs -= (*j) * x[j.index()];
395  coliterator diag=j++; // *diag = a_ii and increment coliterator j from a_ii to a_i+1,i to skip diagonal
396  for (; j != endj; ++j)
397  rhs -= (*j) * x[j.index()];
398  x[i.index()] = rhs / (*diag);
399  }
400  else
401  {
402  for (; j.index()<i.index(); ++j) // iterate over a_ij with j < i
403  (*j).mmv(x[j.index()],rhs); // rhs -= sum_{j<i} a_ij * xnew_j
404  coliterator diag=j++; // *diag = a_ii and increment coliterator j from a_ii to a_i+1,i to skip diagonal
405  for (; j != endj; ++j) // iterate over a_ij with j > i
406  (*j).mmv(x[j.index()],rhs); // rhs -= sum_{j>i} a_ij * xold_j
407  algmeta_itsteps<I-1,typename M::block_type>::dbgs(*diag,x[i.index()],rhs,w); // if I==1: xnew_i = rhs/a_ii
408  }
409  }
410  // next two lines: xnew_i = w / a_ii * (b_i - sum_{j<i} a_ij * xnew_j - sum_{j>=i} a_ij * xold_j) + (1-w)*xold;
411  x *= w;
412  x.axpy(K(1)-w,xold);
413  }
414 
415  template<class X, class Y, class K>
416  static void bsorf (const M& A, X& x, const Y& b, const K& w)
417  {
418  typedef typename M::ConstRowIterator rowiterator;
419  typedef typename M::ConstColIterator coliterator;
420  typedef typename Y::block_type bblock;
421  typedef typename X::block_type xblock;
422  bblock rhs;
423  xblock v;
424 
425  // Initialize nested data structure if there are entries
426  if(A.begin()!=A.end())
427  v=x[0];
428 
429  rowiterator endi=A.end();
430  for (rowiterator i=A.begin(); i!=endi; ++i)
431  {
432  rhs = b[i.index()]; // rhs = b_i
433  coliterator endj=(*i).end(); // iterate over a_ij with j < i
434  coliterator j=(*i).begin();
435  if constexpr (IsNumber<typename M::block_type>())
436  {
437  for (; j.index()<i.index(); ++j)
438  rhs -= (*j) * x[j.index()]; // rhs -= sum_{j<i} a_ij * xnew_j
439  coliterator diag=j; // *diag = a_ii
440  for (; j!=endj; ++j)
441  rhs -= (*j) * x[j.index()]; // rhs -= sum_{j<i} a_ij * xnew_j
442  v = rhs / (*diag);
443  x[i.index()] += w*v; // x_i = w / a_ii * (b_i - sum_{j<i} a_ij * xnew_j - sum_{j>=i} a_ij * xold_j)
444  }
445  else
446  {
447  for (; j.index()<i.index(); ++j)
448  (*j).mmv(x[j.index()],rhs); // rhs -= sum_{j<i} a_ij * xnew_j
449  coliterator diag=j; // *diag = a_ii
450  for (; j!=endj; ++j)
451  (*j).mmv(x[j.index()],rhs); // rhs -= sum_{j<i} a_ij * xnew_j
452  algmeta_itsteps<I-1,typename M::block_type>::bsorf(*diag,v,rhs,w); // if blocksize I==1: v = rhs/a_ii
453  x[i.index()].axpy(w,v); // x_i = w / a_ii * (b_i - sum_{j<i} a_ij * xnew_j - sum_{j>=i} a_ij * xold_j)
454  }
455  }
456  }
457 
458  template<class X, class Y, class K>
459  static void bsorb (const M& A, X& x, const Y& b, const K& w)
460  {
461  typedef typename M::ConstRowIterator rowiterator;
462  typedef typename M::ConstColIterator coliterator;
463  typedef typename Y::block_type bblock;
464  typedef typename X::block_type xblock;
465  bblock rhs;
466  xblock v;
467 
468  // Initialize nested data structure if there are entries
469  if(A.begin()!=A.end())
470  v=x[0];
471 
472  rowiterator endi=A.beforeBegin();
473  for (rowiterator i=A.beforeEnd(); i!=endi; --i)
474  {
475  rhs = b[i.index()];
476  coliterator endj=(*i).end();
477  coliterator j=(*i).begin();
478  if constexpr (IsNumber<typename M::block_type>())
479  {
480  for (; j.index()<i.index(); ++j)
481  rhs -= (*j) * x[j.index()];
482  coliterator diag=j;
483  for (; j!=endj; ++j)
484  rhs -= (*j) * x[j.index()];
485  v = rhs / (*diag);
486  x[i.index()] += w*v;
487  }
488  else
489  {
490  for (; j.index()<i.index(); ++j)
491  j->mmv(x[j.index()],rhs);
492  coliterator diag=j;
493  for (; j!=endj; ++j)
494  j->mmv(x[j.index()],rhs);
496  x[i.index()].axpy(w,v);
497  }
498  }
499  }
500 
501  template<class X, class Y, class K>
502  static void dbjac (const M& A, X& x, const Y& b, const K& w)
503  {
504  typedef typename M::ConstRowIterator rowiterator;
505  typedef typename M::ConstColIterator coliterator;
506  typedef typename Y::block_type bblock;
507  bblock rhs;
508 
509  X v(x); // allocate with same size
510 
511  rowiterator endi=A.end();
512  for (rowiterator i=A.begin(); i!=endi; ++i)
513  {
514  rhs = b[i.index()];
515  coliterator endj=(*i).end();
516  coliterator j=(*i).begin();
517  if constexpr (IsNumber<typename M::block_type>())
518  {
519  for (; j.index()<i.index(); ++j)
520  rhs -= (*j) * x[j.index()];
521  coliterator diag=j;
522  for (; j!=endj; ++j)
523  rhs -= (*j) * x[j.index()];
524  v[i.index()] = rhs / (*diag);
525  }
526  else
527  {
528  for (; j.index()<i.index(); ++j)
529  j->mmv(x[j.index()],rhs);
530  coliterator diag=j;
531  for (; j!=endj; ++j)
532  j->mmv(x[j.index()],rhs);
533  algmeta_itsteps<I-1,typename M::block_type>::dbjac(*diag,v[i.index()],rhs,w);
534  }
535  }
536  x.axpy(w,v);
537  }
538  };
539  // end of recursion
540  template<typename M>
541  struct algmeta_itsteps<0,M> {
542  template<class X, class Y, class K>
543  static void dbgs (const M& A, X& x, const Y& b, const K& /*w*/)
544  {
545  A.solve(x,b);
546  }
547  template<class X, class Y, class K>
548  static void bsorf (const M& A, X& x, const Y& b, const K& /*w*/)
549  {
550  A.solve(x,b);
551  }
552  template<class X, class Y, class K>
553  static void bsorb (const M& A, X& x, const Y& b, const K& /*w*/)
554  {
555  A.solve(x,b);
556  }
557  template<class X, class Y, class K>
558  static void dbjac (const M& A, X& x, const Y& b, const K& /*w*/)
559  {
560  A.solve(x,b);
561  }
562  };
563 
564  template<int I, typename T1, typename... MultiTypeMatrixArgs>
565  struct algmeta_itsteps<I,MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>> {
566  template<
567  typename... MultiTypeVectorArgs,
568  class K>
572  const K& w)
573  {
576  }
577 
578  template<
579  typename... MultiTypeVectorArgs,
580  class K>
584  const K& w)
585  {
588  }
589 
590  template<
591  typename... MultiTypeVectorArgs,
592  class K>
596  const K& w)
597  {
600  }
601 
602  template<
603  typename... MultiTypeVectorArgs,
604  class K
605  >
609  const K& w)
610  {
613  }
614  };
615 
616  // user calls
617 
619  template<class M, class X, class Y, class K>
620  void dbgs (const M& A, X& x, const Y& b, const K& w)
621  {
623  }
625  template<class M, class X, class Y, class K, int l>
626  void dbgs (const M& A, X& x, const Y& b, const K& w, BL<l> /*bl*/)
627  {
629  }
631  template<class M, class X, class Y, class K>
632  void bsorf (const M& A, X& x, const Y& b, const K& w)
633  {
635  }
637  template<class M, class X, class Y, class K, int l>
638  void bsorf (const M& A, X& x, const Y& b, const K& w, BL<l> /*bl*/)
639  {
641  }
643  template<class M, class X, class Y, class K>
644  void bsorb (const M& A, X& x, const Y& b, const K& w)
645  {
647  }
649  template<class M, class X, class Y, class K, int l>
650  void bsorb (const M& A, X& x, const Y& b, const K& w, BL<l> /*bl*/)
651  {
653  }
655  template<class M, class X, class Y, class K>
656  void dbjac (const M& A, X& x, const Y& b, const K& w)
657  {
659  }
661  template<class M, class X, class Y, class K, int l>
662  void dbjac (const M& A, X& x, const Y& b, const K& w, BL<l> /*bl*/)
663  {
665  }
666 
667 
670 } // end namespace
671 
672 #endif
static void dbjac(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:566
static void dbgs(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:482
static void bsorb(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:538
static constexpr size_type N()
Return the number of matrix rows.
Definition: multitypeblockmatrix.hh:62
static void bsorf(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:511
void bltsolve(const M &A, X &v, const Y &d)
block lower triangular solve
Definition: gsetc.hh:173
WithDiagType
Definition: gsetc.hh:47
void bsorb(const M &A, X &x, const Y &b, const K &w)
SSOR step.
Definition: gsetc.hh:644
void ubltsolve(const M &A, X &v, const Y &d)
unit block lower triangular solve
Definition: gsetc.hh:186
void dbjac(const M &A, X &x, const Y &b, const K &w)
Jacobi step.
Definition: gsetc.hh:656
void dbgs(const M &A, X &x, const Y &b, const K &w)
GS step.
Definition: gsetc.hh:620
WithRelaxType
Definition: gsetc.hh:52
void bdsolve(const M &A, X &v, const Y &d)
block diagonal solve, no relaxation
Definition: gsetc.hh:335
void butsolve(const M &A, X &v, const Y &d)
block upper triangular solve
Definition: gsetc.hh:200
void bsorf(const M &A, X &x, const Y &b, const K &w)
SOR step.
Definition: gsetc.hh:632
void ubutsolve(const M &A, X &v, const Y &d)
unit block upper triangular solve
Definition: gsetc.hh:213
@ nodiag
Definition: gsetc.hh:49
@ withdiag
Definition: gsetc.hh:48
@ norelax
Definition: gsetc.hh:54
@ withrelax
Definition: gsetc.hh:53
Definition: allocator.hh:9
A Vector class to support different block types.
Definition: multitypeblockvector.hh:56
A Matrix class to support different block types.
Definition: multitypeblockmatrix.hh:44
compile-time parameter for block recursion depth
Definition: gsetc.hh:43
@ recursion_level
Definition: gsetc.hh:44
Definition: gsetc.hh:67
static void butsolve(const M &A, X &v, const Y &d, const K &w)
Definition: gsetc.hh:88
static void bltsolve(const M &A, X &v, const Y &d, const K &w)
Definition: gsetc.hh:69
static void bltsolve(const M &A, X &v, const Y &d, const K &w)
Definition: gsetc.hh:112
static void butsolve(const M &A, X &v, const Y &d, const K &w)
Definition: gsetc.hh:118
static void bltsolve(const M &A, X &v, const Y &d, const K &)
Definition: gsetc.hh:127
static void butsolve(const M &A, X &v, const Y &d, const K &)
Definition: gsetc.hh:132
static void bltsolve(const M &, X &v, const Y &d, const K &w)
Definition: gsetc.hh:140
static void butsolve(const M &, X &v, const Y &d, const K &w)
Definition: gsetc.hh:146
static void bltsolve(const M &, X &v, const Y &d, const K &)
Definition: gsetc.hh:155
static void butsolve(const M &, X &v, const Y &d, const K &)
Definition: gsetc.hh:160
Definition: gsetc.hh:292
static void bdsolve(const M &A, X &v, const Y &d, const K &w)
Definition: gsetc.hh:294
static void bdsolve(const M &A, X &v, const Y &d, const K &w)
Definition: gsetc.hh:314
static void bdsolve(const M &A, X &v, const Y &d, const K &)
Definition: gsetc.hh:323
Definition: gsetc.hh:373
static void bsorb(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:459
static void bsorf(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:416
static void dbjac(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:502
static void dbgs(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:376
static void dbgs(const M &A, X &x, const Y &b, const K &)
Definition: gsetc.hh:543
static void dbjac(const M &A, X &x, const Y &b, const K &)
Definition: gsetc.hh:558
static void bsorf(const M &A, X &x, const Y &b, const K &)
Definition: gsetc.hh:548
static void bsorb(const M &A, X &x, const Y &b, const K &)
Definition: gsetc.hh:553
static void dbgs(const MultiTypeBlockMatrix< T1, MultiTypeMatrixArgs... > &A, MultiTypeBlockVector< MultiTypeVectorArgs... > &x, const MultiTypeBlockVector< MultiTypeVectorArgs... > &b, const K &w)
Definition: gsetc.hh:569
static void dbjac(const MultiTypeBlockMatrix< T1, MultiTypeMatrixArgs... > &A, MultiTypeBlockVector< MultiTypeVectorArgs... > &x, const MultiTypeBlockVector< MultiTypeVectorArgs... > &b, const K &w)
Definition: gsetc.hh:606
static void bsorf(const MultiTypeBlockMatrix< T1, MultiTypeMatrixArgs... > &A, MultiTypeBlockVector< MultiTypeVectorArgs... > &x, const MultiTypeBlockVector< MultiTypeVectorArgs... > &b, const K &w)
Definition: gsetc.hh:581
static void bsorb(const MultiTypeBlockMatrix< T1, MultiTypeMatrixArgs... > &A, MultiTypeBlockVector< MultiTypeVectorArgs... > &x, const MultiTypeBlockVector< MultiTypeVectorArgs... > &b, const K &w)
Definition: gsetc.hh:593