convolution_timestep_omp() [1/4]

template<typename T >

void cntr::convolution_timestep_omp	(	int	omp_num_threads,
		int	tstp,
		herm_matrix< T > &	C,
		herm_matrix< T > &	A,
		herm_matrix< T > &	Acc,
		function< T > &	ft,
		herm_matrix< T > &	B,
		herm_matrix< T > &	Bcc,
		T	beta,
		T	h,
		int	SolveOrder
	)

Returns convolution \(C = A\ast f\ast B\) at a given time step

Purpose

Computes contour convolution C=A*f*B of the objects with class 'herm_matrix' at a given time step 't=nh'. Works for a scalar and square matrices. openMP parallelized version.

Parameters

omp_num_threads	[int] The number of openMP threads to be used. Set to the number of threads in the current team.
tstp	[int] time step
C	[herm_matrix] Matrix to which the result of the convolution on Matsubara axis is given
A	[herm_matrix] contour Green's function
Acc	[herm_matrix] hermitian conjugate of A
ft	[herm_matrix] function \(f\)
B	[herm_matrix] contour Green's function
Bcc	[herm_matrix] hermitian conjugate of B
beta	inversed temperature
h	time step interval
SolveOrder	[int] integrator order

Definition at line 5231 of file cntr_convolution_impl.hpp.

                                                           {
    int ntmin = (tstp == -1 || tstp > SolveOrder ? tstp : SolveOrder);
    if (tstp < -1)
        return;
    int size1 = A.size1();
    std::complex<T> *fttemp;
    assert(C.size1()==size1);
    assert(ft.size1()==size1);
    assert(B.size1()==size1);
    assert(Bcc.size1()==size1);
    assert(Acc.size1()==size1);
    assert(SolveOrder>=0 && SolveOrder <=5);
    assert(C.ntau()>=SolveOrder);
    assert(C.nt()>=ntmin);
    assert(A.nt()>=ntmin);
    assert(Acc.nt()>=ntmin);
    assert(B.nt()>=ntmin);
    assert(Bcc.nt()>=ntmin);
    assert(ft.nt()>=ntmin);
    assert(C.ntau()==A.ntau());
    assert(C.ntau()==Acc.ntau());
    assert(C.ntau()==B.ntau());
    assert(C.ntau()==Bcc.ntau());
    C.set_timestep_zero(tstp);
    fttemp = (tstp == -1 ? ft.ptr(-1) : ft.ptr(0));

    switch (size1) {
    case 1:
        incr_convolution_omp<T, herm_matrix<T>, 1>(omp_num_threads, tstp, CPLX(1, 0), C, A,
                                                   Acc, ft.ptr(-1), fttemp, B, Bcc,
                                                   integration::I<T>(SolveOrder), beta, h);
        break;
    case 2:
        incr_convolution_omp<T, herm_matrix<T>, 2>(omp_num_threads, tstp, CPLX(1, 0), C, A,
                                                   Acc, ft.ptr(-1), fttemp, B, Bcc,
                                                   integration::I<T>(SolveOrder), beta, h);
        break;
    case 3:
        incr_convolution_omp<T, herm_matrix<T>, 3>(omp_num_threads, tstp, CPLX(1, 0), C, A,
                                                   Acc, ft.ptr(-1), fttemp, B, Bcc,
                                                   integration::I<T>(SolveOrder), beta, h);
        break;
    case 4:
        incr_convolution_omp<T, herm_matrix<T>, 4>(omp_num_threads, tstp, CPLX(1, 0), C, A,
                                                   Acc, ft.ptr(-1), fttemp, B, Bcc,
                                                   integration::I<T>(SolveOrder), beta, h);
        break;
    case 5:
        incr_convolution_omp<T, herm_matrix<T>, 5>(omp_num_threads, tstp, CPLX(1, 0), C, A,
                                                   Acc, ft.ptr(-1), fttemp, B, Bcc,
                                                   integration::I<T>(SolveOrder), beta, h);
        break;
    case 6:
        incr_convolution_omp<T, herm_matrix<T>, 6>(omp_num_threads, tstp, CPLX(1, 0), C, A,
                                                   Acc, ft.ptr(-1), fttemp, B, Bcc,
                                                   integration::I<T>(SolveOrder), beta, h);
        break;
    case 8:
        incr_convolution_omp<T, herm_matrix<T>, 8>(omp_num_threads, tstp, CPLX(1, 0), C, A,
                                                   Acc, ft.ptr(-1), fttemp, B, Bcc,
                                                   integration::I<T>(SolveOrder), beta, h);
        break;
    default:
        incr_convolution_omp<T, herm_matrix<T>, LARGESIZE>(
            omp_num_threads, tstp, CPLX(1, 0), C, A, Acc, ft.ptr(-1), fttemp, B, Bcc,
            integration::I<T>(SolveOrder), beta, h);
        break;
    }
}