Re: Sparse matrix problem

[John W. Eaton]

On  3-Jan-2007, David Bateman wrote:

| 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;
| }

Would someone who has time please report these problems to the GCC
maintainers and follow up until the problem is fixed?

Thanks,

jwe