RE: Re: Sparse matrix problem

[John W. Eaton]

On  3-Jan-2007, address@hidden wrote:

| > Ok, this was a lot more work than I though it would be, but here is the
| > patch. Basically it allows scalars to be stored in a sparse matrix
| > container and to be treated as if they were scalars. This is necessary
| > as Mathworks chose the path of letting the user store scalars as sparse
| > matrices rather than automatically reconverting these to scalars. In any
| > case the attached patch makes octave matlab compatible for scalars
| > stored as sparse matrices for all operators...
| 
| But it does not solve the fact that (note element-wise division):
| 
| [1+i 2-i] ./ [0 0-i]
| 
| gives under Octave
| 
| [NaN + NaNi     1 +   2i]
| 
| and under Matlab
| 
| [Inf +    Infi   1.0000 + 2.0000i]
| 
| Does it?

What happens when you compile and run the following program?

  #include <complex>
  #include <iostream>

  int
  main (void)
  {
    std::cout << std::complex<double> (1.0, 1.0) / 0.0 << std::endl;
    return 0;
  }

With g++ 4.1.2 on my system, it prints (nan,nan).  That seems like a
bug to me, but maybe it is standard conforming (and if it is, wtf?!?)?

If it is a bug, then the bug should be fixed in the compiler, not
Octave.

If it is standard conforming behavior, then we should find a way to
work around the problem in Octave.  Probably this will mean using our
own complex class.  Unfortunately, then the bug will only be fixed for
Octave code that uses the special class, and people using std::complex
in their own code (.oct or .mex files) will still have the buggy
behavior.

In the current <complex> header that is part of libstcd++ in the GCC
source tree, I find:

  // 26.2.5/15
  // XXX: This is a grammar school implementation.
  template<typename _Tp>
    template<typename _Up>
    complex<_Tp>&
    complex<_Tp>::operator/=(const complex<_Up>& __z)
    {
      const _Tp __r =  _M_real * __z.real() + _M_imag * __z.imag();
      const _Tp __n = std::norm(__z);
      _M_imag = (_M_imag * __z.real() - _M_real * __z.imag()) / __n;
      _M_real = __r / __n;
      return *this;
    }

  [...]

  //@{
  ///  Return new complex value @a x divided by @a y.
  template<typename _Tp>
    inline complex<_Tp>
    operator/(const complex<_Tp>& __x, const complex<_Tp>& __y)
    {
      complex<_Tp> __r = __x;
      __r /= __y;
      return __r;
    }
    
  template<typename _Tp>
    inline complex<_Tp>
    operator/(const complex<_Tp>& __x, const _Tp& __y)
    {
      complex<_Tp> __r = __x;
      __r /= __y;
      return __r;
    }

so maybe we need a better than grammar school implementation here.

jwe