Open VikasChidananda opened 1 year ago
amontoison
commented
1 year ago Hi @VikasChidananda!
symamd computes a fill-in reducing permutation but it's based on an heuristic. Do you always have the same permutation in Julia and in OCTAVE / MATLAB?
What is the function that you called in OCTAVE / MATLAB?
VikasChidananda
commented
1 year ago 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)
The equivalent function in OCTAVE is also called symamd
amontoison
commented
1 year ago In AMD.jl, we directly call the C routines of SuiteSparse.
For MATLAB / OCTAVE they probably use the *.m files of the this folder.
Maybe the default options are different.
Does the fill-in of the factors is different after a factorization?
VikasChidananda
commented
1 year ago Sorry for the late reply. I realized the above is not the best representation of minimal working example. In fact yes, the fill are not comparable either! Consider the following code:
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: diag(Lv, -1) for both cases the fill-ins are different.
Julia
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-04Hi @VikasChidananda, can you post the number of non-zeros of the factor Lv? The goal of the permutation is to reduce the fill-in of the factors (sparsity of the factors). It's normal that we have different values in Lv because we have different permutation in Julia / Octave.
VikasChidananda
commented
1 year ago I understand. Thanks for clarifying. Here is the output for both:
julia> nnz(Lv)
241>> nnz(Lv)
ans = 242Ok, so the two permutations give almost the same fill-in. It's what in want in practice.
If you want to compare the reorderings, nnz is the most relevant metric.
dpo
commented
1 year ago 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 octave (e.g., by writing one of the matrices to file and reading it in at the other end)?
dpo
commented
1 year ago 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 stats array. It could help debug the issue.
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
VikasChidananda
commented
1 year ago 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 octave (e.g., by writing one of the matrices to file and reading it in at the other end)?
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.
dpo
commented
1 year ago @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!