(this post is long than I thought it would be.

contents:

re: zero-copy reshape (and its myriad uses)

re: sparse matrices

re: matlab compatability

re: MPI and PVM

)

re: zero-copy reshape

>reshape as currently implemented requires a matrix copy,

>but that would be easy to fix.

I agree that a copy is inconvenient. Usually it is not an issue

since the other ops are more expensive than the copy.

But, one thing that I think could be really cool on this topic

is to abstract the block of doubles allocated to a matrix/vector

from the way the block is indexed, supporting multiple

views for the same data, e.g.:

B = A'

or

B = reshape(A,prod(size(A)),1)

would allocate no new memory for B, just create a new 'view'. A new

block only gets allocated when one writes to an entry of B. This

kind of view business should be compatable with the basic BLAS ops

since they allow transpose bits and strides to be set appropriately.

I noticed that GSL supports this trick, to the extent that

v = A(1,:)

is just another view. I'm not sure if they handle

B = A([1 3 4],:)

as a view, and I don't see how one could do a gemm on such a B

without forcing a tmp copy.

BTW I am not an advocate of GSL-ing octave. I've looked at some

of the GSL sources, and I'm not as impressed with the implementation

as I am with the interface.

>There is something to be said for convenient data

>structures, especially if it is convenient to operate

>with slices. The results are going to be more readable

IMHO I'd change 'especially' to 'only'. Sometimes the operations

we think we want aren't the right ones, and when we figure out what

ops we want, we figure out better representations for our data.

>than tricky indexing and reshaping operations. It's

>not clear to me that Matlab does the right thing --- not

>while squeeze/reshape are required to transform a

>slice using a matrix operation.

Definitely. I think it is very hard to work out what one

really should have basic ops on multidimensional data. One

ugly, but sensible, thing one should have is something like

C = CONTRACT(A,B,[2 3],[4 1])

where C(a,b,c) = SUM_{xy} A(a,x,y)*B(y,b,c,x)

but that's kind of nasty. On the other hand one also can do

