Few Questions

[RoyiAvital]

Hello,
Your article and method looks really great.
I was wondering, when you mention MATLAB's direct solver, which function do you mean?
If you mean the mldivide, what was the memory configuration of your computer?

[dilipkay]

Hi - thanks for your interest. Yes, I did mean the mldivide. I can't remember the full details of the machine I used (it was during my PhD) but I think it was a 64-bit quad-core Xeon 3 GHz with 32GB of memory (on Redhat Linux).

[RoyiAvital]

I see.
Because it doesn't seem a normal computer with 4 / 8 GB of memory would be able to handle matrices of the size you use on the article.

Unless I'm missing something.

Do you think today the fastest solver for this kind of problems would be Conjugate Gradient with your Preconditioning?

Thank You.

[dilipkay]

Yes, 4 to 8GB was enough for our method but not for mldivide. "eigs" would run out of memory for smaller problems than for system solves. Even my 32GB machine ran out of memory with the very largest problems (for mldivide).

Regarding speed, no - our method is not the fastest for all nature of Laplacians. If you have a homogenous Poisson problem, there are much faster ways of solving it (e.g. using FFT's). But our method is fast, and is generally applicable to a range of Laplacians. But if you were to use our preconditioner, then using it with Conjugate Gradient would be the way to do it.

[RoyiAvital]

Hi,
I'm talking on inhomogeneous Laplacian matrices as you described in the article.
I would like to have a method which can handle very large scale matrices in a regular memory configurations (4 - 8 GB in total, 2.5 - 6.5 free) and be effective.

What do you think?

[dilipkay]

Then you have come to the right place ;-). I have tested on matrices upto 14 million in unknowns, and it fit into that kind of memory. By "very large scale" what size were you thinking of?

[RoyiAvital]

Hi,
Thought about something like 40 millions unknowns.

Do you have interface to any C based PCG solver?

[dilipkay]

Hi - sorry for the delayed reply. I had a deadline.

40 M unknowns is a lot. How many entries in your matrix? This will decide whether it is possible to fit into 8GB. For this work I was just using MATLAB's PCG solver. But it's very easy to write a PCG solver - you just need to have some dot products and that's pretty much it (assuming you have the call to the preconditioner available as a separate piece of C code).

[RoyiAvital]

@dilipkay ,

It seems your code generate the following error:

error C2733: 'mxUnshareArray': second C linkage of overloaded function not allowed

I use MATLAB R2018a on Windows (Visual Studio 2017).
It used to work on previous versions of MATLAB.

[99991]

I had the same issue when using octave. Maybe this fix also works in MATLAB.

EDIT(easier fix):
Replace the line mxUnshareArray(A, true); in /hsc/mex_funs/hsc_sparsify_and_compensate.cpp with A = mxDuplicateArray(A);

I am using octave 4.2.2 installed with apt because version 5.2.2 installed with snap is currently broken (/snap/octave/22/usr/include/x86_64-linux-gnu/c++/7/bits/os_defines.h:39:10: fatal error: features.h: No such file or directory).

This is the output of running demo_2D.m.

N = 65536; Coarse = 32768; Fine = 32768; sparsify and compensate 0.018711s
Level 1 coarse variables = 50.000000%

Computing prolongation matrix 0.001443 s
N = 32768; Coarse = 16384; Fine = 16384; sparsify and compensate 0.008693s
Level 2 coarse variables = 50.000000%

Computing prolongation matrix 0.000626 s
N = 16384; Coarse = 8192; Fine = 8192; sparsify and compensate 0.004913s
Level 3 coarse variables = 50.000000%

Computing prolongation matrix 0.000367 s
N = 8192; Coarse = 4096; Fine = 4096; sparsify and compensate 0.002004s
Level 4 coarse variables = 50.000000%

Computing prolongation matrix 0.000106 s
N = 4096; Coarse = 2048; Fine = 2048; sparsify and compensate 0.001117s
Level 5 coarse variables = 50.000000%

Computing prolongation matrix 0.000074 s
N = 2048; Coarse = 1024; Fine = 1024; sparsify and compensate 0.000504s
Level 6 coarse variables = 50.000000%

Computing prolongation matrix 0.000022 s
N = 1024; Coarse = 512; Fine = 512; sparsify and compensate 0.000310s
Level 7 coarse variables = 50.000000%

Computing prolongation matrix 0.000011 s
pcg: converged in 4 iterations. pcg: the initial residual norm was reduced 2.83409e+06 times.
pcg: converged in 4 iterations. pcg: the initial residual norm was reduced 2.24464e+06 times.
Setting up HSC
Setting threshold to Infininity !!
Level 1 unknowns 65536 nonzeros 586756 avg. bw 8.953186
Level 2 unknowns 32768 nonzeros 292866 avg. bw 8.937561
Level 3 unknowns 16384 nonzeros 145924 avg. bw 8.906494
Level 4 unknowns 8192 nonzeros 72706 avg. bw 8.875244
Level 5 unknowns 4096 nonzeros 36100 avg. bw 8.813477
Level 6 unknowns 2048 nonzeros 17922 avg. bw 8.750977
Level 7 unknowns 1024 nonzeros 8836 avg. bw 8.628906

HSC hierarchy setup in 0.223915 s.; Galerkin 0.056541s; Triangle sp. 0.036252s
Size of smallest level 512
HSC solver setup time taken 0.226021 
Solve time 0.0674973
Computing Condition Number...
Original Condition Number 1150710.135466; preconditioned system 2.692451

I am not sure if everything is working correctly. On the one hand, pcg converges in just 4 iterations, which is pretty good, but on the other hand, the output contains the line Setting threshold to Infininity !!. Is that an error? At least it seems to be important as suggested by the two exclamation marks.

Also I am not sure how to interpret the results of demo_2D.m. The doc string says that it performs colorization, but the result vector x is of size 65536 2. I would have expected something with 3 channels (RGB) instead. The first layer reshaped to 256 256 looks like this:

[HTDerekLiu]

Thanks for this super interesting paper. I encountered the same issue with mxUnshareArray. I tried to compile it with R2020a, and I got

path/hsc/mex_funs/hsc_sparsify_and_compensate.cpp:29:17: error: functions that differ only in their
return type cannot be overloaded
extern "C" bool mxUnshareArray(mxArray *array_ptr, bool noDeepCopy);
           ~~~~ ^
/Applications/MATLAB_R2020a.app/extern/include/matrix.h:1173:41: note: previous declaration is here
LIBMMWMATRIX_PUBLISHED_API_EXTERN_C int mxUnshareArray(mxArray *pa, int level);

It seems like mxUnshareArray has a different signature now. It seems like MATLABR 2020a needs a different solution from the one above. Could you give me some guidance about how I can fix it?

Best,

[bob1321-eng]

In ''hsc_code/mex_funs/hsc_sparsify_and_compensate.cpp'',the following two parts need to be modified：
1.The function mxunsharearray is no longer overloaded，Revise extern "C" bool mxUnshareArray(mxArray *array_ptr, bool noDeepCopy) to extern "C" int mxUnshareArray(mxArray *array_ptr, int noDeepCopy);
2.Replace mxUnshareArray(A, true) with A = mxDuplicateArray(A);