Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
Proposed algorithm
We an algorithm for the reciprocal scaling of a complex vector
Xby a complex numberA. We use the notationARandAIto denote the real and imaginary parts ofA, respectively. The algorithm is as follows:AIis zero, then call{CS,ZD}RSCL(currently in LAPACK) usingAR.ARis zero, then ifAIis in the safe range, callSCALwith the numberCMPLX( ZERO, -ONE / AI ). Otherwise, do the proper scaling by a power of two similarly toCSRSCL.Aare nonzero, then we computeURandUIsuch thatCMPLX( ONE / UR, - ONE / UI )is the reciprocal ofA, i.e.,URandUIare always different from zero. NaNs only appear if either:ARorAIis a NaN.ARandAIare both infinite, in which case it makes sense to propagate a NaN.Now, we have three cases:
URandUIare both in the safe range, then callSCALwith the numberCMPLX( ONE / UR, - ONE / UI ).URorUIare smaller thanSFMIN, it means that bothURandUIare small numbers. In fact, we can prove that they should be both smaller thanSFMIN(1+1/EPSILON)Therefore, it makes sense to scalexbyCMPLX( SFMIN / UR, - SFMIN / UI ), and then scale the result byONE / SFMIN. UseSCALto do both scalings.URorUIare greater thanONE / SFMIN, then we must check a few things:ARorAIis infinite, thenURandUIare either both infinite or both NaNs. Therefore, we scale byCMPLX( ONE / UR, - ONE / UI )usingSCAL.URorUIis infinite andARandAIare finite, then the algorithm generated infinite numbers that can be avoided. In this case, recompute scaled versions ofURandUIusingSFMIN. Then, scalexbySFMIN, and then scale the result byCMPLX( ONE / UR, - ONE / UI ). UseSCALto do both scalings.xbySFMIN, and then scale the result byCMPLX( ONE / (SFMIN*UR), - ONE / (SFMIN*UI) ). UseSCALto do both scalings.Some comments about the algorithm:
ARandAIto the result. Moreover, it propagates NaNs ifARandAIare both infinite. It does not generate NaNs in other cases.Ais zero.A. This is done by treating cases whereURorUIare infinite numbers.A.