Programs: Sparse matrix multiplication // charpoly

[Gonzalo Tornaria <tornaria@mail.ma.utexas.edu>]

Here are the GP programs; they don't check the dimensions of the
matrices in any way, and in fact there's no way to know the dimensions
of a "studied matrix" (only one dimension is known).

The functions "matrix_sparseleft" and "matrix_sparseright" study any
matrix for multiplication on the left/right.

The functions "sparseleft_mul" and "sparseright_mul", multiply studied
matrices by regular ones. One could define multiplication like:

  mul_l(m1,m2) = sparseleft_mul(matrix_sparseleft(m1), m2);
  mul_r(m1,m2) = sparseright_mul(m1, matrix_sparseright(m2));

Lastly, the function "matrix_sparsecharpoly" computes the charpoly of
(any) matrix using this. It is correct in any case, but it should be
a factor worse than the PARI function for dense matrices. I believe the
factor is constant, and not too big (going to 1 as the size of the
matrix goes to infinity, and in practice seems to be less than 1.10 for
100x100).

Note that I have "inlined" the function "sparseright_mul", and replaced
the matrix that is computed in each iteration by a vector of vectors
(the columns).

The random sparse matrices I'm constructing are given by

DENSE(n)=matrix(n,n, i,j, (random(10000)+1));
SPARSE(n)=matrix(n,n, i,j, (random(10000)+1)*random(n)<4));
SPARSE1(n)=matrix(n,n, i,j, random(n)<4));

Gonzalo

/*******************************************************************************/

/* Studies a matrix to improve multiplication on the left */
matrix_sparseleft(m) =
  { local(n, k, row);
    n=matsize(m)[1];
    k=matsize(m)[2];
    vector(n, i,
      row=[];
      for(j=1,k,
         if(m[i,j]!=0, row=concat(row,[[j,m[i,j]]]))
      );
      row
    )
  }

/* Studies a matrix to improve multiplication on the right */
matrix_sparseright(m) =
  { local(n, k, col);
    n=matsize(m)[2];
    k=matsize(m)[1];
    vector(n, i,
      col=[];
      for(j=1,k,
         if(m[j,i]!=0, col=concat(col,[[j,m[j,i]]]))
      );
      col
    )
  }

/* Multiplies a "sparseleft" by a matrix */
sparseleft_mul(sl, m) =
  { local(n, k, v0, v);
    n=length(sl);
    k=matsize(m)[2];
    v0=vector(k,i,0);
    v=vector(n,i,
     sum(j=1,length(sl[i]), m[sl[i][j][1],]*sl[i][j][2], v0));
    matrix(n,k, i,j, v[i][j])
  }

/* Multiplies a matrix by a "sparseright" */
sparseright_mul(m, sr) =
  { local(n, k, v0, v);
    n=length(sr);
    k=matsize(m)[1];
    v0=vector(k,i,0)~;
    v=vector(n,i,
     sum(j=1,length(sr[i]), m[,sr[i][j][1]]*sr[i][j][2], v0));
    matrix(k,n, i,j, v[j][i])
  }

/* Computes the characteristic polynomial of a matrix */
matrix_sparsecharpoly(m) =
  { local(n, v0, Fk, pk, sr, sri);
    n=matsize(m)[1];
    v0=vector(n,i,0);
    Fk=vector(n,i,v0);
    pk=1; /* k=0 */
    sr=matrix_sparseright(m);
    Pol(concat([1],
    vector(n,k,
      /* Fk = Fk + pk */
      for(i=1,n, Fk[i][i]+=pk);
      /* Fk = Fk * m */
      Fk=vector(n,i,
        sri=sr[i];
	sum(j=1, length(sri), Fk[sri[j][1]]*sri[j][2], v0)
      );
      /* pk = -trace(Fk)/k */
      pk=-sum(i=1,n,Fk[i][i])/k;
    )));
  }

/*******************************************************************************/

-- 
GM/CS/S d? a-- C++(+++)  UL+++(++++)  P++>+++  L+++>++++ E---
W-(+) N+(++) w--- O M-(--) V-(--) PGP(--) b++(+++) G++ e>++++

... A mathematician is a machine for converting coffee into theorems.