exploit_structure.jsl

exploit_structure.jl

using SparseArrays
using LinearAlgebra
using BandedMatrices
using SpecialMatrices
using ToeplitzMatrices
using Random
using BlockArrays
using FFTW
BLAS.set_num_threads(1)

n = 5000;
y = rand(n);

# Diagonal Matrix
a = rand(n);
A = diagm(a);
A_d = Diagonal(a);
Asp = spdiagm(a);

@time x = A\y;
@time x = A_d\y;
@time x = Asp\y;
@time x = y./a;

@benchmark x = A\y setup=(A=diagm(a))
@benchmark x = A_d\y setup=(A_d = Diagonal(a))
@benchmark x = Asp\y setup=(Asp = spdiagm(a))
@benchmark x = y./a setup=(a = rand(n))


# Banded Matrix
Abd = BandedMatrix(-1=> rand(n-1), 0=> rand(n), 1=>rand(n-1));
Asp = spdiagm(-1=> rand(n-1), 0=> rand(n), 1=>rand(n-1));
A = Array(Abd);
@time x = A\y;
@time x = Abd\y;
@time x = Asp\y;

@benchmark x = A\y setup=(A = Array(Abd)) seconds=300  # 5 minutes because the computation is slow
@benchmark x = Abd\y setup=(Abd = BandedMatrix(-1=> rand(n-1), 0=> rand(n), 1=>rand(n-1)))
@benchmark x = Asp\y setup=(Asp = spdiagm(-1=> rand(n-1), 0=> rand(n), 1=>rand(n-1)))


# Non sparse possible
# Circulant Matrix
Acr = Circulant(rand(n));
A = Array(Acr);
@time x = A\y;
@time x = Acr\y;

@benchmark x = A\y setup=(A = Array(Acr)) seconds=300  # 5 minutes because the computation is slow
@benchmark x = Acr\y setup=(Acr = Circulant(rand(n)))


# Vandermonde matrix
Av = Vandermonde(rand(n));
A = Array(Av);
@time x = A\y;
@time x = Av\y;


# Matrix with structure
n = 5000;
n_ext = Int(round(n*1.01));
y_1 = rand(n);
y_2 = rand(n_ext-n);
y = mortar([y_1,y_2]);
a = rand(n);

A = Diagonal(a);
B = rand(n, n_ext-n);
C = rand(n_ext-n,n);
D = rand(n_ext-n,n_ext-n);

AA = Array(mortar(reshape([A,C,B,D],2,2)));
y = Array(y);

AA_sp = sparse(AA);


@time x = AA\y;
@time x = AA_sp\y;
@time begin
tmp1 = y_2-C*(A\y_1);
tmp2 = D-C*(A\B);
x_2 = tmp2\tmp1;
x_1 = A\(y_1-B*x_2);
x = mortar([x_1,x_2]);
end


@benchmark AA\y setup=(AA=Array(mortar(reshape([A,C,B,D],2,2)))) seconds=120
@benchmark x=AA_sp\y setup=(AA_sp=sparse(AA))
@benchmark begin
tmp1 = y_2-C*(A\y_1);
tmp2 = D-C*(A\B);
x_2 = tmp2\tmp1;
x_1 = A\(y_1-B*x_2);
x = mortar([x_1,x_2]);
end setup=(A=Diagonal(a))

n = 1000;
y = rand(n);
A = rand(n,n);
luA=lu(A);
@benchmark x=A\y setup=(A=rand(n,n))
@benchmark lu(A) setup=(A=rand(n,n))
@benchmark luA\y setup=(luA=lu(A))