[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Sparse matrix problem
From: |
David Bateman |
Subject: |
Re: Sparse matrix problem |
Date: |
Wed, 03 Jan 2007 18:50:48 +0100 |
User-agent: |
Thunderbird 1.5.0.7 (X11/20060921) |
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;
}
- Re: Sparse matrix problem, David Bateman, 2007/01/03
- Re: Sparse matrix problem, Michael Goffioul, 2007/01/03
- Re: Re: Sparse matrix problem, Dmitri A. Sergatskov, 2007/01/04
- Re: Sparse matrix problem, David Bateman, 2007/01/04
- Re: Re: Sparse matrix problem, John W. Eaton, 2007/01/04
- Re: Re: Sparse matrix problem, Dmitri A. Sergatskov, 2007/01/04