symamd is not comparable with OCTAVE/MATLAB about amd.jl HOT 11 OPEN

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?

from amd.jl.

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

from amd.jl.

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?

from amd.jl.

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-04

from amd.jl.

Hi @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.

from amd.jl.

I understand. Thanks for clarifying.
Here is the output for both:

julia> nnz(Lv)
241

>> nnz(Lv)
ans = 242

from amd.jl.

Ok, 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.

from amd.jl.

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)?

from amd.jl.

One possible difference is that we don't set default options explicitly:

• https://github.com/DrTimothyAldenDavis/SuiteSparse/blob/dev/CCOLAMD/MATLAB/csymamdmex.c#L90 vs
• https://github.com/JuliaSmoothOptimizers/AMD.jl/blob/main/src/COLAMD.jl#L169
  (should be the same but we should double check).

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

from amd.jl.

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.

from amd.jl.

@VikasChidananda I wonder if #61 fixes your issue.

from amd.jl.