Re: [eigen] Re: Problem with conjugate gradient.

[Billy Araujo <billyaraujo@xxxxxxxxx>]

Hi Gael,

Thank you very much for your answer. I always forget the system is close to singular. Although it is symmetric.

This a matrix I obtain from a 3D poisson equation. 

Regards,

Billy.

On Sun, Mar 22, 2015 at 1:42 PM, Gael Guennebaud <gael.guennebaud@xxxxxxxxx> wrote:

Moreover your matrix is singular (rcond=1e-18) so there is no unique solution. Actually, enforcing Eigen::ConjugateGradient to work in non-symmetric mode with ILUT give me another but as accurate solution:

x^T = [0.199853   0.14473  0.103595 0.0624614 0.0073383  0.153749  0.130412  0.103595 0.0767788 0.0534419  0.131143  0.119665  0.103595 0.0875256 0.0760479  0.120043  0.113539  0.103595 0.0936523 0.0871481  0.115452   0..11086  0.103595 0.0963306 0.0917393]

relative error: (A*x-b).norm()/b.norm() == 1.32655e-15

Respective code:

  ConjugateGradient<SpMat,Lower|Upper,IncompleteLUT<double> > cg_ilut(A);
  x = cg_ilut.solve(b);


cheers,

gael

On Sun, Mar 22, 2015 at 2:27 PM, Gael Guennebaud <gael.guennebaud@xxxxxxxxx> wrote:
>
> ILU is not compatible with Eigen::ConjugateGradient because the later is for symmetric problems while ILU does not preserve symmetry.
>
> cheers,
> gael
>
> On Sat, Mar 21, 2015 at 6:33 PM, Billy Araujo <billyaraujo@xxxxxxxxx> wrote:
>>
>> Just a minor correction:
>>
>> Laspack gives 25 values:
>>
>> [1] 0.096709627074711618 double
>> [2] 0.041738859730487146 double
>> [3] 0.00067511123237074346 double
>> [4] -0.040388637958504281 double
>> [5] -0.095359407662660600 double
>> [6] 0.050941442361153874 double
>> [7] 0.027508948358719330 double
>> [8] 0.00067511042278063161 double
>> [9] -0.026158728561469002 double
>> [10] -0.049591224744352222 double
>> [11] 0.028279326736092902 double
>> [12] 0.016767643972775401 double
>> [13] 0.00067510968097885340 double
>> [14] -0.015417425476270046 double
>> [15] -0.026929109622622786 double
>> [16] 0.017154602993941025 double
>> [17] 0.010635301124244559 double
>> [18] 0.00067510921541819892 double
>> [19] -0.0092850833115109520 double
>> [20] -0.015804386045630902 double
>> [21] 0..012553916220840610 double
>> [22] 0.0079532300763286365 double
>> [23] 0.00067510900276445361 double
>> [24] -0.0066030125504719286 double
>> [25] -0.011203699327577501 double
>>
>> Steps to reproduce in Matlab (quick reference):
>>
>> A = mmread('A.mat')
>> b = mmread('b.vec')
>> x = A\b
>>
>>
>>
>>
>> On Sat, Mar 21, 2015 at 5:26 PM, Billy Araujo <billyaraujo@xxxxxxxxx> wrote:
>>>
>>>
>>> I have solved a system using Laspack using conjugate gradient and ILU preconditioner (tol. 1E-12) and it gives the following solution:
>>>
>>> [1] 0.096709627074711618 double
>>> [2] 0.041738859730487146 double
>>> [3] 0.00067511123237074346 double
>>> [4] -0.040388637958504281 double
>>> [5] -0.095359407662660600 double
>>> [6] 0.050941442361153874 double
>>> [7] 0.027508948358719330 double
>>> [8] 0.00067511042278063161 double
>>> [9] -0.026158728561469002 double
>>> [10] -0.049591224744352222 double
>>> [11] 0.028279326736092902 double
>>> [12] 0.016767643972775401 double
>>> [13] 0.00067510968097885340 double
>>> [14] -0.015417425476270046 double
>>> [15] -0.026929109622622786 double
>>> [16] 0.017154602993941025 double
>>> [17] 0.010635301124244559 double
>>> [18] 0.00067510921541819892 double
>>> [19] -0.0092850833115109520 double
>>> [20] -0.015804386045630902 double
>>> [21] 0.012553916220840610 double
>>> [22] 0.0079532300763286365 double
>>> [23] 0.00067510900276445361 double
>>> [24] -0.0066030125504719286 double
>>>
>>> I have attached A and b files and in Matlab the solution of x = A\b gives:
>>>
>>>    0.0963
>>>     0.0411
>>>    -0.0000
>>>    -0.0411
>>>    -0.0963
>>>     0.0502
>>>     0.0268
>>>    -0.0000
>>>    -0.0268
>>>    -0.0502
>>>     0.0275
>>>     0.0161
>>>    -0.0000
>>>    -0.0161
>>>    -0.0275
>>>     0.0164
>>>     0.0099
>>>    -0.0000
>>>    -0.0099
>>>    -0.0164
>>>     0.0119
>>>     0.0073
>>>    -0.0000
>>>    -0.0073
>>>    -0.0119
>>>
>>> In Eigen however the solution of the same system using ConjugateGradient and ILU preconditioner is:
>>>
>>>             ConjugateGradient<SparseMatrix<Real>, Lower, IncompleteLUT<Real> > solver(A);
>>>
>>>             solver.setTolerance(1E-12);
>>>             x = solver.solve(b);
>>>
>>> x:
>>>
>>> 0.552865
>>> 0.497742
>>> 0.456608
>>> 0.415474
>>> 0.360351
>>> 0.506762
>>> 0.483425
>>> 0.456608
>>> 0.429791
>>> 0.406454
>>> 0.484156
>>> 0.472678
>>> 0.456608
>>> 0.440538
>>> 0.429060
>>> 0.473055
>>> 0.466551
>>> 0.456608
>>> 0.446665
>>> 0.440161
>>> 0.468464
>>> 0.463873
>>> 0.456608
>>> 0.449343
>>> 0.444752
>>>
>>> This is very different from Matlab and Laspack. What am I doing wrong? Am I calling Eigen correctly?
>>>
>>> Please find in attachment the files to reproduce the solution. Steps:
>>>  
>>>
>>
>