C = reshape( reshape(A,na,nx*ny)*reshape( reshape(B,ny*nb*nc,nx)' ,nx*ny,nb*nc) ,na,nb,nc)

which is also kind of nasty, but at least it introduces no new ops

to the language.

And yes, you have to have the cock-eyed view of

matrices like I have , so you are right that this is for the

sophisticated user and not for everyone. Maybe being forced

to do things like this makes one sophisticated. Maybe I'd just

too proud of myself for tricks like this.

If zero-copy reshape/transpose were implemented, and some support

for inline functions were around, then perhaps one could just

write .m scripts to support N-D arrays and sparse matrices.

re: sparsity

For sparse matrices, I agree that idxop is definitely for sophisticates

and wrapping it up might require making some sort of special 'struct'

with some kind of overloading. I'm kind of partial to some special

OO struct type which has a special field which says 'I have functions

in me', and then when someone does A*B and the interpreter sees that one

of them is a struct, and instead of the usual 'I don't know what to do with

structs' error message, it can check for the OO flag and look in the struct

fields for 'times' to pull out the appropriate functions. This is a

nice way to take advantage of the fact that functions are first-class

(or nearly first class) objects in octave, and I think it is better than

MATLABs "greedy overloading by directory name" crap.

function s = sparse(A)

[i,j,a] = find(A);

s.OO = true;

s.i = i; s.j = j; s.a = a;

s.times = sparse_times;

s.plus = sparse_plus;

...

endfunction

BTW just one last thing on idxop. The thing that really necessitates it

is the fact that += doesn't do 'the right thing' (which may not be right

to some people) when the LHS is multiply indexed:

v = zeros(2,1);

v([1 2 2]) += 5

currently sets v = [5 5] instead of [5 10]. A sparse multiply could then

be

[i,j,a] = find(A);

y = zeros(size(A,1),1)

y(i) += a.*x(j)

I had this argument before with some people on this list (and I'm content tha

my position may not be the right one), and I am aware that it would require

some serious fiddling with octave sources to do it my way. Then there is

the issue of not having a max= or a min= operator for the other idxops,

though one might co-opt &= and |= for these services.

re: MATLAB compatability

The bit about rewriting code on the web is a really good one, for

which I don't have a good answer. MATLAB has allowed people to

write some obnoxious things, and for octave to support that code

it needs to support some of these obnoxious things. Perhaps, it

can be done with little change to octave other than some mods to

the basic operations then wrapped up in .m's. That would be nice.

My point is that octave presents an opportunity not just to be

compatable with MATLAB, but to go that extra mile and do things

smarter than MATLAB. (BTW another personal example of this is that

octave includes the LAPACK sylvester equation solver, which, to my

shock, MATLAB doesn't, and I was solving sylvester equations a lot

at one point and probably will need to again in the future). I'd

like to see new features put in because it is a good idea, whether

from MATLAB or elsewhere, with a MATLAB-like interface when

applicable.

re: MPI and PVM

>Has anyone looked at the PVM project for scilab?

Three years ago I wrote some PVM DLD functions for octave. The mkoctfile

script has become way more handy and useful for these purposes (either

that or I have learned how to use mkoctfile). I had a lot of fun creating

matrices on one octave process and passing them to another octave process.

It was a fun toy, but it was a hack. If someone would like to unhack it

it is at

http://www.eskimo.com/~ripper/research/pvmoct.htmlI had thought of doing the same for MPI, but that whole ARGV thing made

me stop.

-r

-----Original Message-----

From: Paul Kienzle [mailto:

[hidden email]]

Sent: Saturday, March 01, 2003 7:32 AM

To: Lippert, Ross A.

Cc: octave-maintainers mailing list

Subject: Re: Octave 3.0

Lippert, Ross A. wrote:

>I think it is important to separate those items which are on this list

>for the sake of MATLAB compatibility and those items on this list which

>are for the general improvement of the package. For example, MATLAB's

>N-dimensional arrays are a convenient storage class, but they don't

>interface with the mathematical operations very well. Many times people

>think they want an N-dimensional array because they are dealing with,

>say a 3D mesh, but I think one quickly finds that it is more worthwhile

>to keep these values in a 1D vector and the does 'reshape' as necessary.

>

reshape as currently implemented requires a matrix copy,

but that would be easy to fix.

There is something to be said for convenient data

structures, especially if it is convenient to operate

with slices. The results are going to be more readable

than tricky indexing and reshaping operations. It's

not clear to me that Matlab does the right thing --- not

while squeeze/reshape are required to transform a

slice using a matrix operation.

>* Maybe I am just shooting my mouth off, but is there anyone out there who

>makes any serious use of N-d arrays in their MATLAB work?

>

>** On the other hand, if this is a move to further increase octave's .m

>compatability, then it is necessary, however awful.

>

>As for sparse matrices, I have gotten a lot of mileage out of a DLD function

>Paul Kienzle and I came up with called idxop which does the following:

> y = idxop(idx,x,['sum'|'prod'|'max'|'min'])

>where idx is an index vector whose length is equal to the number of rows of x

>and the resulting y is essentially given by

> y = zeros(max(idx),size(x,2))+[0|1|-inf|inf]

> for i=1:size(x,1),

> y(idx(i),:) = [plus|times|max|min] ( y(idx(i),:), x(i,:) )

> endfor

>This single function allows for a variety of sparse operations to occur on

>regular vectors/matrices without the need for a sparse data structure. The

>advantage of doing things in terms of sparse ops and not sparse data is that

>quite often the operations one wishes to perform are a mix of sparse and

>dense ops (e.g. you have a graph laplacian on a disconnected graph and you

>wish to form the dense laplacian for each graph component and take its

>eigenvalues). It is my belief that in practise, one spends quite a bit

>of time doing 'sparse' and 'full' when using MATLAB sparse matrix functions

>just so they can do a sparse op on their data at some point in the compute.

>

>Again * and ** apply.

>

* applies until you can provide recipes for

common 3-D and sparse operations that

unsophisticated users can apply --- yes you

can program everything in machine code but

macro assemblers sure help a lot. The octave

wiki is a good forum for doing this:

http://www.scarymath.org/octave** applies while there exists interesting

code on the web that you want to use

without having to rewrite.

>Private functions would be a nice way to clean up the ever growing octave

>namespace, whether MATLAB allows this currently or not.

>

Not enough I expect. To get that you need something

like namespace support.

>MPI support would raise a lot of eyebrows. It is something my PhD advisor

>was once trying to shoe-horn into a customized version of MATLAB, but with

>little support from the mathworks. Such a feature would really differentiate

>octave.

>

Has anyone looked at the PVM project for scilab?