MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
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
37using namespace std;
38using 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]. */
52void RationalApproximation_AAA(const Vector &val, const Vector &pt,
54 real_t 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
70 Array<real_t> c_i;
71
72 // mean of the value vector
73 Vector R(val.Size());
74 real_t 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 real_t tmp_max = 0;
83 for (int j = 0; j < size; j++)
84 {
85 real_t 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<real_t> 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());
118 Add(C,Ctemp,-1.0,A);
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,'N','A');
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]. */
175void ComputePolesAndZeros(const Vector &z, const Vector &f, const Vector &w,
176 Array<real_t> & poles, Array<real_t> & zeros, real_t &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 */
246 Array<real_t> & zeros, Array<real_t> & 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 real_t tmp_numer=1.0;
263 for (int j=0; j<zsize; j++)
264 {
265 tmp_numer *= poles[i]-zeros[j];
266 }
267
268 real_t 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<real_t> & coeffs, Array<real_t> & poles,
297 real_t lmax = 1000.,
298 real_t 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 {
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 real_t eps = std::numeric_limits<real_t>::epsilon();
324
325 if (abs(alpha - 0.33) < eps)
326 {
327 coeffs = Array<real_t> ({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<real_t> ({-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<real_t>({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<real_t> ({-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)
350 {
351 alpha = 0.5;
352 }
353 coeffs = Array<real_t>({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<real_t>({-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#else
372 MFEM_CONTRACT_VAR(print_warning);
373#endif
374
375 Vector x(npoints);
376 Vector val(npoints);
377 real_t dx = lmax / (real_t)(npoints-1);
378 for (int i = 0; i<npoints; i++)
379 {
380 x(i) = dx * (real_t)i;
381 val(i) = pow(x(i),1.-alpha);
382 }
383
384 // Apply triple-A algorithm to f(x) = x^{1-a}
385 Array<real_t> z, f;
386 Vector w;
387 RationalApproximation_AAA(val,x,z,f,w,tol,max_order);
388
389 Vector vecz, vecf;
390 vecz.SetDataAndSize(z.GetData(), z.Size());
391 vecf.SetDataAndSize(f.GetData(), f.Size());
392
393 // Compute poles and zeros for RA of f(x) = x^{1-a}
394 real_t scale;
395 Array<real_t> zeros;
396 ComputePolesAndZeros(vecz, vecf, w, poles, zeros, scale);
397
398 // Remove the zero at x=0, thus, delivering a RA for f(x) = x^{-a}
399 zeros.DeleteFirst(0.0);
400
401 // Compute partial fraction approximation of f(x) = x^{-a}
402 PartialFractionExpansion(scale, poles, zeros, coeffs);
403}
void DeleteFirst(const T &el)
Delete the first entry with value == 'el'.
Definition array.hpp:847
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition array.hpp:697
int Size() const
Return the logical size of the array.
Definition array.hpp:144
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Definition array.hpp:769
T * GetData()
Returns the data.
Definition array.hpp:118
Class for Singular Value Decomposition of a DenseMatrix.
Definition densemat.hpp:956
DenseMatrix & RightSingularvectors()
Return right singular vectors.
void Eval(DenseMatrix &M)
Evaluate the SVD.
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
void Mult(const real_t *x, real_t *y) const
Matrix vector multiplication.
Definition densemat.cpp:220
void RightScaling(const Vector &s)
RightScaling: this = this * diag(s);.
Definition densemat.cpp:427
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition densemat.hpp:105
void InvLeftScaling(const Vector &s)
InvLeftScaling this = diag(1./s) * this.
Definition densemat.cpp:414
void UseExternalData(real_t *d, int h, int w)
Change the data array and the size of the DenseMatrix.
Definition densemat.hpp:80
void LeftScaling(const Vector &s)
LeftScaling this = diag(s) * this.
Definition densemat.cpp:401
void GetRow(int r, Vector &row) const
static bool Root()
Return true if the rank in MPI_COMM_WORLD is zero.
static bool IsInitialized()
Return true if MPI has been initialized.
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition operator.hpp:66
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition operator.hpp:72
Vector data type.
Definition vector.hpp:80
void SetDataAndSize(real_t *d, int s)
Set the Vector data and size.
Definition vector.hpp:175
real_t Normlinf() const
Returns the l_infinity norm of the vector.
Definition vector.cpp:850
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
real_t Sum() const
Return the sum of the vector entries.
Definition vector.cpp:1283
const real_t alpha
Definition ex15.cpp:369
void PartialFractionExpansion(real_t scale, Array< real_t > &poles, Array< real_t > &zeros, Array< real_t > &coeffs)
Definition ex33.hpp:245
void ComputePolesAndZeros(const Vector &z, const Vector &f, const Vector &w, Array< real_t > &poles, Array< real_t > &zeros, real_t &scale)
Definition ex33.hpp:175
void ComputePartialFractionApproximation(real_t &alpha, Array< real_t > &coeffs, Array< real_t > &poles, real_t lmax=1000., real_t tol=1e-10, int npoints=1000, int max_order=100)
Definition ex33.hpp:295
void RationalApproximation_AAA(const Vector &val, const Vector &pt, Array< real_t > &z, Array< real_t > &f, Vector &w, real_t tol, int max_order)
Definition ex33.hpp:52
void mfem_error(const char *msg)
Definition error.cpp:154
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
bool IsFinite(const real_t &val)
Definition vector.hpp:507
std::function< void(const Vector &, Vector &)> f_vec(bool grad_div_problem)
Definition lor_mms.hpp:68
float real_t
Definition config.hpp:43
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30
void Add(const DenseMatrix &A, const DenseMatrix &B, real_t alpha, DenseMatrix &C)
C = A + alpha*B.