Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

construct_partially_matrix_free for MPI #97

Open
AlexGKim opened this issue Jul 3, 2023 · 2 comments
Open

construct_partially_matrix_free for MPI #97

AlexGKim opened this issue Jul 3, 2023 · 2 comments

Comments

@AlexGKim
Copy link

AlexGKim commented Jul 3, 2023

I would like to check whether MPI StructuredMatrix construction using construct_partially_matrix_free is implemented.

I added the following test to examples/dense/testStructuredMPI.cc , which compiles but terminates with an error.

` if (world.is_root())
cout << endl << endl
<< "====================================" << endl
<< " Compression, partially matrix-free " << endl
<< "====================================" << endl;
for (auto type : types) {
options.set_type(type);
try {
auto Toeplitz_block =
[&Toeplitz](const std::vectorstd::size_t& I,
const std::vectorstd::size_t& J,
DistributedMatrix& B) {
for (std::size_t j=0; j<J.size(); j++)
for (std::size_t i=0; i<I.size(); i++)
B(i, j) = Toeplitz(I[i], J[j]);
};

        // Construct a StructuredMatrix using both a (fast)
  // matrix-vector multiplication and an element extraction
  // routine. This is mainly usefull for HSS construction, which
  // requires element extraction for the diagonal blocks and
  // random projection with the matrix-vector multiplication for
  // the off-diagonal compression.
  auto H = structured::construct_partially_matrix_free<double>
    (&grid, n, n, Tmult2d, Toeplitz_block, options);

   // print_info(world, H.get(), options);
   // check_accuracy(A2d, H.get());
   // factor_and_solve(nrhs, A2d, H.get());
   // preconditioned_solve(world, Tmult1d, H.get());
   // test_shift(nrhs, A2d, H.get());
} catch (std::exception& e) {
    if (world.is_root())
      cout << get_name(type) << " failed: " << e.what() << endl;
  }
}`
@pghysels
Copy link
Owner

pghysels commented Jul 3, 2023

Try replacing

B(i, j) = Toeplitz(I[i], J[j]);

with

B.global(i, j, Toeplitz(I[i], J[j]));

the global routine treats i and j as global indices in the distributed matrix B, and it will only assign the value if i, j is local to B.

Or equivalently, use:

if (B.is_local(i, j))
  B.global(i, j) = Toeplitz(I[i], J[j]);

then the check is done separately, so you don't need to compute the value if it's not locally required.

@AlexGKim
Copy link
Author

AlexGKim commented Jul 3, 2023

This works, thanks!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants