SF_CHEBYCHEV

Description

Approximation of \(f(A)\,v\), the action of a function of a linear operator \(A\) on a seed vector \(v\), by a truncated Chebyshev series. The operator is never needed explicitly: only its action on a vector is used, through a user-supplied matrix-vector product (MatVec).

The operator is first mapped onto \([-1, 1]\),

\[\tilde{A} = \frac{A - b}{a},\]

where \(a > 0\) is the scale and \(b\) the shift: the spectrum of \(A\) must lie in \([b - a,\, b + a]\). The function is then expanded as

\[f(A) \simeq \sum_{m=0}^{M} c_m\, T_m(\tilde{A}),\]

where \(T_m\) is the Chebyshev polynomial of the first kind and degree \(m\). Applied to the seed vector, the polynomials are generated by the recurrence

\[T_0 v = v, \qquad T_1 v = \tilde{A} v, \qquad T_{m+1} v = 2 \tilde{A}\, T_m v - T_{m-1} v,\]

so that only \(M\) products of \(A\) with a vector are needed.

The module provides four generic interfaces:

cheb_coefficients(): compute the coefficients \(c_m\) of a given function, optionally damped by the Jackson kernel;
cheb_set(): store the scale, the shift and the coefficients in the module memory;
cheb_projection(): evaluate the expansion, projected on selected components or vectors;
cheb_moments(): compute the individual terms \(T_m(\tilde{A})\,v\) of the expansion, projected in the same way, which do not depend on the coefficients.

Every interface exists for real(8) and for complex(8) data.

Note

The scale, the shift and the coefficients are kept in module-level variables shared by all the calls: the module is not thread safe. Invalid input is reported through the optional argument info of cheb_set(), cheb_projection() and cheb_moments(); if it is absent the program stops. cheb_coefficients() always stops on invalid input.

Example, with the spectrum of \(A\) in \([-4, 4]\) (gauss is a real function of one real variable and MatVec the matrix-vector product of \(A\)):

real(8) :: coeff(0:99), w(3)
call cheb_coefficients(gauss, 4d0, 0d0, coeff, kind="jackson")
call cheb_projection(MatVec, vseed, [1, 2, 3], w)

Quick access

Routines:: cheb_coefficients(), cheb_moments(), cheb_projection(), cheb_set()

Used modules

sf_constants: contains mathematical, numerical and physical constants
- pi
- zero
sf_iotools: SciFortran module for reading and writing data
- to_lower()
- str()

Subroutines and functions

interface sf_chebychev/cheb_coefficients(fun, a, b, coeff[, kind, tol, activate])

Chebyshev coefficients of a function, real or complex valued.

Computes the coefficients \(c_m\), \(m = 0, \dots, M\), of the expansion

\[f(x) \simeq \sum_{m=0}^{M} c_m\, T_m\!\left(\frac{x - b}{a}\right), \qquad x \in [b - a,\, b + a],\]

where \(M + 1\) is the size of the output array coeff.

The coefficients are obtained by Chebyshev-Gauss quadrature with \(n_q\) nodes,

\[c_m = \frac{2 - \delta_{m0}}{n_q} \sum_{j=1}^{n_q} f(a \cos\theta_j + b)\, \cos(m \theta_j), \qquad \theta_j = \frac{\pi (j - 1/2)}{n_q}.\]

Starting from \(n_q = \max(64,\, 2(M + 1))\), the number of nodes is doubled until two successive sets of coefficients differ by less than tol, relative to \(\max(1, \max_m |c_m|)\); the iteration is aborted after 128 doublings.

If kind selects the Jackson kernel, the coefficients are damped to reduce the Gibbs oscillations of the truncated series, \(c_m \to g_m c_m\), with

\[g_m = \frac{(M + 2 - m) \cos\frac{\pi m}{M + 2} + \sin\frac{\pi m}{M + 2} \cot\frac{\pi}{M + 2}}{M + 2}.\]

Unless activate is false, the scale, the shift and the coefficients are also stored in the module memory (see cheb_set()), ready for cheb_projection().

The function fun takes one real(8) argument; it returns a real(8) value for real coefficients and a complex(8) value for complex coefficients.

The program stops if kind is unknown, if tol is not positive, if a is not positive, if coeff is empty, or if the quadrature cannot be converged.

Parameters:

fun [external] – Real/complex function of one real variable, sampled on \([b - a,\, b + a]\)
a [real, in] – Scale: half-width of the interval \([b - a,\, b + a]\), must be positive
b [real, in] – Shift: center of the interval \([b - a,\, b + a]\)
coeff (1 + ) [real/complex, out] – Chebyshev coefficients \(c_0, \dots, c_M\); the size fixes the order \(M\) = size(coeff) - 1

Options:

kind [character(len=*)] – Damping kernel, case insensitive: "jackson" (or "j") or "none" (or "n"). Default "none"
tol [real, in] – Relative tolerance of the quadrature, must be positive. Default 1d-12
activate [logical, in] – If .true., a, b and coeff are stored in the module memory with cheb_set(). Default .true.

interface sf_chebychev/cheb_set(a, b, coeff[, coeff, info])

Store the scale, the shift and, optionally, the coefficients of the expansion in the module memory.

The scale \(a\), the shift \(b\) and the coefficients \(c_m\) of the expansion

\[f(A) \simeq \sum_{m=0}^{M} c_m\, T_m\!\left(\frac{A - b}{a}\right)\]

