58 int size = val.
Size();
59 MFEM_VERIFY(pt.
Size() == size,
"size mismatch");
63 for (
int i = 0; i < size; i++) { J[i] = i; }
76 for (
int i = 0; i<R.
Size(); i++) { R(i) = mean_val; }
78 for (
int k = 0; k < max_order; k++)
83 for (
int j = 0; j < size; j++)
85 real_t tmp = abs(val(j)-R(j));
102 for (
int j = 0; j < size; j++)
104 C_tmp[j] = 1.0/(pt(j)-pt(idx));
107 int h_C = C_tmp.
Size();
113 f_vec.SetDataAndSize(
f.GetData(),
f.Size());
122 int w_Am = A.
Width();
124 for (
int i = 0; i<h_Am; i++)
127 for (
int j = 0; j<w_Am; j++)
133#ifdef MFEM_USE_LAPACK
151 for (
int i = 0; i<J.
Size(); i++)
186 for (
int i = 1; i<=m; i++)
194#ifdef MFEM_USE_LAPACK
198 for (
int i = 0; i<evalues.
Size(); i++)
211 for (
int i = 1; i<=m; i++)
214 E(0,i) = w(i-1) *
f(i-1);
219#ifdef MFEM_USE_LAPACK
223 for (
int i = 0; i<evalues.
Size(); i++)
234 scale = w *
f / w.
Sum();
248 int psize = poles.
Size();
249 int zsize = zeros.
Size();
260 for (
int i=0; i<psize; i++)
263 for (
int j=0; j<zsize; j++)
265 tmp_numer *= poles[i]-zeros[j];
269 for (
int k=0; k<psize; k++)
271 if (k != i) { tmp_denom *= poles[i]-poles[k]; }
273 coeffs[i] *= tmp_numer / tmp_denom;
298 real_t tol=1e-10,
int npoints = 1000,
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");
307 bool print_warning =
true;
312#ifndef MFEM_USE_LAPACK
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"
323 const real_t eps = std::numeric_limits<real_t>::epsilon();
325 if (abs(
alpha - 0.33) < eps)
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});
336 else if (abs(
alpha - 0.99) < eps)
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});
349 if (abs(
alpha - 0.5) > eps)
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});
365 mfem::out <<
"=> Using precomputed values for alpha = "
366 <<
alpha <<
"\n" << std::endl;
372 MFEM_CONTRACT_VAR(print_warning);
378 for (
int i = 0; i<npoints; i++)
381 val(i) = pow(x(i),1.-
alpha);
void DeleteFirst(const T &el)
Delete the first entry with value == 'el'.
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
int Size() const
Return the logical size of the array.
int Append(const T &el)
Append element 'el' to array, resize if necessary.
T * GetData()
Returns the data.
Vector & EigenvaluesRealPart()
Class for Singular Value Decomposition of a DenseMatrix.
DenseMatrix & RightSingularvectors()
Return right singular vectors.
void Eval(DenseMatrix &M)
Evaluate the SVD.
Data type dense matrix using column-major storage.
void Mult(const real_t *x, real_t *y) const
Matrix vector multiplication.
void RightScaling(const Vector &s)
RightScaling: this = this * diag(s);.
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
void InvLeftScaling(const Vector &s)
InvLeftScaling this = diag(1./s) * this.
void UseExternalData(real_t *d, int h, int w)
Change the data array and the size of the DenseMatrix.
void LeftScaling(const Vector &s)
LeftScaling this = diag(s) * this.
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().
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
void SetDataAndSize(real_t *d, int s)
Set the Vector data and size.
real_t Normlinf() const
Returns the l_infinity norm of the vector.
int Size() const
Returns the size of the vector.
real_t Sum() const
Return the sum of the vector entries.
void PartialFractionExpansion(real_t scale, Array< real_t > &poles, Array< real_t > &zeros, Array< real_t > &coeffs)
void ComputePolesAndZeros(const Vector &z, const Vector &f, const Vector &w, Array< real_t > &poles, Array< real_t > &zeros, real_t &scale)
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)
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)
real_t f(const Vector &p)
void mfem_error(const char *msg)
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
bool IsFinite(const real_t &val)
std::function< void(const Vector &, Vector &)> f_vec(bool grad_div_problem)
void Add(const DenseMatrix &A, const DenseMatrix &B, real_t alpha, DenseMatrix &C)
C = A + alpha*B.