search

acroy / QuDOS.jl

Quantum Dynamics of Open Systems in Julia
Other
8 stars 3 forks source link

Roadmap/Wishlist #1

Open acroy opened 10 years ago

acroy commented 10 years ago

This issue is supposed to give a rough idea about the features/functionality I am planning to implement and would like to have (not necessarily in this package). Any help and comments are greatly appreciated. Apart from the points below there is of course still a lot of work to do regarding testing, optimization and documentation.

(The references given above are supposed to give examples for the respective method, they are by no means complete)

Jutho commented 10 years ago

Very nice. I very much like that you have implemented the Krylov method for the matrix exponential (times vector). I implemented this at some point as well, and it would be good to have this available in some basic Julia package. I agree that IterativeSolvers.jl would be a good place, although it would require some refactoring to fit it within the framework that is used in that package (although I am not certain that their framework is ideal for maximal efficiency, but it might not be too late to suggest some changes).

acroy commented 10 years ago

Thanks. Originally I was thinking to integrate expmv into ODE.jl, but maybe IterativeSolvers.jl is a better place. nulls might also be of general interest (although I am not entirely satisfied with it).

acroy commented 10 years ago

@Jutho : Thinking about it, if there is some interest in this function I guess expmv could also be part of a separate package (Expokit.jl?). Then people can test it and it could be refactored/optimized to fit into IterativeSolvers.jl or whatever.

Jutho commented 10 years ago

That sounds like a nice idea. I have no experience whatsoever with licensing questions etc. I guess one could email the authors and ask whether this is allowed and under what license ‘expokit.jl’ can be distributed.

I would at some point like to contribute to IterativeSolvers.jl and hopefully reshape it a little bit and add some eigensolver methods. Then I would also like to add expmv. However, I don’t currently have the time for this.

On 24 Jul 2014, at 12:20, acroy notifications@github.com wrote:

@Jutho : Thinking about it, if there is some interest in this function I guess expmv could also be part of a separate package (Expokit.jl?). Then people can test it and it could be refactored/optimized to fit into IterativeSolvers.jl or whatever.

— Reply to this email directly or view it on GitHub.

acroy commented 10 years ago

My knowledge about licensing is also very limited, but Expokit comes with a copyright notice which seems to be largely compatible with the MIT license. Putting this notice into LICENSE.md should be sufficient. But I can also write to Roger Sidje and ask him explicitly. To me it is a bit unclear if an implementation in Julia actually qualifies as a modification of Expokit like it would for GPLed code ... we will see.

Once it is set up I can add you as collaborator (if you want). Maybe we can also find someone who implements phiv which I found useful in the past.

acroy commented 10 years ago

More or less done: Expokit.jl. I added you as a collaborator. We might want to change the interface slightly before we announce it (for instance I think the current version is overtyped)?

Jutho commented 10 years ago

Looks great. Thanks a lot. Did you contact/inform the authors of Expokit? I might look at the overtyping somewhere this week, as I might already need to make use of this :-).

acroy commented 10 years ago

Yes, I contacted Roger and he basically said that we are "free to use Expokit in any way we want". So I chose the MIT license for this. It would be great if you could look into the typing, since I am also not sure if the function is type-stable the way it is now (well, I am more or less certain it is not ;-) I won't be able to do much for 1-2 weeks since we are moving, so feel free to push things forward ... (I also set up Travis which hopefully works after the next push to the repo)

i2000s commented 10 years ago

In practice, I think a better way to implement the expmv() function is to use the Al-Mohy's algorithm (A. H. Al-Mohy and N. J. Higham. Computing the action of the matrix exponential, with an application to exponential integrators. MIMS EPrint 2010.30, The University of Manchester, 2010; to appear in SIAM J. Sci. Comput.). The MATLAB implementation and a rewrite in C can be found here and here. I personally use this algorithm in MATLAB, and it is very fast.

Jutho commented 10 years ago

I am not an export on this, and I just quickly skimmed through the paper. From Table 6.4 in the paper, I see that this algorithm actually uses more matrix vector products than expv, but that the matlab implementation of expv takes more time because of the Arnoldi iteration. This is because Matlab is absolutely unsuitable to write something like Lanczos or Arnoldi efficiently. In Julia, on the other hand, the cost of the Arnoldi iteration should be negligible compared to the matrix vector product, provided it is well implemented (i.e. avoids allocation in loops etc).

@acroy, I know I promised to look at your expokit.jl implementation, but I haven’t had the time yet to look at it in detail. I remember that it looked very good, but that I thought that a few more allocations could be avoided. I need an exponential integrator at the moment, but only for a Hermitian matrix, so I took my old expv implementation using Lanczos iteration instead.

On 10 Aug 2014, at 23:17, Qi notifications@github.com wrote:

In practice, I think a better way to implement the expmv() function is to use the Al-Mohy's algorithm (A. H. Al-Mohy and N. J. Higham. Computing the action of the matrix exponential, with an application to exponential integrators. MIMS EPrint 2010.30, The University of Manchester, 2010; to appear in SIAM J. Sci. Comput.). The MATLAB implementation and a rewrite in C can be found here and here. I personally use this algorithm in MATLAB, and it is very fast.

— Reply to this email directly or view it on GitHub.

acroy commented 10 years ago

Certainly interesting. It might be interesting to compare implementations in Julia of the two methods, but I tend to agree with Jutho. Unfortunately I wont have time to look into it. @i2000s : Would you mind to open a new issue (here or at Expokit.jl) so that this suggestion doesn't get lost?

@Jutho : Don't worry, there is no need to rush. If you find the time, maybe you could test Expokit.expv against your expv?

acroy commented 10 years ago

@Jutho: I just pushed some changes to a new branch at Expokit.jl which should reduce the allocations. Maybe you could have a brief look?

Jutho commented 10 years ago

The problem might be that there are still a lot of constructions like vm[:,j] which create copies. I think slice(vm,:,j) should be better. Adding 2 vectors also allocates a new vector. Check out Base.LinAlg.axpy! which is extremely useful for writing Arnoldi iterations with little extra overhead.

On 11 Aug 2014, at 11:46, acroy notifications@github.com wrote:

@Jutho: I just pushed some changes to a new branch at Expokit.jl which should reduce the allocations. Maybe you could have a brief look?

— Reply to this email directly or view it on GitHub.

acroy commented 10 years ago

Thanks! Did you look at the opt branch? There I tried to avoid vm[:,j] etc. But you are right, axpy! gave even more speedup.