[boost] [uBlas] [Tensor] multiprecision and FFT?

[Janek Kozicki via Boost]

Cem Bassoy via Boost said: (by the date of Wed, 6 Jan 2021 10:32:15 +0100)

> Please consider to use and contribute to *Boost.uBlas*


Hi,

Does it work with boost::multiprecision types? In particular I am
interested in types listed here:

http://yade-dem.org/doc/HighPrecisionReal.html


Does it support FFT ? In particular I am interested in replacing my
very crude NDimTable class with your code:

https://gitlab.com/cosurgi/trunk/-/blob/addQuantumMechanics_FixEnum_FixRebase1/lib/base/NDimTable.hpp

I am about to start refactoring this part of the code, when I noticed your post.

best regards
Janek


--
Janek Kozicki, PhD. DSc. Arch. Assoc. Prof.
Gdańsk University of Technology
Faculty of Applied Physics and Mathematics
Department of Theoretical Physics and Quantum Information
--
http://yade-dem.org/
http://pg.edu.pl/jkozicki (click English flag on top right)

_______________________________________________
Unsubscribe & other changes: http://lists.boost.org/mailman/listinfo.cgi/boost

[Cem Bassoy via Boost]

Am Mi., 6. Jan. 2021 um 23:19 Uhr schrieb Janek Kozicki via Boost <
bo...@lists.boost.org>:

> Cem Bassoy via Boost said: (by the date of Wed, 6 Jan 2021 10:32:15
> +0100)
>
> > Please consider to use and contribute to *Boost.uBlas*
>
>
> Hi,
>
> Does it work with boost::multiprecision types? In particular I am
> interested in types listed here:
>
> http://yade-dem.org/doc/HighPrecisionReal.html
>
>

It does work with *boost::multiprecision*:
Have a look @ https://godbolt.org/z/cMKc9T

>
> Does it support FFT ?

You can get the pointer (like for vector) to the underlying contiguous
memory (you do not own the memory):

using format_t = boost::numeric::ublas::column_major;

using tensor_t = boost::numeric::ublas::tensor<float,format_t>;
auto A = tensor_t(shape{3,4,2},2);auto ap = A.data();

You can also use the standard c++ library for convenience:

std::for_each(A.begin(), A.end(), [](auto& a){ ++a; });

If you do not want to use the Einstein-notation, you can as well use either
the prod function:

// C3(i,l1,l2) = A(i,j,k)*T1(l1,j)*T2(l2,k);

q = 3u;

tensor_t C3 = prod(prod(A,matrix_t(m+1,n[q-2],1),q-1),matrix_t(m+2,n[q-1],1),q);

"prod" uses internally a C-like interface which will not allocate memory
at all.

ttm(m, p,

c.data(), c.extents().data(), c.strides().data(),

a.data(), a.extents().data(), a.strides().data(),

bb, nb.data(), wb.data());

Cheers,
Cem

[Janek Kozicki via Boost]

Cem Bassoy via Boost said: (by the date of Thu, 7 Jan 2021 12:43:47 +0100)

> Have a look @ https://godbolt.org/z/cMKc9T

wow, This is very good. Thank you.

> You can get the pointer (like for vector) to the underlying contiguous
> memory (you do not own the memory):

This is exactly what I need. FFT does in-place transform by modifying the memory.

> using format_t = boost::numeric::ublas::column_major;

FFT does its work on row_major format. You have this one so
I assume that you also have the other one ;)

However I needed to write a custom memory allocator (also very crude):

https://gitlab.com/cosurgi/trunk/-/blob/addQuantumMechanics_FixEnum_FixRebase1/lib/base/NDimTable.hpp#L96
https://gitlab.com/cosurgi/trunk/-/blob/addQuantumMechanics_FixEnum_FixRebase1/lib/base/FFTW3_Allocator.hpp

Can I provide a custom allocator to tensor_t ?

In fact maybe during past 5 years some better C++ interface for
libfftw appeared. I will have to check :)
maybe even boost library has an FFT interface? :)

If not then I will need to refactor this crude FFTW3_Allocator memory allocator.

