A new dense p-adic lift algorithm has improved the efficiency of LinearAlgebra[LinearSolve] and LinearAlgebra[Modular][IntegerLinearSolve]. See Enhancements to LinearAlgebra for details.

Maple 13 includes a fast modular algorithm for solving linear systems over cyclotomic fields. For example, if a linear system involves $\mathrm{RootOf}\left(f\,x\right)$ where $f\left(x\right)$ is a cyclotomic polynomial of degree > 1, then the fast code will be used. See also numtheory[cyclotomic].

As a comparison, consider the following example:

r := RootOf(x^4+x^3+x^2+x+1);

$r≔\mathrm{RootOf}\left({\mathrm{\_Z}}^4+{\mathrm{\_Z}}^3+{\mathrm{\_Z}}^2+\mathrm{\_Z}+1\right)$

cf := ()->randpoly(r, degree = 3):

nv := 50:

vars := seq(x[i],i=1..nv):

sys := {seq(randpoly([vars],degree=1,coeffs=cf),i=1..nv)}:

tt := time():

solve(sys,{vars}):

time()-tt;

$0.212$

In Maple 12, this same problem took almost 7 times as long.

with(LinearAlgebra):

N := 1000:
M := RandomMatrix(N,N,generator=-0.5..0.5, outputoptions=[datatype=float[8]]):

L := Matrix(M,shape=triangular[lower]):

M := Matrix(L.Transpose(L)):

for i from 1 to N do M[i,i]:=i; end do:

for i from 1 to N do M[i,i] := i; end do:

st := time():

LUDecomposition(M,method=Cholesky):

Warning, Matrix may be nonsymmetric; using only data from lower triangle

time() - st;

$0.02$

with(LinearAlgebra):

N := 1000:
M := RandomMatrix(N,N,generator=-0.5..0.5, outputoptions=[datatype=float[8]]):

st := time():

LUDecomposition(M):

time()-st;

$0.03$

with(LinearAlgebra):

m := 10000:
n := 100:
M := RandomMatrix(m,n,generator=-1.0..1.0, outputoptions=[datatype=float[8]]):

st := time():

SingularValues(M,output=[U,S,Vt],thin=true):

time()-st;

$1.78$

p1 := randpoly([x,y,z],degree=20,terms=50):

tt := time():

p2 := modp(Expand(p1 ^ 3),13):

time()-tt;

$0.013$

tt := time():

modp(Divide(p2,p1),13);

$\mathrm{true}$

time()-tt;

$0.014$

alias(p=RootOf(z^3+t*z^2+s*t,z)):

alias(q=RootOf(z^2+2)):

g := (y+p*t+q)*(t*x+s*y*p+q):

a := t*x+s*y+p*s*t+p^2+q:

b := x+s*x*y-t*p^2+p+t+q*s:

A := evala(Expand(g*a)):

B := evala(Expand(g*b)):

tt := time():

evala(Gcd(A, B));

time()-tt;

$0.5$

vars := [seq(x[i],i=0..19)];

$\mathrm{vars}≔\left[x_0\,x_1\,x_2\,x_3\,x_4\,x_5\,x_6\,x_7\,x_8\,x_9\,x_{10}\,x_{11}\,x_{12}\,x_{13}\,x_{14}\,x_{15}\,x_{16}\,x_{17}\,x_{18}\,x_{19}\right]$

f := add(i*vars[i], i=1..20);

$f≔x_0+2x_1+3x_2+4x_3+5x_4+6x_5+7x_6+8x_7+9x_8+10x_9+11x_{10}+12x_{11}+13x_{12}+14x_{13}+15x_{14}+16x_{15}+17x_{16}+18x_{17}+19x_{18}+20x_{19}$

g := expand(f^6):

(nops,length)(g);

$177100,10915266$

t := time():
(lcoeff,tcoeff)(g, vars);
time() - t;

$1,64000000$

$0.005$

p := proc(A::Vector,B::Vector,C::Vector,m::integer,n::integer)
    option hfloat;
    local i;

    for i from m to n
    do
        C[i] := A[i]^2+B[i]^2;
    end do;
end proc:
n := 10^7:
m := n/2:
A := LinearAlgebra:-RandomVector( n, outputoptions=[datatype=hfloat] ):
B := LinearAlgebra:-RandomVector( n, outputoptions=[datatype=hfloat] ):
C := Vector( 1..n, datatype=hfloat ):

tt := time[real]():
p(A,B,C,1,m):
p(A,B,C,m+1,n):
time[real]() - tt;

$33.089$

tt := time[real]():
id1 := Threads:-Create( p(A,B,C,1,m) ):
p(A,B,C,m+1,n):
Threads:-Wait( id1 ):
time[real]() - tt;

$22.873$

A := LinearAlgebra:-RandomMatrix(1000):

time(zip(`*`,A,2));

$0.608$

time( A *~ 2 );

$0.096$

B := LinearAlgebra:-RandomMatrix(1000,outputoptions=[datatype=float[8]]):

time(zip(`*`,B,2));

$1.684$

time( B *~ 2 );

$0.017$

time(map(`sin`,B));

$5.456$

time( sin~(B));

$0.116$

N := 500:

n := 100:

image := Array(1..N,1..N,datatype=float[8],order=C_order):

container := Array(1..n,1..n,datatype=float[8],order=C_order):

st := time():

for i from 1 to 1000 do
  ImageTools:-GetSubImage(image,10,10,n,n,output=container);
end do:

time() - st;

$0.220$

with(LinearAlgebra):

gcTest := proc()
          local A, n, resultList;

          for n to 30 do
              A := RandomMatrix(n,n);
              resultList:= LUDecomposition(A,output=['P','L','U1','R']);
          end do;
end proc:

kernelopts(gcfreq=10^6);

$\left[10000000\,0.1\right]$

gc():
start := time():
gcTest():
t1 := time() - start;

$\mathrm{t1}≔0.924$

kernelopts(gcfreq=10^7);

$1000000$

gc():
start := time():
gcTest():
t2:= time() - start;

$\mathrm{t2}≔0.860$