Inverting a matrix object?

General Discussion
1

Hi there,

I’m not sure if this has been answered before, and probably seems like a trivial question, but suppose that I have the following code

varf a(u,v) = int2d(Th)(dx(u)*dx(v) + dy(u)*dy(v));
matrix K = a(Vh,Vh);

and that I would like to invert the matrix K. I’ve tried inv(K) from fflapack and K^-1 and both did not work. I also skimmed through the documentation but with no luck. How do I invert a matrix object in general?

Thanks in advance!

2

Why do you need to invert the matrix? If you are using the matrix to solve a linear system you should definitely not invert it.

You cannot invert a matrix with fflapack – you will need to first convert it to a dense matrix, i.e. real[int,int].
https://doc.freefem.org/examples/developers.html#matrix-inversion

3

Thank you for the message and sorry for the late reply. I’m doing a time dependent simulation where I need the inverse of K at each iteration. So instead of solving K^-1*u at each step, I thought I’d just compute Kinv=K^-1 and then do Kinv*u directly to speed things up.

Also, how do you convert to a dense matrix? I tried the following incredibly slow solution

int col = K.n;
int row = K.m;
real[int,int] Kin(row,col); 
for(int i=0; i<col; i++){	
	for(int j=0; j<row; j++){	
		Kin(j,i)=K(j,i);	
	}				
}					
inv(Kin);

Not sure why copying takes thing long for a 4096\times 4096 matrix!

4

Computing Kinv and then Kinv*u is in no way faster than computing K^-1*u. The L and U factors of the decomposition of K are computed only once, and then it’s just successive forward and backward solves, which are basically as fast as matrix-vector products.
Doing for loop in FreeFEM is slow, you should replace your loop by simply doing Kin = K;.

5

So if I understand correctly, FreeFEM stores the LU-factorization of K after the first use of K^-1*u?

Also, Kin = K generates an error as K is of matrix type and Kin is of type real[int][int].

6

You understand correctly, but it’s not specific to FreeFEM, it’s almost all sparse exact factorization solvers that behave like this (SuiteSparse, MUMPS, PARDISO…).

It does not fail with a real[int, int], which is the type you used in your original snippet.

7

Sorry I meant real[int,int] in my previous post. Also I confirm that

varf a(u,v) = int2d(Th)(dx(u)*dx(v) + dy(u)*dy(v));
matrix K = a(Vh,Vh);
real[int,int] Kin(K.m,K.n);
Kin = K;

creates an error.

8

Right, this does not work, you can simply do for[I,J,C : K] Kin(I,J) = C;.

1 Like