New sparse class

classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

New sparse class

David Bateman-3
Dear All,

I've been working on rewriting Andy Adler's sparse class from octave-forge
recently, in a way that would be acceptable for inclusion in octave. That
is following the existing conventions and without the use of malloc/free.
John and Andy have seen the code and John is considering the inclusion of
this code in the core of octave.

At the moment this code is a seperate package to octave, rather than a
patch, due to the fact that this makes my compile time more reasonable
during the initial development. To allow joint development of this patch
John has placed in on the octave CVS server. For those that are interested
the code can be found at

http://www.octave.org/cgi-bin/viewcvs.cgi/octave-sparse/

or through the CVS server by doing

cvs -d :ext:[hidden email]:/cvs -z 9 checkout octave-sparse

At this point the basic functionality is complete and tested, but there
are still some missing features, these being

 * There are no solvers included yet, so the "/" and "\" operators do
   not work. This is a matter of adding the hooks to the external solver
   code, probably UMFPACK

 * I've added some functions using "dispatch" like reshape, diag, min, max,
   atan2. However the basic 1-arg mapper functions like real, abs, etc are
   handled by the octave class octave_mapper, and rather than use dispatch
   for these 51 functions, it would be better to fix the octave_mapper class
   in octave itself so that the mapper functions are directly available to
   external classes. Therefore the mapper functions are not yet available

However, I've attempted to be more compatiable with Matlab than Andy's
code which might lead to some wierd issues. For example the return type
of some operators will be different than Andy's code. Also the "./" operator
for two sparse matrices doesn't include the NaN values as implied by the
division by the zero elements of the sparse matrix. This is how Matlab
treats this even though it is inconsistent with full matrices.

One place where this code goes beyond Andy's code is the number of mixed
operators and the assign operator. There are many more mixed operators
than Andy's code, and this is now pretty much complete. As for the assignment
operator, things like

a = speye(10)
a(1) = 2
a(1,:) = 1
a(1,:) = rand(1,10)
a(1,:) = []

etc, now work as you would expect, rather than converting to a full
matrix first.

If you would like to look at the code, see the CVS server above, and please
feel free to comment, or even add stuff if you are motivated :-)

Regards
David

--
David Bateman                                [hidden email]
Motorola CRM                                 +33 1 69 35 48 04 (Ph)
Parc Les Algorithmes, Commune de St Aubin    +33 1 69 35 77 01 (Fax)
91193 Gif-Sur-Yvette FRANCE

The information contained in this communication has been classified as:

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


Reply | Threaded
Open this post in threaded view
|

Re: New sparse class

Paul Thomas-10

>Also the "./" operator
> for two sparse matrices doesn't include the NaN values as implied by the
> division by the zero elements of the sparse matrix. This is how Matlab
> treats this even though it is inconsistent with full matrices.
>

R12 seems to do the right thing

>> b = sparse( [1 2] , [2 1] , [5 0] , 2 , 2)
b =
  (1,2)        5
>> a = sparse( [1 2] , [2 1] , [5 5] , 2 , 2)
a =

   (2,1)        5
   (1,2)        5
>>  a ./ b
Warning: Divide by zero.
ans =
   (1,1)      NaN
   (2,1)      Inf
   (1,2)        1
   (2,2)      NaN
>> full(a) ./ full(b)
Warning: Divide by zero.
ans =
   NaN     1
   Inf   NaN

Paul T


Reply | Threaded
Open this post in threaded view
|

Re: New sparse class

David Bateman-3
Dapr├Ęs Paul Thomas <[hidden email]> (le 02/11/2004):

>
> >Also the "./" operator
> > for two sparse matrices doesn't include the NaN values as implied by the
> > division by the zero elements of the sparse matrix. This is how Matlab
> > treats this even though it is inconsistent with full matrices.
> >
>
> R12 seems to do the right thing
>
> >> b = sparse( [1 2] , [2 1] , [5 0] , 2 , 2)
> b =
>   (1,2)        5
> >> a = sparse( [1 2] , [2 1] , [5 5] , 2 , 2)
> a =
>
>    (2,1)        5
>    (1,2)        5
> >>  a ./ b
> Warning: Divide by zero.
> ans =
>    (1,1)      NaN
>    (2,1)      Inf
>    (1,2)        1
>    (2,2)      NaN
> >> full(a) ./ full(b)
> Warning: Divide by zero.
> ans =
>    NaN     1
>    Inf   NaN

Err, ok, then its a bug on my part... In fact the operation I meant was
".^", which Matlab R12 when used with a matrix and a scalar gives.

>> a = speye(2) .^ 0      

a =

   (1,1)        1
   (2,2)        1

>> eye(2) .^ 0

ans =

     1     1
     1     1


There are quite a few strange and undocumented misfeatures of this kind
in the Matlab sparse class.

Regards
David

--
David Bateman                                [hidden email]
Motorola CRM                                 +33 1 69 35 48 04 (Ph)
Parc Les Algorithmes, Commune de St Aubin    +33 1 69 35 77 01 (Fax)
91193 Gif-Sur-Yvette FRANCE

The information contained in this communication has been classified as:

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


Reply | Threaded
Open this post in threaded view
|

Re: New sparse class

Paul Thomas-10
Uuugggghhh!

> >> a = speye(2) .^ 0      
>
> a =
>
>    (1,1)        1
>    (2,2)        1
>
It also misses speye(2).^-N

I don't think that we should emulate that!

Paul