Discrete matrices are not strictly symmetric

General Discussion
1

Dear all,
Hi! Recently, I found that discrete matrices are not strictly symmetric, but in theory it should be symmetric. How can I get the strictly symmetric one?

varf a([uu,vv],[w,s])=
int2d(Th)(
2mu(dx(uu)dx(w)+dy(vv)dy(s)+ ((dx(vv)+dy(uu))(dx(s)+dy(w)))/2 )
+ lambda(dx(uu)+dy(vv))*(dx(w)+dy(s))
)

  • on(1,uu=0,vv=0)
    ;

varf b([uu,vv],[w,s])=
int2d(Th)(uuw + vvs) ;

matrix A= a(Vh,Vh,solver=“SPARSESOLVER”);
matrix B= b(Vh,Vh,solver=CG,eps=1e-20);

1、I tried to At=(A+A’)/2, but It doesn’t work.
2、I tried to use MatScale, but It is useless.
3、I tried this,
matrix A= a(Vh,Vh,solver=“SPARSESOLVER”);
matrix B= b(Vh,Vh,solver=CG,eps=1e-20);
A = (A+(A’));
Mat At=A;
At *= 0.5;
B = (B+B’);
Mat Bt=B;
Bt *= 0.5;

But, Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range

Are there other methods???

Best regards,
Baishi.

2

What do you mean by “not strictly symmetric”?

3

Hi, prj, thanks for your reply!

I export the discrete matrices to MATLAB, and Not Strictly Symmetric means that

          >> max(max(A-A.'))
          ans =
               1.136868377216160e-13
          
          or norm(A-A')!=0
4

It’s nearly symmetric up to machine epsilon. You are likely loosing precision while exporting the matrix to MATLAB.

5

Sorry, I am not sure that.
I export them to SLEPc as binary file, however, can’t use many eps, such as lobpcg.

   [0]PETSC ERROR: The solver 'lobpcg' cannot be used for non-symmetric problems

So I want to know how to process manually makes it symmetric in FreeFem?

6

You should not use LOBPCG, much less robust than Krylov–Schur, but you need to add the option -eps_gen_hermitian.

7

Thanks for your reminding, and I still want to know how to process manually makes it symmetric in FreeFem?

8

It is symmetric, up to machine epsilon.

10 
mesh Th = square(10, 10);
varf vPb(u, v) = int2d(Th)(dx(u)*dx(v) + dy(u)*dy(v));
fespace Vh(Th, P1);
matrix A = vPb(Vh, Vh);
matrix B = A+A';
B = 0.5 * B;
11

That doesn’t work…
2d-beam-eigen.edp (4.2 KB)
Segmentation fault (core dumped)

    matrix AA= a(Vh,Vh,solver="SPARSESOLVER"); 
    matrix C = AA+AA';
    matrix A = 0.5 * C;
    matrix BB= b(Vh,Vh,solver=CG,eps=1e-20); 
    C = BB+BB';
    matrix B = 0.5 * C;
12

That does work… 2d-beam-eigen.edp (2.7 KB)

13

Thanks, but another problem occurs.
Segmentation fault (core dumped) when call the subroutine Eigenvalue

I don’t know what happen :sweat:

I want to compare the results of Eigvalue(ARPACK) in FreeFem with ARPACK in SLEPc

14

I don’t recommend the use of ARPACK (here is a prime example of why), and I don’t know how to fix this.

15

Well, thank you very much for your help today

16

Hi, prj, there still exist the problem that the binary file without processing was considered non-symmetric by SLEPc. Because of this, I can’t use many eps :sob:

17

Why are you not calling SLEPc through FreeFEM? In which case it won’t complain if you give the proper option I mentioned earlier.