are kept in the module memory, from which cheb_projection() and cheb_moments() read them. cheb_coefficients() does it automatically; call cheb_set() directly to use coefficients computed elsewhere, or to change the scale and the shift only.

call cheb_set(a, b [, coeff, info])

With coefficients, they replace the stored ones. Complex coefficients are flagged as complex if at least one of them has a non-zero imaginary part: the real version of cheb_projection() then refuses to run.
Without coefficients (real version only), just the scale and the shift are set and any stored expansion is discarded. This is enough for cheb_moments().

Invalid input (non-positive scale, or empty coefficients) leaves the module memory unchanged: info is set to 1 or, if it is absent, the program stops.

Parameters:

a [real, in] – Scale of the map \(\tilde{A} = (A - b)/a\), must be positive
b [real, in] – Shift of the map \(\tilde{A} = (A - b)/a\)

Options:

coeff (1 + ) [real/complex, in] – Chebyshev coefficients \(c_0, \dots, c_M\) to store. If absent, any stored coefficients are discarded
info [integer, out] – Exit code: 0 on success, 1 if a is not positive or coeff is empty (the stored state is then unchanged). If absent, a failure stops the program

interface sf_chebychev/cheb_projection(matvec, vseed, rows, values, u[, info, numop])

Evaluate \(f(A)\,v\) with the stored Chebyshev expansion, projected on selected components or vectors.

Using the scale \(a\), the shift \(b\) and the coefficients \(c_m\) stored in the module memory (see cheb_coefficients() and cheb_set()), it computes

\[w = P f(A) v \simeq \sum_{m=0}^{M} c_m\, P\, T_m(\tilde{A})\, v, \qquad \tilde{A} = \frac{A - b}{a},\]

for the seed vector \(v\) (vseed), with \(M\) products of \(A\) with a vector, supplied by the callback MatVec. The projection \(P\) is selected by the third argument:

an integer array rows: \(w_k = [f(A)v]_{r_k}\), the components \(r_k\) of the result;
a matrix u: \(w_k = u_k^\dagger f(A) v\), the scalar product with the column \(u_k\) of \(U\) (conjugated for complex vectors).

call cheb_projection(MatVec, Vseed, rows, values [, info, NumOp])
call cheb_projection(MatVec, Vseed, U,    values [, info, NumOp])

With the real version the stored coefficients must be real.

On exit info is 0 on success, 1 for invalid input (inconsistent sizes, rows out of range, empty seed, or complex coefficients with real vectors) and 3 if the scale or the coefficients have not been set. If info is absent, a failure stops the program.

Parameters:

matvec [external] – Matrix-vector product subroutine, vout = A*vin, on vectors of length Nloc
vseed (•) [real/complex, in] – Seed vector \(v\), of size \(N\)
rows (•) [integer, in] – Indices, between 1 and \(N\), of the components of \(f(A)v\) to return
values (•) [real/complex, out] – Components \(w_k = [f(A)v]_{r_k}\), one for each entry \(r_k\) of rows
u (•, •) [complex/real, in] – Matrix of shape \(N \times K\) whose columns \(u_k\) are contracted with \(f(A)v\)

Options:

info [integer, out] – Exit code: 0 success, 1 invalid input, 3 scale or coefficients not set. If absent, a failure stops the program
numop [integer, out] – Number of matrix-vector products performed

interface sf_chebychev/cheb_moments(matvec, vseed, rows, moments, u[, info, numop])

Chebyshev moments of the seed vector, projected on selected components or vectors.

Computes the moments

\[\mu_{k,m} = [P\, T_m(\tilde{A})\, v]_k, \qquad \tilde{A} = \frac{A - b}{a}, \qquad m = 0, \dots, M,\]

for the seed vector \(v\) (vseed), where \(M + 1\) is the number of columns of moments. The projection \(P\) is selected as in cheb_projection():

an integer array rows: \(\mu_{k,m} = [T_m(\tilde{A}) v]_{r_k}\);
a matrix u: \(\mu_{k,m} = u_k^\dagger T_m(\tilde{A}) v\) (conjugated for complex vectors).

call cheb_moments(MatVec, Vseed, rows, moments [, info, NumOp])
call cheb_moments(MatVec, Vseed, U,    moments [, info, NumOp])

The moments do not depend on the coefficients: they are computed once and can be combined later with the coefficients of any function, \(P f(A) v \simeq \sum_m c_m \mu_{\cdot,m}\). The scale and the shift, however, must have been set (cheb_set(), even without coefficients).

On exit info is 0 on success, 1 for invalid input (inconsistent sizes, rows out of range, empty seed) and 3 if the scale and the shift have not been set. If info is absent, a failure stops the program.

Parameters:

matvec [external] – Matrix-vector product subroutine, vout = A*vin, on vectors of length Nloc
vseed (•) [real/complex, in] – Seed vector \(v\), of size \(N\)
rows (•) [integer, in] – Indices, between 1 and \(N\), of the components of \(T_m v\) to extract
moments (•, 1 + ) [real/complex, out] – Moments moments(k,m) = \([T_m(\tilde{A}) v]_{r_k}\); the first dimension must be the size of rows, and \(M\) = size(moments,2) - 1
u (•, •) [real/complex, in] – Matrix of shape \(N \times K\) whose columns \(u_k\) are contracted with \(T_m v\)

Options:

info [integer, out] – Exit code: 0 success, 1 invalid input, 3 scale and shift not set. If absent, a failure stops the program
numop [integer, out] – Number of matrix-vector products performed