This is a quick programming exercise to correct the following problem in Maple:

for n from 8 to 12 do 
  A := Matrix(2^n, 2^n, storage='sparse'):  # zero matrix
  print( 2^n, time(LinearAlgebra:-Transpose(A)) );
end do:

The problem is that the LinearAlgebra:-Transpose command is not sparse. That is, the time it takes is proportional to the overall size of the matrix and not to the number of non-zero entries - even when your matrix uses sparse storage. In this post we will look at what is required to program a new Transpose command which can handle much larger matrices.

The first thing we need is to get information out of the matrix. Run the following commands:

A := Matrix([[1,0,3,0],[2,0,0,1],[0,0,0,1]],storage=sparse);
op(A);

and you should see a sequence: 3, 4, {(1, 1) = 1, (2, 1) = 2, (1, 3) = 3, (2, 4) = 1, (3, 4) = 1}, etc. This is misleading: op(1,A) gives the sequence of dimensions (3,4) and op(2,A) gives you the set of elements {(1,1)=1, ...}. You can read all about it on the help page ?Matrix.

We are not going to worry too much about the shape and indexing functions; our function will work only for sparse rectangular matrices. To transpose the matrix we will need to go through the list of elements (a,b)=c and return (b,a)=c. We start by writing a small procedure:

reverselhs := proc(eqn)
  local a,b,c;
   a,b := lhs(eqn);
   c := rhs(eqn);
   (b,a) = c;
end proc:

The problem with this is that it is slow. We could eliminate a line and the local variable c by using rhs(eqn) in the last line, but we can get an even better improvement if we can use an inline procedure. Given an equation eqn := (a,b)=c, op(eqn) returns the sequence a,b,c. Thus we can write an inline procedure which takes a,b, and c as input and outputs (b,a)=c, and compose it with the op command when mapping it onto the set. This is what I mean:

reverse := proc(a,b,c) option inline; (b,a)=c end proc:
map(reverse@op, op(2,A));
op(2,A);  # compare

The only thing left is to quickly build the sparse matrix. We can use the table initialization option of the Matrix constructor to make this quick:

S := map(reverse@op, op(2,A));
Matrix(4,3,table(S));

Here is the completed command. Don't forget the 'storage'='sparse' option.

reverse := proc(a,b,c) option inline; (b,a)=c end proc:

transpose := proc(A::Matrix)
  local n, m, S;

   n, m := op(1,A);
   S := map(reverse@op, op(2,A));
   Matrix(m, n ,table(S), 'storage'='sparse');
end proc:

And now a little something for the experts out there:

reverse := proc(a,b,c) option inline; (b,a)=c end proc:
transpose := proc(A) option inline:
  Matrix(op(1,A)[2], op(1,A)[1], map(reverse@op, op(2,A)), 'storage'='sparse')
end proc:

You can make it work for Vectors too (keeping it inline). Have fun :)