# LightKrylov (v0.1.0-beta) : Pre-release **URL:** \ **Category:** Announcements\ **Created:** [October 29, 2024, 1:58pm UTC](https://fortran-lang.discourse.group/t/lightkrylov-v0-1-0-beta-pre-release/8766 "2024-10-29T13:58:05Z")\ **Posts on this page:** 20\ **Page:** 1
**Author:** ![loiseaujc](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/loiseaujc/32/4582_2.png) [@loiseaujc](https://fortran-lang.discourse.group/u/loiseaujc)\ **Post date:** [October 29, 2024, 1:58pm UTC](https://fortran-lang.discourse.group/t/lightkrylov-v0-1-0-beta-pre-release/8766/1 "2024-10-29T13:58:05Z")
Hej everyone! I’m happy to announce that we’ve just officially released `LightKrylov v0.1.0-beta`! The code can be found [here](https://github.com/nekStab/LightKrylov) and the associated documentation [here](https://nekstab.github.io/LightKrylov/). It is a joint effort with two postdoctoral researchers in my lab (Simon Kern and Ricardo Frantz) and is funded partially thanks to an _Agence Nationale pour la Recherche_ project (the French equivalent of the NSF basically). I had already discussed quickly about `LightKrylov` on [this topic](https://fortran-lang.discourse.group/t/generic-or-kind-agnostic-linear-system-solvers/7574/1) but for those of you who have no idea what it is, here is the _too long/didn’t read_ summary: > `LightKrylov` is a lightweight implementation of standard Krylov-based techniques relying on modern `Fortran` features, in particular its `abstract type` capabilities. Having no dependencies other than `stdlib`, it exposes `abstract vector` and `abstract linop` types which you can easily extend and adapt to the particular data structure and/or parallelization paradigm you use. Having extending these two types to accomodate your particular applications, `LightKrylov` lets you solve linear systems, compute leading eigenvalues or singular values of your linear operator fairly easily using standard Krylov techniques such `gmres`, `arnoldi` or `lanczos`. The code is still a bit rough around the edges and we would be more than happy if any of you could give it a test drive and let us know any issue you would encounter or any improvement you’d like to see happen. We have a couple of examples at the moment, including computing the leading eigenpairs of the linearized Ginzburg-Landau operator, or using a Newton-Krylov algorithm to compute an unstable periodic orbit of the Roessler system in the chaotic regime and the associated Lyapunov vectors. We plan to add new examples shortly, in particular how to solve in parallel a Helmholtz equation (the `Hello World` of PDE) using the conjugate gradient method (with the preconditioned version being added quite soon). Note that we have already started to incorporate `LightKrylov` into [`neklab`](https://github.com/nekStab/neklab), a bifurcation and linear stability analysis toolbox we’re writing for the massively parallel spectral element solver `Nek5000`. Preliminary results on a slightly older version indicated that `LightKrylov` is competitive in terms of computational performances with `Arpack` for computing the leading eigenvalues of the linearized Navier-Stokes operator for the incompressible two-dimensional cylinder flow (roughly a million degrees of freedom on a dozen cores) while requiring far fewer lines of code and a much easier integration into the existing code base. `LightKrylov` is also the base package upon which we are building [`LightROM`](https://github.com/nekStab/LightROM), another package that will eventually makes its way into `neklab` to solve really high-dimensional Lyapunov and Riccati equations for linear optimal feedback control or estimation (i.e. LQR and LQG) using [Dynamical Low-Rank Approximation](https://epubs.siam.org/doi/10.1137/050639703). We’ll make a dedicated annoucement once `LightROM` has reached a sufficient level of maturity though. Finally, the low-level routines of `LightKrylov` rely quite extensively on `stdlib` and in particular `stdlib_linalg` and `stdlib_linalg_lapack`. As such, I’d like to personally thank @FedericoPerini, @hkvzjal and @jeremie.vandenplas for their amazing work on this! ---
**Author:** ![FedericoPerini](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/federicoperini/32/1750_2.png) [@FedericoPerini](https://fortran-lang.discourse.group/u/FedericoPerini)\ **Post date:** [October 29, 2024, 2:35pm UTC](https://fortran-lang.discourse.group/t/lightkrylov-v0-1-0-beta-pre-release/8766/2 "2024-10-29T14:35:19Z")
This is amazing @loiseaujc, shout out to you and your lab team! I have a personal interest in eigenvalues for large sparse problems (for chemical explosive mode analysis), in combustion it’s a very hard problem for packages like ARPACK because of the stiffness: the explosive modes typically have O(1) real components, while the most damping modes have magnitudes \>O(10^15) in the negative real range, so ARPACK can fail capturing the positive real modes because they are only found when the Krylov subspace is now almost full. On the other hand, the dense approach is painstakingly slow. So as soon as I have some time, I’ll take a look to see if LightKrylov can have applications there. On another note, shout out also to @jeremie.vandenplas and @hkvzjal for the linear algebra capability, it’s amazing to see that it’s already being effectively deployed. Please do not hesitate to contribute to it by reporting bugs, requesting features, etc. ---
**Author:** ![loiseaujc](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/loiseaujc/32/4582_2.png) [@loiseaujc](https://fortran-lang.discourse.group/u/loiseaujc)\ **Post date:** [October 29, 2024, 3:25pm UTC](https://fortran-lang.discourse.group/t/lightkrylov-v0-1-0-beta-pre-release/8766/3 "2024-10-29T15:25:57Z")
> [@FedericoPerini](#): > > Please do not hesitate to contribute to it by reporting bugs, requesting features, etc. Sure! Sorry about letting you down with the `norm` / `opnorm` stuff. Summer has been a bit overwhelming with the kids and I had a handful of new classes at the begining of this semester I had to prepare which prevented me from doing any serious work outside of `LightKrylov` essentially. We do have a handful of functions which could eventually make their ways into `stdlib_linalg` though. These include: - Matrix exponential using Pade approximation similar to `expm` in `SciPy`. - Schur factorization based on `gees`. - Re-ordered Schur factorization (i.e. `ordschur` in `Matlab`). I’m also working on a side project simultaneously: [`SpecialMatrices`](https://github.com/loiseaujc/SpecialMatrices/tree/main). It’ll essentially be very similar to [`SpecialMatrices.jl`](https://github.com/JuliaLinearAlgebra/SpecialMatrices.jl) and [`ToeplitzMatrices.jl`](https://github.com/JuliaLinearAlgebra/ToeplitzMatrices.jl), i.e. providing special types and computational routines which can be used as drop-in replacement for `stdlib_linalg` when working with `Diagonal`, `Bidiagonal`, `Tridiagonal`, `SymTridiagonal`, `Circulant`, `Toeplitz` or `Hankel` matrices. I basically plan to overload `matmul`, `solve`, `eig` / `eigh`, and a couple other things for each of these types to make it very easy to use. > [@FedericoPerini](#): > > the explosive modes typically have O(1) real components, while the most damping modes have magnitudes \>O(10^15) in the negative real range This is a fairly typical problem actually. Have you tried using a shift-invert strategy (albeit it does need a good preconditioned iterative solver)? Another alternative, which is typically what we use for Navier-Stokes is to use a time-stepping approach. Say the matrix you are interested in is A. You can associate to it a linear-time invariant dynamical system \dot{x} = Ax. Its solution reads x(\tau) = \exp(\tau A) x\_0. The eigenvalues \mu\_i of the exponential propagator are related to those of A via the exp-transform \mu\_i = \exp(\tau \lambda\_i) (and vice versa \lambda\_i = \log(\mu\_i) / \tau). Hence, eigenvalues corresponding to very stable modes get maps to zero, while those closer to the imaginary axis get maps close to the unit circle. If you then use `arnoldi` on the exponential propagator rather than A iteself, the first eigenvalues to converge are actually the ones you’re looking for. Note that you never need to form \exp(\tau A), you simply need its action onto a vector. This can easily be computed using your favorite time-integration scheme and this where iterative eigenvalue solvers shine and what you would typically encapsulate in the extended `abstract linop` in `LightKrylov`. How many Krylov vectors need to be generated to converge the eigenvalues obviously depends on \tau which itself impact how time-consuming computing \exp(\tau A)x is, but a rough estimate is to choose \tau compatible with the dominant time-scale (either the damping time or the period if it is oscillatory) of the leading eigenvalues you want to converge. If you look at the Ginzburg-Landau example in the `LightKrylov`, this is actually how we’ve done it. Likewise, this time-stepper appraoch is our go-to approach when dealing with the linearized Navier-Stokes operator. We have written a review paper on this [here](https://arxiv.org/pdf/2301.12940). We can setup a call sometime if you want to know more. ---
**Author:** ![Abhi\_dev](https://avatars.discourse-cdn.com/v4/letter/a/2bfe46/32.png) [@Abhi\_dev](https://fortran-lang.discourse.group/u/Abhi_dev)\ **Post date:** [October 29, 2024, 3:42pm UTC](https://fortran-lang.discourse.group/t/lightkrylov-v0-1-0-beta-pre-release/8766/4 "2024-10-29T15:42:29Z")
Is there a way to track the progress of stdlib\_linalg module. What are the features that are implemented, what are currently in progress and the roadmap so that we can contribute ? ---
**Author:** ![rsci](https://avatars.discourse-cdn.com/v4/letter/r/e480ec/32.png) [@rsci](https://fortran-lang.discourse.group/u/rsci)\ **Post date:** [October 29, 2024, 5:43pm UTC](https://fortran-lang.discourse.group/t/lightkrylov-v0-1-0-beta-pre-release/8766/5 "2024-10-29T17:43:21Z")
I took a brief look at the code, and it looks great! Regarding massively parallel computing, are there any plans to implement communication-avoiding variants of the CG or gmres solvers? ---
**Author:** ![FedericoPerini](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/federicoperini/32/1750_2.png) [@FedericoPerini](https://fortran-lang.discourse.group/u/FedericoPerini)\ **Post date:** [October 30, 2024, 9:01am UTC](https://fortran-lang.discourse.group/t/lightkrylov-v0-1-0-beta-pre-release/8766/6 "2024-10-30T09:01:02Z")
Wow, thanks a lot @loiseaujc, this is great insight, thank you! It’s a very interesting idea, I had been using [EXPOKIT](https://www.maths.uq.edu.au/expokit/download.html) in the past for this purpose, but only for dense matrices. It seems like there are also routines for sparse matrices using spMv, I’ll definitely dive deeper into `LightKrylov` when I have more time. BTW I appreciate that by integrating stdlib into your work, you haven’t found huge issues yet 🤣 Also: @hkvzjal and I have been working on the matrix norms and I should open a PR for that shortly. The Schur decomposition would definitely be another great one for the immediate/short term. Best of luck with `LightKrylov`, it’s great to see math abstractions coming to life in Modern Fortran, and thanks again for sharing your insights! ---
**Author:** ![FedericoPerini](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/federicoperini/32/1750_2.png) [@FedericoPerini](https://fortran-lang.discourse.group/u/FedericoPerini)\ **Post date:** [October 30, 2024, 9:02am UTC](https://fortran-lang.discourse.group/t/lightkrylov-v0-1-0-beta-pre-release/8766/7 "2024-10-30T09:02:06Z")
Welcome to the forum @Abhi_dev, thanks for the interest! I would think that the best way to open a discussion for further features is to open an issue on the Github page, or if you have a clean implementation, you could open a PR directly. ---
**Author:** ![loiseaujc](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/loiseaujc/32/4582_2.png) [@loiseaujc](https://fortran-lang.discourse.group/u/loiseaujc)\ **Post date:** [October 30, 2024, 11:30am UTC](https://fortran-lang.discourse.group/t/lightkrylov-v0-1-0-beta-pre-release/8766/8 "2024-10-30T11:30:10Z")
> [@FedericoPerini](#): > > Wow, thanks a lot @loiseaujc, this is great insight, thank you! It’s a very interesting idea, I had been using [EXPOKIT](https://www.maths.uq.edu.au/expokit/download.html) in the past for this purpose, but only for dense matrices. It seems like there are also routines for sparse matrices using spMv, I’ll definitely dive deeper into `LightKrylov` when I have more time. My pleasure! In general, if you want to solve Ax = \lambda x using Krylov-based techniques but the eigenvalues of interest are not on the outer rim of the spectrum, you can then replace it with f(A)x = f(\lambda) x where f is an invertible function which maps the eigenvalues of interest to the outer rim of the eigenspectrum. These will then be the first to converge. The critical point however is to find a function f which can easily be evaluated. - If you have a really good preconditioner for A, defining f(\lambda) = \dfrac{1}{\lambda - \sigma} where \sigma \in \mathbb{C} is your shift is probably the best option. It corresponds to the standard shift-invert strategy. It has some limitations though: not only do you need a good preconditioner, but your code also need to handle complex arithmetic. Additionally, if you change A you may have to completely rethink your preconditioning strategy. It definitely is worth the time if you can, but it’ll come at the price of extended code development time. - Using f(\lambda) = \exp(\tau \lambda) implies that you actually work with the exponential propagator \exp(\tau A). As I’ve said, you never want to actually compute this matrix (which will be dense even if A is sparse and requires \mathcal{O}(n^3) operations anyway). This is where time-integration schemes come into play to approximate it. - Say you use a simple first order explicit Euler scheme with constant time step \Delta t = \tau/n. Then, you’re actually approximating the exponential with the following polynomial function f(\lambda) \simeq \hat{f}(\lambda) = (1 + \delta t \lambda)^n. - If you use a first order implicit Euler scheme, then \hat{f}(\lambda) = \dfrac{1}{(1 - \delta t \lambda)^n}. - If you use a semi-implicit scheme, then you approximate f with the rational function \hat{f}(\lambda) = \dfrac{(1 + \frac{\delta t}{2} \lambda)^n}{(1 - \frac{\delta t}{2} \lambda)^n}. - If you use an adaptive Runge-Kutta scheme, you get yet another polynomial approximation of the exponential, etc. Another alternative are indeed Krylov-based exponential integrators such as the ones in `EXPOKIT`. If my memory serves me right, it has been proven that such Krylov-based integrators are actually near-optimal in the sense that, for a polynomial approximation and a given accuracy, they actually minimize the number of matrix-vector products you need to do to evaluate \exp(\tau A) x. Note however that the implementation in `LightKrylov` is still pretty experimental and nowhere as feature-complete as `EXPOKIT`. We do have plans later on to actually incorporate time-step adaptive capabilities balancing how expensive the matrix-vector product is with how good the Krylov approximation is. I’m not that familiar with `EXPOKIT` but I wouldn’t be surprised if they already had such features implemented. In any case, at the end of the day, the choice of the time-integration scheme is not the primal factor from my point of view. The major benefit is that, if you already have a (linearized) simulation code where you spent some time optimizing the time-integration, then using this time-stepping approach to do linear stability analysis requires actually very little modification to the code base. These ideas have been developed since the 1980’s by Laurette Tuckerman and others, and the _time-stepper approach to linear stability_ has become the go-to in the hydrodynamic stability community. It was actually the main theme of my own PhD thesis. ---
**Author:** ![loiseaujc](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/loiseaujc/32/4582_2.png) [@loiseaujc](https://fortran-lang.discourse.group/u/loiseaujc)\ **Post date:** [October 30, 2024, 12:09pm UTC](https://fortran-lang.discourse.group/t/lightkrylov-v0-1-0-beta-pre-release/8766/9 "2024-10-30T12:09:30Z")
At the moment we do not have such plans, even though it is somewhere in the back of our minds. Our main goal for now is to make sure that the standard algorithms in `LightKrylov` work as expected, can easily be incorporated into an existing code-base with minimal effort and have great baseline performances. I would also have to look more closely at these algorithms and see exactly how they can be implemented using the `abstract type` route we’ve decided to take. Note however that for solving Ax = b with A symmetric positive definite, we are currently working with a colleague on a high-performing implementation of the Chebyshev-accelerated Jacobi (and block-Jacobi) method, in particular for the finite-difference approximation of the Laplace operator in 2 and 3 dimensions. While the convergence rate might be slightly less good than conjugate gradient, the main benefit of the Chebyshev-Jacobi method is that it does not require inner product whatsoever which actually is the main bottleneck for parallel CG as it requires a communication which synchronizes the processes. You do need halo exchanges though, but these are also present in CG when evaluating Ax but can be done with non-blocking communications. What we expect is that, despite the slightly less good convergence rate, the wall-clock time-to-solution of the Chebyshev-Jacobi be better than CG precisely because of the absence of synchronizing inner products. I’ll keep everyone updated once we have the code ready. ---
**Author:** ![rsci](https://avatars.discourse-cdn.com/v4/letter/r/e480ec/32.png) [@rsci](https://fortran-lang.discourse.group/u/rsci)\ **Post date:** [October 30, 2024, 2:38pm UTC](https://fortran-lang.discourse.group/t/lightkrylov-v0-1-0-beta-pre-release/8766/10 "2024-10-30T14:38:37Z")
In our atmospheric model code, we have already conducted some preliminary comparisons between CG and Chebyshev methods. The asymptotic convergence of both is similar, but Chebyshev iterations scale much better. However, Chebyshev methods require accurate estimates of the upper bound of the spectrum; even a small underestimation can lead to a noticeable drop in the convergence rate or even divergence of the method. We also have a geometric multigrid method that demonstrates excellent convergence, but its parallel efficiency is limited by the scalability of the coarse levels, moreover constructing an effective smoother implicitly requires some knowelgde about the matrix spectrum (one needs to optimize convergence for the [\lambda\_s, \lambda\_{max}] part of the spectrum, where \lambda\_s is the boundary between slowly and fast-varying eigenvectors). Therefore, I am considering trying communication-avoiding Krylov methods, which, in theory, allow for an arbitrary number k of iterations of the method within a single global communication. In practice, it seems that using large values of k leads to numerical stability issues (since one needs to operate on A^k v vectors without orthogonalization), so some additional care is needed. It’s great that this forum allows us to discuss not only Fortran but also relevant scientific issues:) Good luck with your project, and I look forward to hearing about your progress. ---
**Author:** ![hkvzjal](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/hkvzjal/32/3055_2.png) [@hkvzjal](https://fortran-lang.discourse.group/u/hkvzjal)\ **Post date:** [October 30, 2024, 3:20pm UTC](https://fortran-lang.discourse.group/t/lightkrylov-v0-1-0-beta-pre-release/8766/11 "2024-10-30T15:20:08Z")
Nice work @loiseaujc I’ll definitely take a close look! > [@rsci](#): > > are there any plans to implement communication-avoiding variants of the CG or gmres solvers? > [@loiseaujc](#): > > Chebyshev-Jacobi method is that it does not require inner product whatsoever which actually is the main bottleneck for parallel CG as it requires a communication which synchronizes the processes. Indeed, my understanding is also that trying to “asynchronize” methods with a strong dependency on an accurate computation of the inner-product (CG or gmres) is not a good idea as convergency can not be guaranteed. If this a domain of interest, you might want to follow the works of Prof. Frédéric Magoulès [https://www.researchgate.net/profile/Frederic-Magoules/publications](https://www.researchgate.net/profile/Frederic-Magoules/publications) who has extensively worked on robust asynchronous solvers for which convergence can be demonstrated and ensured at user-defined tolerance [Google Scholar](https://scholar.google.com/scholar?hl=en&as_sdt=0%2C5&q=Frederic+Magoules+Asynchronous&btnG=). ---
**Author:** ![loiseaujc](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/loiseaujc/32/4582_2.png) [@loiseaujc](https://fortran-lang.discourse.group/u/loiseaujc)\ **Post date:** [October 30, 2024, 4:32pm UTC](https://fortran-lang.discourse.group/t/lightkrylov-v0-1-0-beta-pre-release/8766/12 "2024-10-30T16:32:50Z")
> [@rsci](#): > > The asymptotic convergence of both is similar, but Chebyshev iterations scale much better. Cool! Glad to know my understanding is corroborated numerically! > [@rsci](#): > > However, Chebyshev methods require accurate estimates of the upper bound of the spectrum; even a small underestimation can lead to a noticeable drop in the convergence rate or even divergence of the method. That is indeed a major limitation compared to CG. It is also one of the reasons we’re focusing on the Laplace operator on standard domains (very academic admittedly): the spectral radius for finite-difference approximations on hyper-cubes can be computed analytically. For more complicated domains, we actually have simple pre-processing technique where we run the standard Jacobi method for a number of steps and infer its spectral radius from the convergence history and use that as the upper bound. This simple heuristic works suprisingly well, albeit you do need to have this pre-processing step. Doing so make sense if you need to solve Ax = b for multiple right-hand side vectors as the divergence-free projection step in an incompressible Navier-Stokes simulation. If A changes at every iteration, this is no longer practical. > [@rsci](#): > > Therefore, I am considering trying communication-avoiding Krylov methods, which, in theory, allow for an arbitrary number k of iterations of the method within a single global communication. In practice, it seems that using large values of k leads to numerical stability issues (since one needs to operate on A^k v vectors without orthogonalization), so some additional care is needed. > [@hkvzjal](#): > > Indeed, my understanding is also that trying to “asynchronize” methods with a strong dependency on an accurate computation of the inner-product (CG or gmres) is not a good idea as convergency can not be guaranteed. Yep. As far as I know, the main problem of CA-GMRES is the loss of orthogonality of the Krylov vectors. An alternative we’re considering to implement in `LightKrylov` though is the use polynomial-Arnoldi preconditioners. You’re still using a high-degree polynomial but maintain at the same time the orthogonality of the Krylov basis. You can look at some papers by Mark Embree if you want to know more. This is definitely at the top of our list of things we’d like to include in `LightKrylov`. (edit) PS: I will actually take some time tomorrow after I’ve discussed with Simon to come up with a list of the top features we’d like to include in `LightKrylov`, just to give ideas if any one wants to give it a go. ---
**Author:** ![CRquantum](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/crquantum/32/730_2.png) [@CRquantum](https://fortran-lang.discourse.group/u/CRquantum)\ **Post date:** [November 1, 2024, 11:15pm UTC](https://fortran-lang.discourse.group/t/lightkrylov-v0-1-0-beta-pre-release/8766/13 "2024-11-01T23:15:07Z")
If you use EXPOKIT, note that there are some small bugs in some of its code which can cause integer overflow and lead to wrong results. For example, in dgpadm.f, around line 87 ![image](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/2X/5/513dfc724c3e057ae5ab250be065930bba604a2d.png) Originally it is ```auto scale = t / DBLE(2**ns) ``` that ` 2**ns` will cause integer overflow if ns\>30 if the default integer type is integer 4, then the DBLE(2\*\*ns) will be something like NAN. You need to change `DBLE(2**ns)` to something like `2.0d0**ns` Without this fix, every time if ns\>30, dgpadm will just output wrong results like NAN or something. Basically any `DBLE()` in those code in EXPOKT need to be carefully checked to prevent integer overflow issues. ---
**Author:** ![FedericoPerini](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/federicoperini/32/1750_2.png) [@FedericoPerini](https://fortran-lang.discourse.group/u/FedericoPerini)\ **Post date:** [November 2, 2024, 11:24am UTC](https://fortran-lang.discourse.group/t/lightkrylov-v0-1-0-beta-pre-release/8766/14 "2024-11-02T11:24:50Z")
Great catch @CRquantum, thanks! Do you know if the package is still being maintained somehow? It would be nice to use a modernized (and possibly bugfixed) version ---
**Author:** ![CRquantum](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/crquantum/32/730_2.png) [@CRquantum](https://fortran-lang.discourse.group/u/CRquantum)\ **Post date:** [November 3, 2024, 2:29am UTC](https://fortran-lang.discourse.group/t/lightkrylov-v0-1-0-beta-pre-release/8766/15 "2024-11-03T02:29:17Z")
Thanks @FedericoPerini Perhaps it is best to directly contact the EXPOKIT author Roger B. Sidje? > **[Roger B. Sidje – Mathematics](https://math.ua.edu/people/roger-b-sidje/)** > **[Roger B Sidje](https://scholar.google.com/citations?user=tGKL1jUAAAAJ&hl=fr)** > > University of Alabama - Cité(e) 2 305 fois PS. Besides expokit, I found some latest matrix expontial solvers seems working good too, such as those by Jorge Sastre, > **[Jorge Sastre](https://scholar.google.com/citations?hl=en&user=Gls8MJwAAAAJ&view_op=list_works&sortby=pubdate)** > > Full Professor (Universitat Politècnica de València, Spain) - Citado por 1,334 - Applied Mathematics - Music Technology - Video processing > **[Two Taylor Algorithms for Computing the Action of the Matrix Exponential on a...](https://www.mdpi.com/1999-4893/15/2/48)** > > The action of the matrix exponential on a vector eAtv, A∈Cn×n, v∈Cn, appears in problems that arise in mathematics, physics, and engineering, such as the solution of systems of linear ordinary differential equations with constant coefficients.... > **[Boosting the computation of the matrix exponential](https://scholar.google.com/citations?view_op=view_citation&hl=en&user=Gls8MJwAAAAJ&cstart=20&pagesize=80&sortby=pubdate&citation_for_view=Gls8MJwAAAAJ:RYcK_YlVTxYC)** > > J Sastre, J Ibáñez, E Defez, Applied mathematics and computation, 2019 - Cited by 39 Usually in his paper he provide his solver’s matlab (can be easily tranlated into Fortran) or Fortran version. ---
**Author:** ![loiseaujc](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/loiseaujc/32/4582_2.png) [@loiseaujc](https://fortran-lang.discourse.group/u/loiseaujc)\ **Post date:** [February 17, 2025, 1:09pm UTC](https://fortran-lang.discourse.group/t/lightkrylov-v0-1-0-beta-pre-release/8766/16 "2025-02-17T13:09:44Z")
Quick update : I’ve just implemented the polynomial arnoldi preconditioner for GMRES to try it out on the 2D Poisson equation. [Here](https://github.com/loiseaujc/polynomial_arnoldi) is the repo. Although some bits and pieces are still a bit rough around the edges, it works fantastically well. Not only is the number of matrix-vector products reduced compared to the unpreconditioned one, but the number of dot products that need to be computed is even further decreased! Haven’t had the chance to try it out in parallel, but I’m pretty sure it’ll have pretty good scaling. I’ll port it to [LightKrylov](https://github.com/nekStab/LightKrylov) shortly. ---
**Author:** ![hkvzjal](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/hkvzjal/32/3055_2.png) [@hkvzjal](https://fortran-lang.discourse.group/u/hkvzjal)\ **Post date:** [February 18, 2025, 8:19am UTC](https://fortran-lang.discourse.group/t/lightkrylov-v0-1-0-beta-pre-release/8766/17 "2025-02-18T08:19:30Z")
> [@loiseaujc](#): > > Not only is the number of matrix-vector products reduced compared to the unpreconditioned one, but the number of dot products that need to be computed is even further decreased Interesting, just out of curiosity, will you compare against a Jacobi and/or Incomplete Cholesky preconditioner? ---
**Author:** ![rsci](https://avatars.discourse-cdn.com/v4/letter/r/e480ec/32.png) [@rsci](https://fortran-lang.discourse.group/u/rsci)\ **Post date:** [February 18, 2025, 9:22am UTC](https://fortran-lang.discourse.group/t/lightkrylov-v0-1-0-beta-pre-release/8766/18 "2025-02-18T09:22:42Z")
Nice work! Just to clarify, are you using the GMRES polynomial as a preconditioner for restarted GMRES? If there’s any reference for the approach you implemented, could you share a link? I’d like to take a closer look! ---
**Author:** ![loiseaujc](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/loiseaujc/32/4582_2.png) [@loiseaujc](https://fortran-lang.discourse.group/u/loiseaujc)\ **Post date:** [February 21, 2025, 10:22am UTC](https://fortran-lang.discourse.group/t/lightkrylov-v0-1-0-beta-pre-release/8766/19 "2025-02-21T10:22:30Z")
The nice thing about this GMRES polynomial preconditioner is that, if you already use another preconditioner (say block-Jacobi or incomplete Cholesky for instance), you can actually use B = A M^{-1} to further improve the already preconditioned system. Unfortunately, this is a pet project for me so I can’t dedicate much time. I do plan however to compare the performances against block-Jacobi preconditioner, Chebyshev-Jacobi acceleration, and possibly a simple geometric multigrid solver if I have time to implement it. The Bartels-Stewart algorithm for solving AX + XB = C is also on the table as you can write the 2D Poisson problem on the unit square into this Sylvester equation form. All of these will only be serial computations though, so whenever I post the result, you’ll need to take it with a grain of salt. ---
**Author:** ![loiseaujc](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/loiseaujc/32/4582_2.png) [@loiseaujc](https://fortran-lang.discourse.group/u/loiseaujc)\ **Post date:** [February 21, 2025, 10:33am UTC](https://fortran-lang.discourse.group/t/lightkrylov-v0-1-0-beta-pre-release/8766/20 "2025-02-21T10:33:34Z")
Indeed, the preconditioner is built from the GMRES polynomial. Provided you have a routine to compute the k-step Arnoldi factorization, you’re good to go. Here are a couple of references: - M. Embree, J. A. Loe, and R. Morgan. _Polynomial preconditioned Arnoldi with stability control_, SIAM Journal on Scientific Computing, vol. 43, 2021 [[link]](https://epubs.siam.org/doi/abs/10.1137/19M1302430). - J. Loe, and R. Morgan. _Toward efficient polynomial preconditioning for GMRES_, Numerical linear algebra with applications, vol. 28, 2022 [[arxiv]](https://onlinelibrary.wiley.com/doi/abs/10.1002/nla.2427) - M. Embree, J. Henningsen, J. Jackson and R. Morgan. _Polynomial Approximation to the inverse of a large matrix_ [[PDF]](https://bpb-us-w2.wpmucdn.com/sites.baylor.edu/dist/e/71/files/2024/11/Polynomial_Inverse_Approximation.pdf) - J. Henngingsen. _A polynomial approximation to the inverse of a matrix_. [[PhD thesis]](https://baylor-ir.tdl.org/server/api/core/bitstreams/a2acc8b6-11df-45e8-9581-e00bdf892ec5/content) The code in the repo is mostly based on Loe & Morgan (second reference). I’m in the process of reading the last two references.