expv.hpp

Go to the documentation of this file.

/// @file expv.hpp
/// @brief Krylov subspace matrix exponential-vector product: compute exp(t*A)*v
///
/// Approximates \f$\exp(tA)v \approx \|v\| Q_m \exp(tH_m)e_1\f$ where
/// \f$AQ_m \approx Q_{m+1}\bar{H}_m\f$ is the Arnoldi relation.
#pragma once
 
#include "core/matrix.hpp"
#include "core/types.hpp"
#include "core/vector.hpp"
#include "kernel/subspace.hpp"
#include "linalg/sparse/sparse.hpp"
#include "operator/operators.hpp"
#include <stdexcept>
#include <utility>
#include <vector>
 
namespace num {
 
namespace detail {
Matrix dense_expm_pade6(const Matrix& A);
}
 
/// @brief Compute \f$\exp(tA)v\f$ for any \f$y=A x\f$ adapter.
template<class Op>
    requires operators::LinearOperator<Op, Vector, Vector>
Vector expv(real t, const Op& A, const Vector& v, int m_max = 30, real tol = 1e-8) {
    const idx n = A.rows();
    if (A.cols() != n || v.size() != n) {
        throw std::invalid_argument("expv: dimension mismatch");
    }
 
    real beta = norm(v);
    if (beta < 1e-300) {
        return Vector(n, 0.0);
    }
 
    std::vector<Vector> V;
    V.reserve(m_max + 1);
 
    Vector v0(n);
    for (idx i = 0; i < n; i++) {
        v0[i] = v[i] / beta;
    }
    V.push_back(std::move(v0));
 
    Matrix H(m_max + 1, m_max, 0.0);
    int m_actual = m_max;
    std::vector<real> h_col(m_max + 1, 0.0);
 
    for (int j = 0; j < m_max; j++) {
        Vector w(n, 0.0);
        A.apply(V[j], w);
 
        const real h_next = kernel::subspace::mgs_orthogonalize(V, w, h_col, j + 1);
        for (int i = 0; i <= j; i++) {
            H(i, j) = h_col[i];
        }
        H(j + 1, j) = h_next;
 
        if (h_next < tol) {
            m_actual = j + 1;
            break;
        }
 
        scale(w, real(1) / h_next);
        V.push_back(std::move(w));
    }
 
    Matrix Hm(m_actual, m_actual, 0.0);
    for (int i = 0; i < m_actual; i++) {
        for (int j = 0; j < m_actual; j++) {
            Hm(i, j) = t * H(i, j);
        }
    }
 
    Matrix E = detail::dense_expm_pade6(Hm);
 
    Vector result(n, 0.0);
    for (int j = 0; j < m_actual; j++) {
        axpy(beta * E(j, 0), V[j], result);
    }
 
    return result;
}
 
Vector expv(real t,
            const SparseMatrix& A,
            const Vector& v,
            int m_max = 30,
            real tol = 1e-8);
 
} // namespace num