Thanks for the suggestion! My matrices originate from graph-theoretic NP-Complete problems such as 3-colorability, independent set, etc.
In particular, for 3-colorability, the matrices contain only 0, -1, or +1. Furthermore, the solution vector is all zeros, and one 1.
Every column contains either 2 or 3 entries. If 2 entries, then the column contains both +1 and -1. If 3 entries, then 3 +1s. I'm not sure how to express the number of entries per row, but its very sparse.
I would be particularly interested in pre-processing algorithms for simplifying this kind of matrix. Perhaps something special with the rank? I would love insight or comments, if anyone has the time.
Best,
Susan

Thank you very much! I had tried LinBox before with NO success, but I see that they have a new release out, as of January 31, 2007! And so this is very exiciting news.
I will certainly give it a try.
Best,
Susan

Hello, all! Thank you very much for your earlier comments! To acer: no, it doesn't work with float[8]. And to roman, I have actually found a very small example, just a 6 x 9, matrix, that illustrates the problem. A maple script of the bad test case follows below. Any insight would be great AND fabulous!
Take care,
Susan
restart;
with(LinearAlgebra);
ls_m := Matrix(6,9 , datatype='integer[1]', storage='sparse');
ls_m[1,1] := 1;
ls_m[2,2] := 1;
ls_m[2,4] := 1;
ls_m[3,5] := 1;
ls_m[4,1] := -1;
ls_m[4,3] := 1;
ls_m[4,4] := -1;
ls_m[4,7] := 1;
ls_m[5,2] := -1;
ls_m[5,5] := -1;
ls_m[5,6] := 1;
ls_m[5,8] := 1;
ls_m[6,3] := -1;
ls_m[6,6] := -1;
ls_m[6,9] := 1;
B := Vector(6,0);
B[6] := 1;
# Gives error:
# Error, (in LinearAlgebra:-LA_Main:-LinearSolve) method
# SparseIterative incompatible with data of type anything
sol := LinearSolve(ls_m, B, method='SparseIterative', free='t');
# Gives error:
# Error, (in LinearAlgebra:-LA_Main:-LinearSolve) method SparseLU
# incompatible with data of type anything
sol := LinearSolve(ls_m, B, method='SparseDirect', free='t');
# works perfectly
sol := LinearSolve(ls_m, B, free='t');

Hello, all! Thank you very much for your earlier comments! To acer: no, it doesn't work with float[8]. And to roman, I have actually found a very small example, just a 6 x 9, matrix, that illustrates the problem. A maple script of the bad test case follows below. Any insight would be great AND fabulous!
Take care,
Susan
restart;
with(LinearAlgebra);
ls_m := Matrix(6,9 , datatype='integer[1]', storage='sparse');
ls_m[1,1] := 1;
ls_m[2,2] := 1;
ls_m[2,4] := 1;
ls_m[3,5] := 1;
ls_m[4,1] := -1;
ls_m[4,3] := 1;
ls_m[4,4] := -1;
ls_m[4,7] := 1;
ls_m[5,2] := -1;
ls_m[5,5] := -1;
ls_m[5,6] := 1;
ls_m[5,8] := 1;
ls_m[6,3] := -1;
ls_m[6,6] := -1;
ls_m[6,9] := 1;
B := Vector(6,0);
B[6] := 1;
# Gives error:
# Error, (in LinearAlgebra:-LA_Main:-LinearSolve) method
# SparseIterative incompatible with data of type anything
sol := LinearSolve(ls_m, B, method='SparseIterative', free='t');
# Gives error:
# Error, (in LinearAlgebra:-LA_Main:-LinearSolve) method SparseLU
# incompatible with data of type anything
sol := LinearSolve(ls_m, B, method='SparseDirect', free='t');
# works perfectly
sol := LinearSolve(ls_m, B, free='t');

Just great! Wonderful! My research would grind to a complete halt if it wasn't for this website and the people who respond to posts!! :) :)
Thank you. This exactly did the trick.

Just great! Wonderful! My research would grind to a complete halt if it wasn't for this website and the people who respond to posts!! :) :)
Thank you. This exactly did the trick.