> using tensor_t = boost::numeric::ublas::tensor<float,format_t>;
> auto A = tensor_t(shape{3,4,2},2);auto ap = A.data();
> You can also use the standard c++ library for convenience:
> std::for_each(A.begin(), A.end(), [](auto& a){ ++a; });
> If you do not want to use the Einstein-notation, you can as well use either

I'm definitely going to use the Einstein notation :)

> the prod function:
> // C3(i,l1,l2) = A(i,j,k)*T1(l1,j)*T2(l2,k);
> q = 3u;
> tensor_t C3 = prod(prod(A,matrix_t(m+1,n[q-2],1),q-1),matrix_t(m+2,n[q-1],1),q);
>
> "prod" uses internally a C-like interface which will not allocate memory
> at all.
>
> ttm(m, p,
> c.data(), c.extents().data(), c.strides().data(),
> a.data(), a.extents().data(), a.strides().data(),
> bb, nb.data(), wb.data());

I have another question about indexing the tensor.

Can I have a custom index redirection? I mean, when I type code to
access zero-th element of 1D array:

tensor[0] = ... ;

I mean to access the element which is located exactly in the middle
of the reserved memory region. This is because libfftw implementation
is not compatible with other parts of the quantum dynamic time
propagation algorithm. This incompatibility is possible because FT is
intrinsically periodic. Currently I have to use this function
shiftByHalf():

https://gitlab.com/cosurgi/trunk/-/blob/addQuantumMechanics_FixEnum_FixRebase1/lib/base/NDimTable.hpp#L1107

Which shifts by half an N-dimensional table. So that libfftw sees the
data which it will process correctly. And doing a single calculation
step is rather slow because after doing FFT it has to be shifted back:

https://gitlab.com/cosurgi/trunk/-/blob/addQuantumMechanics_FixEnum_FixRebase1/lib/base/NDimTable.hpp#L1436

this makes calcuations very slow. If I could use some kind of index
redirection, that will make all iterators treat the data as being
shifted by half:

* zero-th element: in the middle of reserved memory.
* half: at the end of reserved memory
* half+1: at the beginning of reserved memory

then all other parts of my algorithm will not need to be modified,
because they solely work with iterators.

(Also I need to check if it's now possible to tell FFT to treat data
differently, it wasn't possible 5 years ago, then I wouldn't need all
these complications with indexing)

Thank you, your code is very promising.
Janek


--
# Janek Kozicki http://janek.kozicki.pl/

[Cem Bassoy via Boost]

Am Do., 7. Jan. 2021 um 14:24 Uhr schrieb Janek Kozicki via Boost <
bo...@lists.boost.org>:

> Cem Bassoy via Boost said: (by the date of Thu, 7 Jan 2021 12:43:47
> +0100)
>
> > Have a look @ https://godbolt.org/z/cMKc9T
>
> wow, This is very good. Thank you.
>
> > You can get the pointer (like for vector) to the underlying contiguous
> > memory (you do not own the memory):
>
> This is exactly what I need. FFT does in-place transform by modifying the
> memory.
>
> > using format_t = boost::numeric::ublas::column_major;
>
> FFT does its work on row_major format. You have this one so
> I assume that you also have the other one ;)
>
> However I needed to write a custom memory allocator (also very crude):
>
>
> https://gitlab.com/cosurgi/trunk/-/blob/addQuantumMechanics_FixEnum_FixRebase1/lib/base/NDimTable.hpp#L96
>
> https://gitlab.com/cosurgi/trunk/-/blob/addQuantumMechanics_FixEnum_FixRebase1/lib/base/FFTW3_Allocator.hpp
>
> Can I provide a custom allocator to tensor_t ?
>

Yes, you can. We work with *std::vector* as the underlying base type.

template<class T, class F = first_order, class A =
std::vector<T,std::allocator<T>> >

class tensor

Just specialize your tensor type with your own allocator.
You can also provide your own base type instead of *std::vector* but
it must model the sequential container concept.
using other_tensor_t =
tensor<other_numeric_type,first_order,std::vector<T,other_allocator>>;

