amgcl icon indicating copy to clipboard operation
amgcl copied to clipboard

How to deal with multiple right hand sides.

Open g873249824 opened this issue 3 years ago • 8 comments

Hi Professor, I have a basic question. If I have vectors rhss and xs as arguments to solver operator () as (rhss,xs), which are multiple right hand sides (continusly stored) and solutions, I want to handle rhss and xs inside the header file by changing the operator (). Can I declare a vector reference to part of rhss or xs without additional memory storage? or something alternative?

{
vector &rhs=rhss(0,n-1);
vector &x=xs(0,n-1);
//do some iteration cycles for the first rhs.
}
{
vector &rhs=rhss(n,2n-1);
vector &x=xs(n,2n-1);
//do some iteration cycles for the second rhs.
}// maybe i will continue  the iteration cycle for first rhs

How can i implement this in a elegant way for several backends? With best regards.

g873249824 avatar Oct 31 '22 12:10 g873249824

Hey,

You should be able to use amgcl::make_iterator_range as in

for(int i = 0; i < nv; ++i) {
    auto b = make_iterator_range(&rhss[i*n], &rhss[i*(n+1)]);
    auto x = make_iterator_range(&xs[i*n], &xs[i*(n+1)]);
    solve(b, x);
}

ddemidov avatar Oct 31 '22 13:10 ddemidov

Thanks for your help. I have just tried that, i think it works. Sincerely!

g873249824 avatar Oct 31 '22 13:10 g873249824

Sorry to bother you again, I was trying to do something like operator = or - like

vector t = b;
// or vector t= t-b;

I have tried converting it to numa_vector<rhs_type> and functions like backend::copy or backend::axby and I think I tried the wrong usage (because it is an iterator type). What is the right way to do this if I use 'amgcl::make_iterator_range'?

g873249824 avatar Nov 01 '22 02:11 g873249824

The following works for me:

#include <amgcl/backend/builtin.hpp>

int main() {
    std::vector<double> x(10, 1), y(10, 2);

    auto X = amgcl::make_iterator_range(x.begin(), x.end());
    auto Y = amgcl::make_iterator_range(y.begin(), y.end());

    // y = x
    amgcl::backend::copy(X, Y);

    std::cout << "1. y = " << y[0] << std::endl;

    // y = y - x
    amgcl::backend::axpby(-1, X, 1, Y);

    std::cout << "2. y = " << y[0] << std::endl;
}
1. y = 1
2. y = 0

ddemidov avatar Nov 01 '22 04:11 ddemidov

I will try that. Thanks for your time and patience.

g873249824 avatar Nov 01 '22 04:11 g873249824

I don't know if I can implement this for several backends, such as Eigen backend or cuda.

g873249824 avatar Nov 01 '22 12:11 g873249824

  • It should work for cpu-based backends (you can get pointers to the raw data from Eigen vector)
  • I don't think cuda thrust library has range concept, so its not possible with the cuda backend (you would have to copy the current vector chunk into a smaller vector)
  • It could probably work with vexcl backend using slices (https://vexcl.readthedocs.io/en/latest/expressions.html#slicing) but I have not tried this.

ddemidov avatar Nov 01 '22 13:11 ddemidov

Really thanks for your comment. Maybe i will try to copy the data to a smaller vector so I can try to compile for different backends.

g873249824 avatar Nov 01 '22 13:11 g873249824