Type promotion

[pmli]

This pull request is to enable automatic type promotion and correct handling of complex types. This is as discussed in #333, except unlike there, I decided not to include special handling when imaginary part is zero. The main reason is that I noticed that such conditions can be handled in "higher-level" functions (e.g. reductors). Also, it just seems like too much work to do it in all "lower-level" function (e.g. operator and vectorarray methods)...

So far, I only made some changes in vectorarrays and operators. @pymor/pymor-devs @andreasbuhr, please check the changes, and see if additional changes are necessary.

Besides vectorarrays and operators, I noticed one thing in algorithms which should be changed. In line 341 in algorithms/genericsolvers.py, the assignment doesn't work for complex numbers. Do we want to touch LGMRES code?

I haven't noticed any other ComplexWarnings or TypeErrors (like in #333), but in future, we should probably include tests with complex data (and mixtures or float and complex data). Is this something easy to enable in the current testing framework?

There is also a general question of handling complex numbers, like in dot (see #259). Should this be a separate pull request?

[sdrave]

Hi @pmli, I see you have added the dtype argument to VectorArrayInterface.zeros and not only to NumpyVectorArray.zeros. Could you elaborate why/where this is necessary? Almost no PDE-solver will support specifying the dtype (typically, the type of vector you can have will be determined by the discrete function space used in the solver), so specifying it will usually case an error.

[pmli]

Hi @sdrave. I guess it's not necessary.

[pmli]

@sdrave Can VectorSpaceInterface.zeros also be without dtype argument?

[sdrave]

I would say it is perfectly fine to add additional optional arguments to specializations. So, yes. Unless there is a reason to have this for arbitrary backends, I would only add this to NumpyVectorArray.
(Not sure how we want to proceed regarding solving shifted systems (A + sE)x = f for PDE discretizations ..)

[pmli]

That's a good point about solving shifted systems. Do any PDE libraries support complex valued arrays? Is it realistic that they would add it? If not, I guess the alternative is rewriting complex systems as real systems of double the size.

[sdrave]

I did a little research: FEniCS does not support complex vectors, deal.II supports complex vectors but recommends assembling as a real system, DUNE supports complex vectors in general, but this is not very well tested.

[pmli]

So, should I put dtype everywhere, or only in NumpyVectorArray (and we deal with other later)?

[sdrave]

For now, I would say, only add it to NumpyVectorArray.

[sdrave]

If I remember our discussion correctly, zeros should probably have no dtype argument and there should be no dtype=None in the subspace.zeros call?

[pmli]

I guess we agreed on your option 2 in #333, but then limited dtype in zeros to NumpyVectorArray in this thread. Then I realized that it would make sense that BlockVectorArray.zeros also has dtype if it is created from NumpyVectorArrays. Now I see that this code would not work for arbitrary BlockVectorArray.
I'm not sure what to do here. Should I add a check isinstance(subspace, NumpyVectorSpace), or add dtype everywhere, or remove it from everywhere?

[sdrave]

Then I realized that it would make sense that BlockVectorArray.zeros also has dtype if it is created from NumpyVectorArrays.

Could you give an example?

[pmli]

Yes, sure. In sys-mor branch, algorithms.sylvester creates an empty vectorarray from A.source. Operator A can be a BlockOperator, e.g. when transforming a second order system into a first order system (see discretizations.iosys.SecondOrderSystem.to_lti).

[sdrave]

Ok, to ask the other way around: Are there situations where it would be crucial to have a dtype argument for zeros or empty? As far as I can see, there are mainly two possibilities:

1. automatic type promotion is not allowed for some reason
2. automatic type promotion would be to expensive (probably in terms of memory)

[pmli]

Another possibility is implementing mixed precision algorithms, e.g. switching between single and double precision, but that will probably not happen soon, if ever.
I don't know if automatic type promotion could be prohibitively expensive in practice. So, I guess it is not very crucial to have a dtype argument.

[sdrave]

I would assume that for mixed precision algorithms you would want to have good control where which precision is used, so the precision should probably be part of the VectorSpace. - Since copies are usually relatively cheap, I would agree to remove the dtype argument until we see a specific need for the functionality.

[pmli]

Ok, good, I will work on the branch this week.

Could you please take a look at the questions in the first post here?

[sdrave]

Hi @pmli, regarding your questions:

Besides vectorarrays and operators, I noticed one thing in algorithms which should be changed. In line 341 in algorithms/genericsolvers.py, the assignment doesn't work for complex numbers. Do we want to touch LGMRES code?

I would say, we should not put any effort into this for now. If we would change this, we should write proper tests. And since I am not aware of anyone really using the code right now, this feels like a waste of effort.

I haven't noticed any other ComplexWarnings or TypeErrors (like in #333), but in future, we should probably include tests with complex data (and mixtures or float and complex data). Is this something easy to enable in the current testing framework?

Well, the easiest thing to try would be to add some VectorArray fixtures with complex data in pymortests.fixtures.vectorarray. If I haven't overlook anything, type promotion should ensure all tests to pass.

There is also a general question of handling complex numbers, like in dot (see #259). Should this be a separate pull request?

Yes. I believe we are still waiting for @andreasbuhr, to fully work out a proposal ...

[sdrave]

Hi @pmli, I had another look at this PR. Looks fine to me, I only added type Promotion also to NumpyVectorArrayView. Merge if you believe this is ready.