# Extract of left and right vectors of multidimensional SVD

Hi,

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`

`Mat<complex> MF(Qu, Qf, prodFunc, transpose = prodFuncTranspose);`

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):

`int nsv = 4;`
`real[int] values(nsv);`
`Uvv<complex>[int] [rsolx,rsoly](nsv);`
`Uvvp<complex>[int] [lsolx, lsoly, lsolp](nsv);`
` int nconv = SVDSolve(MF, sparams = "-svd_largest -svd_view_values -svd_type cyclic -svd_nsv "+nsv, values = values, lvectors = lsolx, rvectors = rsolx);`

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:

for (int idSv = 0; idSv < nconv; ++idSv)
{
cout << "MAX " << (rsolx[idSv][]).linfty << " / " << (lsolx[idSv][]).linfty << endl;
}

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?

Best,

Lucas

Could you please provide a full MWE? You can send it in private if it’s closed-source code. Thanks.

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?

Not completely!

There is no domain decomposition attached to `MF`

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)

Best,

Lucas

MF is not indirectly aware of the decomposition

Correct. `C` is, not `MF`.

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.

Correct. `C` is, not `MF` .

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?

Good catch! Should read `(*rarray)(':', i) = cpy;`, not `(*array)(':', i) = cpy;` in SLEPc-code.hpp.

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?

1 Like

Thanks for the quick code correction!

Yes please commit the modified example

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:

`MatMult(MF, rvecMF(:,idSV), temp2);`
`ChangeNumbering(A, temp1, temp2, inverse = true, exchange = true);`
`temp2 = MF * vMF[];`
`temp2 -= temp1;`
`assert(temp2.linfty < 1.e-6);`

It is not allowed to do `MF * vMF[]` ?

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);
``````

I see, thank you for the explanation

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.

Indeed, it’s always best to stick as much as possible to PETSc numbering. This was merely to show the intent of `ChangeNumbering()` & friends

1 Like

Hi!

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.

Best,

Lucas

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.

This should be fixed in this commit, don’t really know what I was thinking of when writing that in the first place…

That works all fine, thanks!