amen_solve2 fails to converge, as opposed to its ttpy counterpart

[aiboyko]

So Ivan believes that it's a distinctive bug in MATLAB version

so on the same matrix, with the same inputs:

x=amen_solve2(A,b,tol,'max_full_size',15000,'kickrank',10,'x0',b,'nswp',15,...
'resid_damp',1 ,'rmax',100);

ttpy
amen_solve: swp=1, max_dx= 1.999E+00, max_res= 1.999E+00, max_rank=12
amen_solve: swp=2, max_dx= 1.939E+00, max_res= 4.579E-01, max_rank=22
amen_solve: swp=3, max_dx= 7.191E-01, max_res= 5.247E+00, max_rank=32
amen_solve: swp=4, max_dx= 4.561E-01, max_res= 2.237E+00, max_rank=42
amen_solve: swp=5, max_dx= 5.151E-01, max_res= 3.403E+00, max_rank=52
amen_solve: swp=6, max_dx= 4.186E-01, max_res= 3.818E+00, max_rank=62
amen_solve: swp=7, max_dx= 4.798E-01, max_res= 2.816E+00, max_rank=72
amen_solve: swp=8, max_dx= 4.060E-01, max_res= 1.799E+00, max_rank=82
amen_solve: swp=9, max_dx= 1.037E-02, max_res= 7.897E-02, max_rank=87
amen_solve: swp=10, max_dx= 7.175E-04, max_res= 1.934E-03, max_rank=97
amen_solve: swp=11, max_dx= 9.770E-05, max_res= 7.761E-04, max_rank=92
amen_solve: swp=12, max_dx= 8.662E-06, max_res= 2.994E-05, max_rank=96
amen_solve: swp=13, max_dx= 1.226E-05, max_res= 2.031E-05, max_rank=97
amen_solve: swp=14, max_dx= 5.315E-06, max_res= 1.115E-05, max_rank=90
amen_solve: swp=15, max_dx= 7.059E-06, max_res= 1.252E-05, max_rank=86

TT-Toolbox
=amen_solve= sweep 1, max_dx: 1.999e+00, max_res: 1.999e+00, max_rank: 12
=amen_solve= sweep 2, max_dx: 1.158e+00, max_res: 7.469e-01, max_rank: 22
=amen_solve= sweep 3, max_dx: 2.800e+00, max_res: 6.261e+01, max_rank: 32
=amen_solve= sweep 4, max_dx: 6.230e-01, max_res: 4.422e+00, max_rank: 42
=amen_solve= sweep 5, max_dx: 5.247e-01, max_res: 1.826e+00, max_rank: 52
=amen_solve= sweep 6, max_dx: 3.155e-01, max_res: 1.740e+00, max_rank: 62
=amen_solve= sweep 7, max_dx: 5.078e-01, max_res: 1.809e+00, max_rank: 72
=amen_solve= sweep 8, max_dx: 3.439e-01, max_res: 2.840e+00, max_rank: 82
=amen_solve= sweep 9, max_dx: 5.985e-01, max_res: 8.537e-01, max_rank: 92
=amen_solve= sweep 10, max_dx: 3.544e-01, max_res: 6.569e-01, max_rank: 102
=amen_solve= sweep 11, max_dx: 3.274e-01, max_res: 4.083e-01, max_rank: 110
=amen_solve= sweep 12, max_dx: 2.944e-01, max_res: 5.549e-01, max_rank: 110
=amen_solve= sweep 13, max_dx: 3.590e-01, max_res: 4.755e-01, max_rank: 110
=amen_solve= sweep 14, max_dx: 3.269e-01, max_res: 1.476e+00, max_rank: 110
=amen_solve= sweep 15, max_dx: 2.507e-01, max_res: 4.531e-01, max_rank: 110

[dolgov]

The rank selection in the residual-based truncation was different: in Matlab, it started from below and increased the rank, while in Fortran it started from the full rank and decreased it.
I made it Fortran way in both cases.
Try to pull and see if it helps.

[aiboyko]

It doesn't seem to help, although I'm pretty sure code files were updated. Unless MATLAB needs to reboot to reload files.

amen_solve2(A,b,tol,'max_full_size',15000,'kickrank',10,'x0',b,'nswp',15,...
'resid_damp',1 ,'rmax',100);

