I try to adapt the example “mf-2d-SLEPc.edp” to compute singular values for a multi-dimensional problem. Namely, I want to compute the singular values of the operator MF
where Qu is a matrix created on PuPuPp and Qf is a matrix created on PuPu where:
func PuPuPp=[P2, P2, P1]; func PuPu=[P2, P2];
The associated spaces are:
fespace Uvvp(Th, PuPuPp); fespace Uvv(Th, PuPu);
I wrote the following lines to compute the singular values and extract the left and right vectors (the example helmholtz-2d-SLEPc-complex.edp helped me here):
SLEPc computes singular values that I can access in the array values. However the left and right vectors are all zero. I can observe that for instance by adding the following lines:
Is there something wrong in the syntax to call SVDSolve to define lvectors and rvectors? Or is there something wrong about how I access them one the computation is done?
Actually by I modified the example “mf-2d-SLEPc.edp” (see attachment) mf-2d-SLEPc.edp (2.3 KB) and I observe the same issue. So it is not due to the fact that my example is multi-dimensional. I guess I am doing something wrong when accessing the vectors once SVDSolve was called.
I see. There is no domain decomposition attached to MF, so you cannot use output vectors stored in FreeFEM numbering (lvectors and rvectors). Instead, you need to use vectors stored in PETSc numbering (larray and rarray), and if you need to go back to FreeFEM numbering, use ChangeNumbering(..., inverse = true). Does that make sense to you?
MF is not indirectly aware of the decomposition made by buildDmesh(Th); ? Is there a way to attach the domain decomposition to MF ?
I roughly understood from “mf-2d-SLEPc.edp” how to use ChangeNumbering. I understand that I should replace lvectors = lvec by larray = lvec where lvec is now real[int,int] of size (Vh.ndof, 10) (as in Composite eigenvectors with SLEPc)
I tried real[int][int] but it doest not seem to be correct.
What you need is real[int, int] (bidimensional array). real[int][int] is for arrays of arrays, but here, all columns (or rows for that matter) are of the same size.
But then using lvectors on C should work? It does not to be the case as you can have observed from the above attached file.
I would like to extend the example “mf-2d-SLEPc.edp” to show how to use larray/rarray and then push it on the FreeFem Git as I think it adds some valuable explanation on how to use SVDSolve. I would like to add three asserts:
same left vectors with and without matrix-free resolution
same right vectors with and without matrix-free resolution
matrix x (left vector) = (singular value) x (right vector)
Please find enclosed what I did so far: mf-2d-SLEPc.edp (2.9 KB) . It is not functional yet: the comparison of left vectors fails. Moreover I struggle to change the numbering of the right vectors. If I understand correctly the matrix B should be used to change the numbering of the right vectors, but the re-numbering only works if I use the matrix A. Could you, please, help me to correct this?
I compiled FreeFem with the proposed correction in SLEPc-code.hpp. The assert on the left vectors now successes. However, it still fail to handle the right vectors, mf-2d-SLEPc.edp (3.1 KB)
I do not understand what I do wrong but I can not renumber the right vectors using the matrix B.
Another good catch! Sorry, I’m not using SVDSolve() much, so thank you for your careful investigation. This part of SLEPc-code.hpp is wrong, rarray can’t be of the same size as array in the rectangular case. That would explain why the ChangeNumbering() with A works (correct dimension) while it doesn’t with B (incorrect dimension). But you are correct, it should be the other way around.
I’m AFK until the early afternoon, but I’ll try to fix this later today. Thank you for bearing with me .
It should be good now. I’ve fixed mf-2d-SLEPc.edp (3.1 KB) accordingly, is it alright on your side if I commit this modified file on the FreeFEM official repository?
Actually there is one more question from my side. The last assert corresponds to check that C v = sigma u in the PETSc world. If I understand correctly, this should be true as well in the FreeFem world, isn’t it ?
Such test in FreeFem world would read: temp1 = MF * vMF[]; temp2 = valuesMF[idSV] * uMF[]; temp2 -= temp1; assert(temp2.linfty < 1.0e-6);
I tried that but that fails, did I do something wrong?
Best,
Lucas
EDIT:
The following assert (that fails) convinced me that something is wrong when performing the matrix vector multiplication in the FreeFem world:
Can’t do MF * because MF is rectangular and I haven’t done the necessary plumbing to handle this case. Must use MatMult instead. If you had a FreeFEM compiled with debugging turned on, you would see an assert error when doing either temp1 = MF * vMF[]; or temp2 = MF * vMF[];, whatever comes next is garbage. (maybe not that simple, let me check and I’ll get back at you)
Edit: you forgot to resize temp2 before MF *, the following seems to work.
MatMult(MF, rvecMF(:,idSV), temp2);
ChangeNumbering(A, temp1, temp2, inverse = true, exchange = true);
temp2.resize(uMF[].n); // should be done automatically, but necessary plumbing not done yet...
temp2 = MF * vMF[];
temp2 -= temp1;
assert(temp2.linfty < 1.e-6);
And just a remark for those like that would be interested (like me ^^), it is possible to stay in the “PETSc world” in prodFunc. In such case, prodFunc would read:
func real[int] prodFunc(real[int]& up) {
real[int] u(A.n);
MatMultTranspose(RPETSc, up, u); // where RPETSc is defined by Mat RPETSc(B, A, R);
MatMult(A, u, up);
return up;
}
I guess this should improve the performance as we avoid two re-numbering.
It seems that the correction you proposed broke something when using EPSolve( ... , array = vec);. I came to this conclusion by comparing the result of the same script executed with the v4.8 and the develop branch. In the last case, vec has wrong entries (either really small or really large) while with the v4.8 I obtain the correct eigenvector shape.
I’m actually not quite sure how that code ever worked… I’ll fix a push after diner, I’ve got it locally but I need to figure out what was previously going on.
.
