MFEM  v4.6.0 Finite element discretization library
ex33.hpp
Go to the documentation of this file.
1 // MFEM Example 33 - Serial/Parallel Shared Code
2 // (Implementation of the AAA algorithm)
3 //
4 // Here, we implement the triple-A algorithm [1] for the rational approximation
5 // of complex-valued functions,
6 //
7 // p(z)/q(z) ≈ f(z).
8 //
9 // In this file, we always assume f(z) = z^{-α}. The triple-A algorithm
10 // provides a robust, accurate approximation in rational barycentric form.
11 // This representation must be transformed into a partial fraction
12 // representation in order to be used to solve a spectral FPDE.
13 //
14 // More specifically, we first expand the numerator in terms of the zeros of
15 // the rational approximation,
16 //
17 // p(z) ∝ Π_i (z - z_i),
18 //
19 // and expand the denominator in terms of the poles of the rational
20 // approximation,
21 //
22 // q(z) ∝ Π_i (z - p_i).
23 //
24 // We then use these zeros and poles to derive the partial fraction expansion
25 //
26 // f(z) ≈ p(z)/q(z) = Σ_i c_i / (z - p_i).
27 //
28 // [1] Nakatsukasa, Y., Sète, O., & Trefethen, L. N. (2018). The AAA algorithm
29 // for rational approximation. SIAM Journal on Scientific Computing, 40(3),
30 // A1494-A1522.
31
32 #include "mfem.hpp"
33 #include <fstream>
34 #include <iostream>
35 #include <string>
36
37 using namespace std;
38 using namespace mfem;
39
40 /** RationalApproximation_AAA: compute the rational approximation (RA) of data
41  @a val [in] at the set of points @a pt [in].
42
43  @param[in] val Vector of data values
44  @param[in] pt Vector of sample points
45  @param[in] tol Relative tolerance
46  @param[in] max_order Maximum number of terms (order) of the RA
47  @param[out] z Support points of the RA in rational barycentric form
48  @param[out] f Data values at support points @a z
49  @param[out] w Weights of the RA in rational barycentric form
50
51  See pg. A1501 of Nakatsukasa et al. [1]. */
52 void RationalApproximation_AAA(const Vector &val, const Vector &pt,
54  double tol, int max_order)
55 {
56
57  // number of sample points
58  int size = val.Size();
59  MFEM_VERIFY(pt.Size() == size, "size mismatch");
60
61  // Initializations
62  Array<int> J(size);
63  for (int i = 0; i < size; i++) { J[i] = i; }
64  z.SetSize(0);
65  f.SetSize(0);
66
67  DenseMatrix C, Ctemp, A, Am;
68  // auxiliary arrays and vectors
69  Vector f_vec;
70  Array<double> c_i;
71
72  // mean of the value vector
73  Vector R(val.Size());
74  double mean_val = val.Sum()/size;
75
76  for (int i = 0; i<R.Size(); i++) { R(i) = mean_val; }
77
78  for (int k = 0; k < max_order; k++)
79  {
80  // select next support point
81  int idx = 0;
82  double tmp_max = 0;
83  for (int j = 0; j < size; j++)
84  {
85  double tmp = abs(val(j)-R(j));
86  if (tmp > tmp_max)
87  {
88  tmp_max = tmp;
89  idx = j;
90  }
91  }
92
93  // Append support points and data values
94  z.Append(pt(idx));
95  f.Append(val(idx));
96
97  // Update index vector
98  J.DeleteFirst(idx);
99
100  // next column in Cauchy matrix
101  Array<double> C_tmp(size);
102  for (int j = 0; j < size; j++)
103  {
104  C_tmp[j] = 1.0/(pt(j)-pt(idx));
105  }
106  c_i.Append(C_tmp);
107  int h_C = C_tmp.Size();
108  int w_C = k+1;
109  C.UseExternalData(c_i.GetData(),h_C,w_C);
110
111  Ctemp = C;
112
113  f_vec.SetDataAndSize(f.GetData(),f.Size());
114  Ctemp.InvLeftScaling(val);
115  Ctemp.RightScaling(f_vec);
116
117  A.SetSize(C.Height(), C.Width());
119  A.LeftScaling(val);
120
121  int h_Am = J.Size();
122  int w_Am = A.Width();
123  Am.SetSize(h_Am,w_Am);
124  for (int i = 0; i<h_Am; i++)
125  {
126  int ii = J[i];
127  for (int j = 0; j<w_Am; j++)
128  {
129  Am(i,j) = A(ii,j);
130  }
131  }
132
133 #ifdef MFEM_USE_LAPACK
134  DenseMatrixSVD svd(Am,false,true);
135  svd.Eval(Am);
137  v.GetRow(k,w);
138 #else
139  mfem_error("Compiled without LAPACK");
140 #endif
141
142  // N = C*(w.*f); D = C*w; % numerator and denominator
143  Vector aux(w);
144  aux *= f_vec;
145  Vector N(C.Height()); // Numerator
146  C.Mult(aux,N);
147  Vector D(C.Height()); // Denominator
148  C.Mult(w,D);
149
150  R = val;
151  for (int i = 0; i<J.Size(); i++)
152  {
153  int ii = J[i];
154  R(ii) = N(ii)/D(ii);
155  }
156
157  Vector verr(val);
158  verr-=R;
159
160  if (verr.Normlinf() <= tol*val.Normlinf()) { break; }
161  }
162 }
163
164 /** ComputePolesAndZeros: compute the @a poles [out] and @a zeros [out] of the
165  rational function f(z) = C p(z)/q(z) from its ration barycentric form.
166
167  @param[in] z Support points in rational barycentric form
168  @param[in] f Data values at support points @a z
169  @param[in] w Weights in rational barycentric form
170  @param[out] poles Array of poles (roots of p(z))
171  @param[out] zeros Array of zeros (roots of q(z))
172  @param[out] scale Scaling constant in f(z) = C p(z)/q(z)
173
174  See pg. A1501 of Nakatsukasa et al. [1]. */
175 void ComputePolesAndZeros(const Vector &z, const Vector &f, const Vector &w,
176  Array<double> & poles, Array<double> & zeros, double &scale)
177 {
178  // Initialization
179  poles.SetSize(0);
180  zeros.SetSize(0);
181
182  // Compute the poles
183  int m = w.Size();
184  DenseMatrix B(m+1); B = 0.;
185  DenseMatrix E(m+1); E = 0.;
186  for (int i = 1; i<=m; i++)
187  {
188  B(i,i) = 1.;
189  E(0,i) = w(i-1);
190  E(i,0) = 1.;
191  E(i,i) = z(i-1);
192  }
193
194 #ifdef MFEM_USE_LAPACK
196  eig1.Eval();
197  Vector & evalues = eig1.EigenvaluesRealPart();
198  for (int i = 0; i<evalues.Size(); i++)
199  {
200  if (IsFinite(evalues(i)))
201  {
202  poles.Append(evalues(i));
203  }
204  }
205 #else
206  mfem_error("Compiled without LAPACK");
207 #endif
208  // compute the zeros
209  B = 0.;
210  E = 0.;
211  for (int i = 1; i<=m; i++)
212  {
213  B(i,i) = 1.;
214  E(0,i) = w(i-1) * f(i-1);
215  E(i,0) = 1.;
216  E(i,i) = z(i-1);
217  }
218
219 #ifdef MFEM_USE_LAPACK
221  eig2.Eval();
222  evalues = eig2.EigenvaluesRealPart();
223  for (int i = 0; i<evalues.Size(); i++)
224  {
225  if (IsFinite(evalues(i)))
226  {
227  zeros.Append(evalues(i));
228  }
229  }
230 #else
231  mfem_error("Compiled without LAPACK");
232 #endif
233
234  scale = w * f / w.Sum();
235 }
236
237 /** PartialFractionExpansion: compute the partial fraction expansion of the
238  rational function f(z) = Σ_i c_i / (z - p_i) from its @a poles [in] and
239  @a zeros [in].
240
241  @param[in] poles Array of poles (same as p_i above)
242  @param[in] zeros Array of zeros
243  @param[in] scale Scaling constant
244  @param[out] coeffs Coefficients c_i */
245 void PartialFractionExpansion(double scale, Array<double> & poles,
246  Array<double> & zeros, Array<double> & coeffs)
247 {
248  int psize = poles.Size();
249  int zsize = zeros.Size();
250  coeffs.SetSize(psize);
251  coeffs = scale;
252
253  // Note: C p(z)/q(z) = Σ_i c_i / (z - p_i) results in an system of equations
254  // where the N unknowns are the coefficients c_i. After multiplying the
255  // system with q(z), the coefficients c_i can be computed analytically by
256  // choosing N values for z. Choosing z_j = = p_j diagonalizes the system and
257  // one can obtain an analytic form for the c_i coefficients. The result is
258  // implemented in the code block below.
259
260  for (int i=0; i<psize; i++)
261  {
262  double tmp_numer=1.0;
263  for (int j=0; j<zsize; j++)
264  {
265  tmp_numer *= poles[i]-zeros[j];
266  }
267
268  double tmp_denom=1.0;
269  for (int k=0; k<psize; k++)
270  {
271  if (k != i) { tmp_denom *= poles[i]-poles[k]; }
272  }
273  coeffs[i] *= tmp_numer / tmp_denom;
274  }
275 }
276
277
278 /** ComputePartialFractionApproximation: compute a rational approximation (RA)
279  in partial fraction form, e.g., f(z) ≈ Σ_i c_i / (z - p_i), from sampled
280  values of the function f(z) = z^{-a}, 0 < a < 1.
281
282  @param[in] alpha Exponent a in f(z) = z^-a
283  @param[in] lmax, npoints f(z) is uniformly sampled @a npoints times in the
284  interval [ 0, @a lmax ]
285  @param[in] tol Relative tolerance
286  @param[in] max_order Maximum number of terms (order) of the RA
287  @param[out] coeffs Coefficients c_i
288  @param[out] poles Poles p_i
289
290  NOTES: When MFEM is not built with LAPACK support, only @a alpha = 0.33,
291  0.5, and 0.99 are possible. In this case, if @a alpha != 0.33 and
292  @a alpha != 0.99, then @a alpha = 0.5 is used by default.
293
294  See pg. A1501 of Nakatsukasa et al. [1]. */
296  Array<double> & coeffs, Array<double> & poles,
297  double lmax = 1000.,
298  double tol=1e-10, int npoints = 1000,
299  int max_order = 100)
300 {
301  MFEM_VERIFY(alpha < 1., "alpha must be less than 1");
302  MFEM_VERIFY(alpha > 0., "alpha must be greater than 0");
303  MFEM_VERIFY(npoints > 2, "npoints must be greater than 2");
304  MFEM_VERIFY(lmax > 0, "lmin must be greater than 0");
305  MFEM_VERIFY(tol > 0, "tol must be greater than 0");
306
307  bool print_warning = true;
308 #ifdef MFEM_USE_MPI
309  if ((Mpi::IsInitialized() && !Mpi::Root())) { print_warning = false; }
310 #endif
311
312 #ifndef MFEM_USE_LAPACK
313  if (print_warning)
314  {
315  mfem::out
316  << "\n" << string(80, '=')
317  << "\nMFEM is compiled without LAPACK."
318  << "\nUsing precomputed values for PartialFractionApproximation."
319  << "\nOnly alpha = 0.33, 0.5, and 0.99 are available."
320  << "\nThe default is alpha = 0.5.\n" << string(80, '=') << "\n"
321  << endl;
322  }
323  const double eps = std::numeric_limits<double>::epsilon();
324
325  if (abs(alpha - 0.33) < eps)
326  {
327  coeffs = Array<double> ({1.821898e+03, 9.101221e+01, 2.650611e+01,
328  1.174937e+01, 6.140444e+00, 3.441713e+00,
329  1.985735e+00, 1.162634e+00, 6.891560e-01,
330  4.111574e-01, 2.298736e-01});
331  poles = Array<double> ({-4.155583e+04, -2.956285e+03, -8.331715e+02,
332  -3.139332e+02, -1.303448e+02, -5.563385e+01,
333  -2.356255e+01, -9.595516e+00, -3.552160e+00,
334  -1.032136e+00, -1.241480e-01});
335  }
336  else if (abs(alpha - 0.99) < eps)
337  {
338  coeffs = Array<double>({2.919591e-02, 1.419750e-02, 1.065798e-02,
339  9.395094e-03, 8.915329e-03, 8.822991e-03,
340  9.058247e-03, 9.814521e-03, 1.180396e-02,
341  1.834554e-02, 9.840482e-01});
342  poles = Array<double> ({-1.069683e+04, -1.769370e+03, -5.718374e+02,
343  -2.242095e+02, -9.419132e+01, -4.031012e+01,
344  -1.701525e+01, -6.810088e+00, -2.382810e+00,
345  -5.700059e-01, -1.384324e-03});
346  }
347  else
348  {
349  if (abs(alpha - 0.5) > eps && print_warning)
350  {
351  alpha = 0.5;
352  }
353  coeffs = Array<double>({2.290262e+02, 2.641819e+01, 1.005566e+01,
354  5.390411e+00, 3.340725e+00, 2.211205e+00,
355  1.508883e+00, 1.049474e+00, 7.462709e-01,
356  5.482686e-01, 4.232510e-01, 3.578967e-01});
357  poles = Array<double>({-3.168211e+04, -3.236077e+03, -9.868287e+02,
358  -3.945597e+02, -1.738889e+02, -7.925178e+01,
359  -3.624992e+01, -1.629196e+01, -6.982956e+00,
360  -2.679984e+00, -7.782607e-01, -7.649166e-02});
361  }
362
363  if (print_warning)
364  {
365  mfem::out << "=> Using precomputed values for alpha = "
366  << alpha << "\n" << std::endl;
367  }
368
369
370  return;
371 #endif
372
373  Vector x(npoints);
374  Vector val(npoints);
375  double dx = lmax / (double)(npoints-1);
376  for (int i = 0; i<npoints; i++)
377  {
378  x(i) = dx * (double)i;
379  val(i) = pow(x(i),1.-alpha);
380  }
381
382  // Apply triple-A algorithm to f(x) = x^{1-a}
383  Array<double> z, f;
384  Vector w;
385  RationalApproximation_AAA(val,x,z,f,w,tol,max_order);
386
387  Vector vecz, vecf;
388  vecz.SetDataAndSize(z.GetData(), z.Size());
389  vecf.SetDataAndSize(f.GetData(), f.Size());
390
391  // Compute poles and zeros for RA of f(x) = x^{1-a}
392  double scale;
393  Array<double> zeros;
394  ComputePolesAndZeros(vecz, vecf, w, poles, zeros, scale);
395
396  // Remove the zero at x=0, thus, delivering a RA for f(x) = x^{-a}
397  zeros.DeleteFirst(0.0);
398
399  // Compute partial fraction approximation of f(x) = x^{-a}
400  PartialFractionExpansion(scale, poles, zeros, coeffs);
401 }
void RationalApproximation_AAA(const Vector &val, const Vector &pt, Array< double > &z, Array< double > &f, Vector &w, double tol, int max_order)
Definition: ex33.hpp:52
double epsilon
Definition: ex25.cpp:141
void PartialFractionExpansion(double scale, Array< double > &poles, Array< double > &zeros, Array< double > &coeffs)
Definition: ex33.hpp:245
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
T * GetData()
Returns the data.
Definition: array.hpp:115
int Size() const
Returns the size of the vector.
Definition: vector.hpp:197
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
void Eval(DenseMatrix &M)
Evaluate the SVD.
Definition: densemat.cpp:4190
STL namespace.
void DeleteFirst(const T &el)
Delete the first entry with value == &#39;el&#39;.
Definition: array.hpp:837
bool IsFinite(const double &val)
Definition: vector.hpp:486
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:2320
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:759
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:154
double Sum() const
Return the sum of the vector entries.
Definition: vector.cpp:1286
void LeftScaling(const Vector &s)
LeftScaling this = diag(s) * this.
Definition: densemat.cpp:353
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:172
std::function< void(const Vector &, Vector &)> f_vec(bool grad_div_problem)
Definition: lor_mms.hpp:68
void GetRow(int r, Vector &row) const
Definition: densemat.cpp:1314
Class for Singular Value Decomposition of a DenseMatrix.
Definition: densemat.hpp:948
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:66
void ComputePartialFractionApproximation(double &alpha, Array< double > &coeffs, Array< double > &poles, double lmax=1000., double tol=1e-10, int npoints=1000, int max_order=100)
Definition: ex33.hpp:295
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:687
void ComputePolesAndZeros(const Vector &z, const Vector &f, const Vector &w, Array< double > &poles, Array< double > &zeros, double &scale)
Definition: ex33.hpp:175
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
void SetDataAndSize(double *d, int s)
Set the Vector data and size.
Definition: vector.hpp:154
DenseMatrix & RightSingularvectors()
Return right singular vectors.
Definition: densemat.hpp:1086
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
const double alpha
Definition: ex15.cpp:369
void RightScaling(const Vector &s)
RightScaling: this = this * diag(s);.
Definition: densemat.cpp:379
Vector data type.
Definition: vector.hpp:58
void InvLeftScaling(const Vector &s)
InvLeftScaling this = diag(1./s) * this.
Definition: densemat.cpp:366
void UseExternalData(double *d, int h, int w)
Change the data array and the size of the DenseMatrix.
Definition: densemat.hpp:80
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
double Normlinf() const
Returns the l_infinity norm of the vector.
Definition: vector.cpp:853
double f(const Vector &p)