NESSi  v1.0.2
The NonEquilibrium Systems Simulation Library

◆ response_convolution() [4/4]

template<typename T >
void cntr::response_convolution ( int  tstp,
CPLX &  cc,
GREEN &  W,
int  a1,
int  a2,
function< T > &  f,
int  b1,
int  b2,
int  SolveOrder,
beta,
h 
)

Definition at line 96 of file cntr_response_convolution_impl.hpp.

97  {
98  int ntau = W.ntau();
99  T dtau = beta / ntau;
100  int sizew = W.size1(), idxw = a1 * sizew + a2, sw = sizew * sizew;
101  int sizef = f.size1(), idxf = b1 * sizef + b2, sf = sizef * sizef;
102  int n1 = (tstp == -1 || tstp >= SolveOrder ? SolveOrder : SolveOrder);
103  CPLX res;
104 
105  assert(W.nt() >= tstp);
106  assert(tstp == -1 || tstp >= SolveOrder);
107  assert(SolveOrder <= ntau);
108  assert(b1 >= 0 && b1 <= sizef-1);
109  assert(b2 >= 0 && b2 <= sizef-1);
110  assert(a1 >= 0 && a1 <= sizew-1);
111  assert(a2 >= 0 && a2 <= sizew-1);
112 
113  if (tstp == -1) {
114  response_integrate<T>(ntau, dtau, res, W.matptr(0), sw, idxw, integration::I<T>(SolveOrder));
115  cc = res * f.ptr(-1)[idxf];
116  } else if (tstp >= SolveOrder) {
117  response_integrate<T>(ntau, dtau, res, W.tvptr(tstp, 0), sw, idxw,
118  integration::I<T>(SolveOrder));
119  cc = res * (CPLX(0, -1.0) * f.ptr(-1)[idxf]);
120  response_integrate<T>(tstp, h, res, W.retptr(tstp, 0), sw, idxw, f.ptr(0), sf, idxf,
121  integration::I<T>(SolveOrder));
122  cc += res;
123  } else {
124  response_integrate<T>(ntau, dtau, res, W.tvptr(tstp, 0), sw, idxw,
125  integration::I<T>(SolveOrder));
126  cc = res * (CPLX(0, -1.0) * f.ptr(-1)[idxf]);
127  // need to extrapolate W(t,t') for t'>t
128  // ** USING THE ASSUMPTION THAT W IS HERMITIAN**
129  // i.e., Wret(t,t') = -Wret(t',t)^*
130  CPLX *wtmp = new CPLX[SolveOrder + 1];
131  for (int i = 0; i <= tstp; i++)
132  wtmp[i] = W.retptr(tstp, i)[a1 * sizew + a2];
133  for (int i = tstp + 1; i <= SolveOrder; i++)
134  wtmp[i] = -conj(W.retptr(i, tstp)[a2 * sizew + a1]);
135  response_integrate<T>(tstp, h, res, wtmp, 1, 0, f.ptr(0), sf, idxf,
136  integration::I<T>(SolveOrder));
137  cc += res;
138  delete[] wtmp;
139  }
140 }