Hello! That solution is perfectly fabulous! But unfortuneatly, my polynomial is of the form:
mon:=c[1]^2+c[1]*c[2]+c[1]*c[3]+c[2]^2+c[2]*c[3]+c[3]+c[1]+c[2]+c[3]+1;
poly := mon*(v[1]^2 + v[2]) + mon*(v[2] + v[2])^2 + mon*(v[1] + v[2]^2);
I would really like the coeffs of all of the v's.
I can do:
C := [coeffs(g, {v[1],v[2]}, 'M')];
But I can't do:
C := [coeffs(g, {v}, 'M')];
Since I don't know that number of v's in advance, I could use a little help!
Let me know what you think if you have a moment,
Susan

Hello! That solution is perfectly fabulous! But unfortuneatly, my polynomial is of the form:
mon:=c[1]^2+c[1]*c[2]+c[1]*c[3]+c[2]^2+c[2]*c[3]+c[3]+c[1]+c[2]+c[3]+1;
poly := mon*(v[1]^2 + v[2]) + mon*(v[2] + v[2])^2 + mon*(v[1] + v[2]^2);
I would really like the coeffs of all of the v's.
I can do:
C := [coeffs(g, {v[1],v[2]}, 'M')];
But I can't do:
C := [coeffs(g, {v}, 'M')];
Since I don't know that number of v's in advance, I could use a little help!
Let me know what you think if you have a moment,
Susan

I really, really appreciated your comments, and I tried all sorts of things regarding the free variables, adjusting the number of Digits, trying float, etc., but alas! I am still not fixed.
Please consider the following:
my_solve := proc(A::Matrix)
local sz, local_A, B, sol;
sz := Dimension(A);
local_A := Matrix(A, datatype=float);
B := Vector(1...sz[1], 1);
sol := LinearSolve(local_A, B);
check := local_A[1,1...sz[2]] . sol;
print(check);
end proc;
Clearly, check should equal 1, because B is a vector of 1s. But you won't, even believe this, but check is coming in at 5 x 10^8!!!!! Or something. The point is, it is not even remotely close in any sense of the word to 1! I *must* be doing something very wrong. Taking the wrong dot product maybe? Maybe sol isn't what I think it is? I tried substituting 1s for the free variables, but still had problems.
I *know* I am doing something silly. Rounding could never account for an 8 orders of magnitude error!
Any comments will be awaited most eagerly.
PS> The matrix I am using is just too big to include here. But it has real entries between 0 and 1.
random_matrix := RandomMatrix(S_card, n, generator=rand(0..1000)/1000.0);
Thank you very much!

I really, really appreciated your comments, and I tried all sorts of things regarding the free variables, adjusting the number of Digits, trying float, etc., but alas! I am still not fixed.
Please consider the following:
my_solve := proc(A::Matrix)
local sz, local_A, B, sol;
sz := Dimension(A);
local_A := Matrix(A, datatype=float);
B := Vector(1...sz[1], 1);
sol := LinearSolve(local_A, B);
check := local_A[1,1...sz[2]] . sol;
print(check);
end proc;
Clearly, check should equal 1, because B is a vector of 1s. But you won't, even believe this, but check is coming in at 5 x 10^8!!!!! Or something. The point is, it is not even remotely close in any sense of the word to 1! I *must* be doing something very wrong. Taking the wrong dot product maybe? Maybe sol isn't what I think it is? I tried substituting 1s for the free variables, but still had problems.
I *know* I am doing something silly. Rounding could never account for an 8 orders of magnitude error!
Any comments will be awaited most eagerly.
PS> The matrix I am using is just too big to include here. But it has real entries between 0 and 1.
random_matrix := RandomMatrix(S_card, n, generator=rand(0..1000)/1000.0);
Thank you very much!

Thank you so much!!! Now everything is working just perfectly, exactly the way I want it! On to the next thing! :) :)
Thank you very much,
Susan

Thank you so much!!! Now everything is working just perfectly, exactly the way I want it! On to the next thing! :) :)
Thank you very much,
Susan