Hi, I'm the primary author of that file. There are two things being discussed here, the divide-by-zero warnings and the numerical errors. The same section of code is related to both:
olderr = np.geterr()['invalid'] # checks the current error handling > np.seterr(invalid='ignore') > K = np.array( > [C * sadp[t] * > (normalized_dp/(cp-cp[t]) - > (normalized_dp[t]/(cp-cp[t])).conjugate()) > for t in np.arange(NB)], dtype=np.complex128) > np.seterr(invalid=olderr) # resets the error handling > for i in xrange(NB): > K[i, i] = 1 > phi = np.linalg.solve(K, g) / NB * TWOPI # Nystrom The divide by zero warnings come from a large element-wise matrix divide. Only the off-diagonal entries of that division are necessary; the diagonal entries are explicitly overwritten in the following step. Some of the diagonal entries are divided by zero during the divide, producing a nan, but the overwrite means those entries don't affect the computation. The code is written to divide the whole matrix for simplicity and speed. The divide is wrapped in flags that are intended to temporarily disable the numpy divide-by-zero warnings. However, it looks like the numpy people change the syntax of those flags fairly often. If we prefer, the division could be changed to either skip the diagonal entries (inefficient) or fix the zero divisors (since they don't affect the result anyway). As far as I can tell, the numerical errors are completely independent of the divide-by-zero warnings. Given that we have determined ATLAS to be the culprit, the problem must be connected to the np.linalg.solve line. To the best of my knowledge that is the only section of code that uses those linear algebra libraries. It's nothing fancy, just a big complex-valued matrix solve. I don't know of any aspects that could be causing those errors. Most of the errors Jeroen showed are *far* beyond numerical fuzz. The algorithm is not *that* sensitive to numerical instabilities--at least not in the regions being tested by the doctests. It's also worth noting that there's a ticket (#11273) that has been in review for roughly two years (most of the delay is my fault) that adds a lot of new functionality to this package. It doesn't change any code in this section, although it does add more documentation saying roughly what I said above about the divide warnings. Let me know if there's anything I can do to help. I've subscribed to this thread so I'll see any responses. Ethan On Friday, January 18, 2013 2:31:18 AM UTC-8, François wrote: > > On 18/01/13 23:25, Volker Braun wrote: > > On Friday, January 18, 2013 9:47:53 AM UTC, Fran�ois wrote: > > > > doctest:1: RuntimeWarning: divide by zero encountered in divide > > It doesn't correlate to any of the error you reported but it adds to > > the feeling that something needs fixing in that file. > > > > > > I would say it is precisely the kind of mistake that one would suspect: > > Somewhere, a numerically unstable division by a very small value occurs. > > This greatly magnifies numerical errors from ATLAS routines. Can you do > > reproduce this in a verbose doctests and tell us where it happens? > > > > Look at: > http://trac.sagemath.org/sage_trac/ticket/11334#comment:56 > > There are line numbers. > -- You received this message because you are subscribed to the Google Groups "sage-devel" group. To post to this group, send email to [email protected]. To unsubscribe from this group, send email to [email protected]. Visit this group at http://groups.google.com/group/sage-devel?hl=en.