=amen_solve= sweep 1, max_dx: 1.999e+00, max_res: 1.999e+00, max_rank: 12
=amen_solve= sweep 2, max_dx: 1.442e+00, max_res: 3.240e-01, max_rank: 22
=amen_solve= sweep 3, max_dx: 2.146e+00, max_res: 5.447e+00, max_rank: 32
=amen_solve= sweep 4, max_dx: 9.203e-01, max_res: 3.143e+00, max_rank: 41
=amen_solve= sweep 5, max_dx: 3.406e-01, max_res: 2.673e+00, max_rank: 51
=amen_solve= sweep 6, max_dx: 4.106e+00, max_res: 1.783e+01, max_rank: 61
=amen_solve= sweep 7, max_dx: 4.812e-01, max_res: 2.279e+00, max_rank: 71
=amen_solve= sweep 8, max_dx: 3.319e-01, max_res: 2.111e+00, max_rank: 81
=amen_solve= sweep 9, max_dx: 4.278e-01, max_res: 1.130e+00, max_rank: 91
=amen_solve= sweep 10, max_dx: 4.316e-01, max_res: 1.017e+00, max_rank: 101
=amen_solve= sweep 11, max_dx: 2.031e-01, max_res: 2.731e-01, max_rank: 110
=amen_solve= sweep 12, max_dx: 1.977e-01, max_res: 6.424e-01, max_rank: 110
=amen_solve= sweep 13, max_dx: 1.929e-01, max_res: 3.149e-01, max_rank: 110
=amen_solve= sweep 14, max_dx: 4.607e-01, max_res: 7.568e-01, max_rank: 110
=amen_solve= sweep 15, max_dx: 4.373e-01, max_res: 1.863e+00, max_rank: 110

[dolgov]

Then send me A and b, I will try to figure it out

[aiboyko]

https://yadi.sk/d/AP9R_qV4qGt5r

[dolgov]

It may take a while...
The matrix is indefinite, so it's not a surprise that amen may not converge.
It's more intriguing why the Fortran version converges.
It's not exactly a bug, you may try any positive definite complex matrix, e.g. M+i*Laplace, it works.
But there seem to be a subtle algorithmic difference... I hope it's in our code, not in different Lapacks.

[dolgov]

OK, I'm starting to believe this is due to Matlab's Lapack...
I've updated the fort_amen_solve_mex.F90, replacing some older version of amen by dtt/ztt_amen_solve that is now used in ttpy.
The interface fort_amen_solve.m does not play a role, since we give an initial guess x0=b.
So we have the same ztt_amen_solve, but invoked from Matlab.
And...

u = fort_amen_solve(A,b,1e-5,'max_full_size',15000,'kickrank',10,'x0',b,'nswp',15);
amen_solve: swp=1, max_dx= 1.999E+00, max_res= 1.999E+00, max_rank=12
amen_solve: swp=2, max_dx= 1.330E+00, max_res= 7.468E-01, max_rank=22
amen_solve: swp=3, max_dx= 1.298E+00, max_res= 1.173E+01, max_rank=32
amen_solve: swp=4, max_dx= 3.981E-01, max_res= 2.700E+00, max_rank=42
amen_solve: swp=5, max_dx= 5.722E-01, max_res= 5.403E+00, max_rank=52
amen_solve: swp=6, max_dx= 4.127E-01, max_res= 2.105E+00, max_rank=62
amen_solve: swp=7, max_dx= 5.570E-01, max_res= 1.729E+00, max_rank=72
amen_solve: swp=8, max_dx= 2.014E+00, max_res= 2.490E+00, max_rank=82
amen_solve: swp=9, max_dx= 4.248E-01, max_res= 1.155E+00, max_rank=92
amen_solve: swp=10, max_dx= 2.744E-01, max_res= 7.353E-01, max_rank=102
amen_solve: swp=11, max_dx= 3.327E-01, max_res= 5.698E-01, max_rank=112
amen_solve: swp=12, max_dx= 2.664E-01, max_res= 8.439E-01, max_rank=122
amen_solve: swp=13, max_dx= 2.933E-01, max_res= 1.006E+00, max_rank=132
amen_solve: swp=14, max_dx= 4.066E-01, max_res= 9.980E-01, max_rank=142
amen_solve: swp=15, max_dx= 1.927E-01, max_res= 3.627E-01, max_rank=152

... it does not converge.

How robustly reproducible is the convergence in Python/Fortran?