In fact maybe during past 5 years some better C++ interface for
> libfftw appeared. I will have to check :)
> maybe even boost library has an FFT interface? :)
>
> If not then I will need to refactor this crude FFTW3_Allocator memory
> allocator.
>
> > using tensor_t = boost::numeric::ublas::tensor<float,format_t>;
> > auto A = tensor_t(shape{3,4,2},2);auto ap = A.data();
> > You can also use the standard c++ library for convenience:
> > std::for_each(A.begin(), A.end(), [](auto& a){ ++a; });
> > If you do not want to use the Einstein-notation, you can as well use
> either
>
> I'm definitely going to use the Einstein notation :)
>

Glad it helps.
Please note that contractions are not optimized. For small dimensions, the
runtime should be tolerable.
For larger dimensions and rank, please let me know - we are working on fast
tensor contractions (e.g. https://github.com/bassoy/ttv)

>
> > the prod function:
> > // C3(i,l1,l2) = A(i,j,k)*T1(l1,j)*T2(l2,k);
> > q = 3u;
> > tensor_t C3 =
> prod(prod(A,matrix_t(m+1,n[q-2],1),q-1),matrix_t(m+2,n[q-1],1),q);
> >
> > "prod" uses internally a C-like interface which will not allocate memory
> > at all.
> >
> > ttm(m, p,
> > c.data(), c.extents().data(), c.strides().data(),
> > a.data(), a.extents().data(), a.strides().data(),
> > bb, nb.data(), wb.data());
>
> I have another question about indexing the tensor.
>
> Can I have a custom index redirection? I mean, when I type code to
> access zero-th element of 1D array:
>
> tensor[0] = ... ;
>

Yes, random access with read and write.
You can access tensor data with a single index like *t[0] *or*
t[(n*m*...*k)/2]*.
You can also use multi-indices for convenience, e.g. *t.at
<http://t.at>(5,2,4,6)*

> I mean to access the element which is located exactly in the middle
> of the reserved memory region. This is because libfftw implementation
> is not compatible with other parts of the quantum dynamic time
> propagation algorithm. This incompatibility is possible because FT is
> intrinsically periodic. Currently I have to use this function
> shiftByHalf():
>
>
> https://gitlab.com/cosurgi/trunk/-/blob/addQuantumMechanics_FixEnum_FixRebase1/lib/base/NDimTable.hpp#L1107
>
> Which shifts by half an N-dimensional table. So that libfftw sees the
> data which it will process correctly. And doing a single calculation
> step is rather slow because after doing FFT it has to be shifted back:
>
>
> https://gitlab.com/cosurgi/trunk/-/blob/addQuantumMechanics_FixEnum_FixRebase1/lib/base/NDimTable.hpp#L1436
>
> this makes calcuations very slow. If I could use some kind of index
> redirection, that will make all iterators treat the data as being
> shifted by half:
>

I am not really sure if I understand what you mean by *index redirection*.
You can always advance pointers and iterators of tensor.

> * zero-th element: in the middle of reserved memory.
>
* half: at the end of reserved memory
> * half+1: at the beginning of reserved memory
>
>

If I understand you correctly, you would like to create two valid ranges.
Let t be a p-order tensor *t* and n, m, k be its dimensions, where p =
n*m*...*k with n,m,...,k > 1

Then *[first0,last0)* with
*auto first0 = std::advance(std::begin(t), p/2); *// auto first0 =
t.begin()+p/2;

*auto last0 = std::end(t);*
should be your first valid range.

Then *[first1,last1)*

*auto first1 = std::begin(t);*
*auto last1 = std::advance(std::begin(t), p/2-1);*
should be your second valid rang.

>
> then all other parts of my algorithm will not need to be modified,
> because they solely work with iterators.
>

If you need you can also work with iterators:
*std::advance(std::begin(t),5)* as you are used with *std::vector* .
Note (for simplified interface with the fft library you can the shifting
and advancing using raw pointers of the tensor *std::advance(t.data(),5) *or
simply *t.data()+5.*

> (Also I need to check if it's now possible to tell FFT to treat data
> differently, it wasn't possible 5 years ago, then I wouldn't need all
> these complications with indexing)
>
> Thank you, your code is very promising.
> Janek
>

Thanks.
Hope I could help you.