Re: Sparse matrix problem

[David Bateman]

John W. Eaton wrote:
> 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
>
>   
The bug for the C99 gcc compiler for

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=24581

seems related to this though is only for mixed real/complex args..
Though as the attached code show the problem equally exists in C99
c-code as well..

D.


-- 
David Bateman                                address@hidden
Motorola Labs - Paris                        +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 6 72 01 06 33 (Mob) 
91193 Gif-Sur-Yvette FRANCE                  +33 1 69 35 77 01 (Fax) 

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary

#include <complex.h>
#include <stdio.h>

int main (void)
{
  double _Complex a = 1. + I;
  double _Complex b = 0.;
  
  printf("(1+i) / (0 + 0i) : %e %e\n", creal(a/b), cimag (a/b));
  printf("(1+i) / 0 : %e %e\n", creal(a/0.), cimag(a/0.));
  printf("1 / 0: %e\n", 1. / 0.);
  return 0;
}