Gonzalo Tornaria on Tue, 20 May 2003 12:00:08 -0500

 Programs: Sparse matrix multiplication // charpoly

```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.
```