Last active
July 6, 2023 22:44
-
-
Save silgon/671930ff2323d588d4ec71a8c442f9fd to your computer and use it in GitHub Desktop.
exploit_structure.jsl
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| 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)) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment