GlobalVector() applied to an eigenvector

[mfem]

GlobalVector() applied to an eigenvector

[mlstowell]

Hello, Mikhail,

I'm curious in what way this fails. But more to the point I do see a couple things I would change.

First, if you are using MPI_COMM_SELF then every process is performing its own separate eigenvalue solve so they each have global eigenvectors. This means that the GlobalVector() call is not actually necessary. This call allocates a new vector and populates it with data from all processors so that each processor will have a copy of the full vector. Also, since this call allocates a new Vector, technically you should free y after copying it to x (assuming to don't simply eliminate y because it doesn't seem to be needed).

Did you try:
Vector x;
modes.GetColumnReference(i, x);
x =lobpcg->GetEigenvector(i);

I hope this is useful,
Mark

[martemyev]

Hello Mark,

Thank you for your reply. I use MPI parallelism on a coarser scale, i.e. I need to solve many eigenproblems - each one within it's own subdomain. The subdomains, nevertheless, can consist of many cells, so the matrices can be quite large and sparse. Therefore, I'd like to use LOBPCG to solve the eigenproblem. The issue is that this solver is implemented for parallel setup only. So, even though I'd like to solve each eigenproblem with one thread, i.e. sequentially, I have to pretend I solve it in parallel to use LOBPCG.

You're right about the deallocation of the vector y, I had it in the code some time ago :)

And, finally, thank you for the suggestion how to fix the issue - that works great!

Cheers,
Mikhail

P.S. But it would be helpful to know why calling GlobalVector() fails. From what I understand it's not prohibited to create a full vector even though only one thread had the whole vector. I mean, it's clearly not necessary, and not efficient, but shouldn't be harmful. But what I saw was segfaults and memory corruption messages, so I'm curious why.

[martemyev]

Hello Mark,

And, finally, thank you for the suggestion how to fix the issue - that works great!

In fact, it doesn't - I checked that there were no segfaults and other nasty stuff, but in fact this approach doesn't fill up the matrix:

Vector x;
modes.GetColumnReference(i, x);
x =lobpcg->GetEigenvector(i);

The reason is because lobpcg->GetEigenvector(i) returns a vector with size = x.size()+1, i.e. one bigger than x, so when x is copy constructed from lobpcg->GetEigenvector(i), it is reallocated and no longer connected to the modes matrix.

Any suggestions?

Thanks,
Mikhail

[mlstowell]

Hi, Mikhail,

I found an error in how the eigenvectors were being allocated. In the file "linalg/hypre.cpp" the line "part[1]++;" in the HypreLOBPCG::SetOperator() method should be removed.

I've made the change in our local version but I'll need to verify a few things before changing the GitHub version.

Thanks for pointing this out. Have a nice day,
Mark

[martemyev]

Hi Mark,

Now it works fine. Thank you for fixing that!

Cheers,
Mikhail

[tzanio]

Hi Mihail,

Can we mark this issue as resolved?

Cheers,
Tzanio

[martemyev]

Yes, of course. I just thought you'd close it after updating the github repository, but it's sure resolved.

Cheers,
Mikhail