restart:The float[8] technique used in this worksheet is extremely fast and is limited to moduli less than 2^25. If you need a larger modulus, it is very easy to modify (let me know), and just a little slower.N:= 2^7: #order of matrix A
n:= 2^12: #number of vectors b[i]
p:= 127: #prime modulusisprime(p);Generate a random example upper triangular nonsingular A.A:= LinearAlgebra:-RandomMatrix(
N, shape= triangular[upper, unit], datatype= float[8]
):Convert to shape= rectangular, which is required by LinearAlgebra:-Modular.A:= Matrix(A, shape= rectangular):Input A to the Modular system.A:= LinearAlgebra:-Modular:-Mod(p, A, float[8]):Generate a list of random side vectors b[k].b:= [seq](
LinearAlgebra:-RandomVector(N, datatype= float[8]),
k= 1..n
):Make a solving procedure that can be mapped into parallel processors. It returns the solution for a single side vector.Solve:= proc(b::Vector, A::Matrix, p::posint)
uses LAM= LinearAlgebra:-Modular;
local B:= LAM:-Mod(p, b, float[8]);
LAM:-BackwardSubstitute(p, A, B);
B
end proc: For reasons that I don't understand, it seems necessary to make a single bogus call to LinearAlgebra:-Modular:-BackwardSubstitute and generate an error before the parallel version will work. This only seems to be necessary once after each restart.LinearAlgebra:-Modular:-BackwardSubstitute(p);Threads:-Map is a parallel processing version of map. All that you need to run the parallel code is...# X:= Threads:-Map(Solve, b, A, p);...the rest of the code in the command below is just there to measure the time.CodeTools:-Usage(
proc() :-X:= Threads:-Map(Solve, b, A, p) end proc()
): It looks like the code parallelizes extremely well:36.*kernelopts(numcpus)/234.;It's actually the best parallelization ratio that I've ever achieved in my own code. But if you compare the time with the single-processor time for the same number of b's, you will see that there is no savings using the parallel code. You would have to use about 2^16 b's before there would be any savings.Verify results (not really necessary):for k to n do
if LinearAlgebra:-Norm(
LinearAlgebra:-Modular:-AddMultiple(
p,
p-1,
LinearAlgebra:-Modular:-Multiply(p, A, X[k]),
LinearAlgebra:-Modular:-Mod(p, b[k], float[8])
)
)
<> 0
then error "Bad solve"
end if
end do;
"All solutions check";