-
Notifications
You must be signed in to change notification settings - Fork 12
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
symamd is not comparable with OCTAVE/MATLAB #60
Comments
Hi @VikasChidananda! |
Hi @amontoison , yes! I just checked for 5 runs, it always gives the same permutation for the above code in OCTAVE and Julia (as in the code output above) |
In AMD.jl, we directly call the C routines of SuiteSparse. Does the fill-in of the factors is different after a factorization? |
Sorry for the late reply. Julia Code using LinearAlgebra
using SparseArrays
using AMD
Re = 1e2
H = 2; #Height of the container
L = 5; #Length of the container
nx = 10;
ny = 7;
dx = L/(nx-1); #Width of space step(x)
dy = H/(ny-1); #Width of space step(y)
dt = 0.01
e = ones(nx-2);
i = ones(ny-1);
Ax = spdiagm(-1 => e[1:end-1], 0 => -2*e, 1 => e[1:end-1]);
Ay = spdiagm(-1 => i[1:end-1], 0 => -2*i, 1 => i[1:end-1]);
Ay[1,1] = -3;
Ay[end,end]=-3;
A = kron(Ay/dy^2, sparse(I, nx-2, nx-2)) + kron(sparse(I, ny-1, ny-1), Ax/dx^2);
N = (nx-2)*(ny-1)
Dv = sparse(I, N, N) - dt*A/Re;
pv = symamd(Dv);
F = lu(Dv[pv, pv]);
Lv, Uv = F.L, F.U; OCTAVE Code Re = 1e2
H = 2; %Height of the container
L = 5; %Length of the container
nx = 10;
ny = 7;
dx = L/(nx-1); %Width of space step(x)
dy = H/(ny-1); %Width of space step(y)
dt = 0.01
e=ones(nx-2,1);
i=ones(ny-1,1);
Ax=spdiags(e*[1 -2 1],-1:1,nx-2,nx-2);
Ay=spdiags(i*[1 -2 1],-1:1,ny-1,ny-1);
Ay(1,1)=-3;
Ay(end,end)=-3;
A=kron(Ay/dy^2,speye(nx-2)) + kron(speye(ny-1),Ax/dx^2);
Dv=speye((nx-2)*(ny-1))-dt*A/Re;
pv=symamd(Dv);
[Lv Uv]=lu(Dv(pv,pv)); and if we check: julia> diag(Lv, -1)
47-element SparseVector{Float64, Int64} with 37 stored entries:
[1 ] = -0.000322815
[2 ] = -2.89751e-7
[3 ] = -0.000323105
[5 ] = -0.000896997
[6 ] = -2.89824e-7
[7 ] = -0.000896997
[8 ] = -2.90272e-7
[9 ] = -0.000323105
[10] = -2.90178e-7
[14] = -0.000322919
[15] = -2.89657e-7
[18] = -0.000322815
⋮
[34] = -5.80357e-7
[35] = -5.80544e-7
[36] = -1.04431e-7
[37] = -8.06051e-7
[38] = -5.79576e-7
[39] = -8.06051e-7
[40] = -5.80357e-7
[42] = -0.000323106
[43] = -5.80357e-7
[44] = -4.67796e-13
[45] = -5.80096e-7
[46] = -7.30028e-22
[47] = -7.23026e-10
OCTAVE >> diag(Lv, -1)
ans =
Compressed Column Sparse (rows = 47, cols = 1, nnz = 38 [81%])
(1, 1) -> -8.9700e-04
(2, 1) -> -2.8992e-07
(3, 1) -> -2.6006e-10
(4, 1) -> -2.9018e-07
(7, 1) -> -3.2292e-04
(8, 1) -> -2.8966e-07
(9, 1) -> -3.2321e-04
(10, 1) -> -3.2321e-04
(11, 1) -> -2.9018e-07
(14, 1) -> -3.2292e-04
(15, 1) -> -2.8966e-07
(18, 1) -> -8.9700e-04
(19, 1) -> -2.8992e-07
(20, 1) -> -2.6006e-10
(21, 1) -> -2.9018e-07
(22, 1) -> -8.9780e-04
(23, 1) -> -2.9018e-07
(26, 1) -> -3.2292e-04
(27, 1) -> -2.8966e-07
(28, 1) -> -3.2321e-04
(30, 1) -> -3.2292e-04
(31, 1) -> -8.9700e-04
(32, 1) -> -2.8992e-07
(33, 1) -> -8.9700e-04
(34, 1) -> -5.8036e-07
(35, 1) -> -3.2321e-04
(36, 1) -> -2.8137e-10
(37, 1) -> -5.8036e-07
(38, 1) -> -5.8010e-07
(39, 1) -> -9.3454e-13
(40, 1) -> -5.8036e-07
(41, 1) -> -5.8010e-07
(42, 1) -> -4.6717e-13
(43, 1) -> -5.8036e-07
(44, 1) -> -1.0446e-07
(45, 1) -> -1.2125e-13
(46, 1) -> -5.8036e-07
(47, 1) -> -8.9780e-04 |
Hi @VikasChidananda, can you post the number of non-zeros of the factor |
I understand. Thanks for clarifying. julia> nnz(Lv)
241 >> nnz(Lv)
ans = 242 |
Ok, so the two permutations give almost the same fill-in. It's what in want in practice. |
I would expect the permutations to be identical though. SYMAMD is deterministic as far as I know and only depends on the sparsity pattern. @VikasChidananda Have you confirmed that the latter is identical in Julia and |
One possible difference is that we don't set default options explicitly:
Another is these options: https://github.com/DrTimothyAldenDavis/SuiteSparse/blob/dev/CCOLAMD/MATLAB/csymamdmex.c#L99 Also we don't return the finally, I notice that we don't call the "report" function in case of failure: https://github.com/DrTimothyAldenDavis/SuiteSparse/blob/dev/CCOLAMD/MATLAB/csymamdmex.c#L163 |
symamd
is not comparable with OCTAVE/MATLAB
symamd
is not comparable with OCTAVE/MATLAB
Yes, I checked the sparsity pattern of Matrices from both. They are identical, with slight difference being OCTAVE stores the rounded values (until 5 decimal points) and Julia stores until (20 decimal points). The permutations of both matrices (let's call it octave_Dv, Julia_Dv) are the same in Julia. It differs from that of OCTAVE. |
@VikasChidananda I wonder if #61 fixes your issue. |
I have been trying to convert some code from MATLAB to Julia. But I observe that symamd operation is not very comparable to MATLAB's.
I am pasting the whole code snippet for comparison:
Julia Code
Output in Julia:
Output in OCTAVE:
Any help/suggestion/insight is greatly appreciated,
Thanks!
The text was updated successfully, but these errors were encountered: