36#ifndef _FOURTHORDERTENSOR_HPP_
37#define _FOURTHORDERTENSOR_HPP_
50template<
unsigned DIM1,
unsigned DIM2,
unsigned DIM3,
unsigned DIM4>
69 return M + DIM1*N + DIM1*DIM2*P + DIM1*DIM2*DIM3*Q;
87 template<
unsigned CONTRACTED_DIM>
99 template<
unsigned CONTRACTED_DIM>
110 template<
unsigned CONTRACTED_DIM>
121 template<
unsigned CONTRACTED_DIM>
132 double&
operator()(
unsigned M,
unsigned N,
unsigned P,
unsigned Q);
152template<
unsigned DIM1,
unsigned DIM2,
unsigned DIM3,
unsigned DIM4>
155 unsigned size = DIM1*DIM2*DIM3*DIM4;
156 mData.resize(size, 0.0);
159template<
unsigned DIM1,
unsigned DIM2,
unsigned DIM3,
unsigned DIM4>
160template<
unsigned CONTRACTED_DIM>
165 std::vector<double>::iterator iter = mData.begin();
166 std::vector<double>::iterator other_tensor_iter = rTensor.
rGetData().begin();
168 for (
unsigned d=0; d<DIM4; d++)
170 for (
unsigned c=0; c<DIM3; c++)
172 for (
unsigned b=0; b<DIM2; b++)
174 for (
unsigned a=0; a<DIM1; a++)
176 for (
unsigned N=0; N<CONTRACTED_DIM; N++)
186 *iter += rMatrix(a,N) * *other_tensor_iter;
194 other_tensor_iter -= CONTRACTED_DIM;
202template<
unsigned DIM1,
unsigned DIM2,
unsigned DIM3,
unsigned DIM4>
203template<
unsigned CONTRACTED_DIM>
208 std::vector<double>::iterator iter = mData.begin();
209 std::vector<double>::iterator other_tensor_iter = rTensor.
rGetData().begin();
211 for (
unsigned d=0; d<DIM4; d++)
213 for (
unsigned c=0; c<DIM3; c++)
215 for (
unsigned b=0; b<DIM2; b++)
217 for (
unsigned N=0; N<CONTRACTED_DIM; N++)
219 for (
unsigned a=0; a<DIM1; a++)
229 *iter += rMatrix(b,N) * *other_tensor_iter;
234 if (N != CONTRACTED_DIM-1)
241 other_tensor_iter -= CONTRACTED_DIM*DIM1;
248template<
unsigned DIM1,
unsigned DIM2,
unsigned DIM3,
unsigned DIM4>
249template<
unsigned CONTRACTED_DIM>
254 std::vector<double>::iterator iter = mData.begin();
255 std::vector<double>::iterator other_tensor_iter = rTensor.
rGetData().begin();
257 for (
unsigned d=0; d<DIM4; d++)
259 for (
unsigned c=0; c<DIM3; c++)
261 for (
unsigned N=0; N<CONTRACTED_DIM; N++)
263 for (
unsigned b=0; b<DIM2; b++)
265 for (
unsigned a=0; a<DIM1; a++)
275 *iter += rMatrix(c,N) * *other_tensor_iter;
281 if (N != CONTRACTED_DIM-1)
289 other_tensor_iter -= CONTRACTED_DIM*DIM1*DIM2;
295template<
unsigned DIM1,
unsigned DIM2,
unsigned DIM3,
unsigned DIM4>
296template<
unsigned CONTRACTED_DIM>
301 std::vector<double>::iterator iter = mData.begin();
302 std::vector<double>::iterator other_tensor_iter = rTensor.
rGetData().begin();
304 for (
unsigned d=0; d<DIM4; d++)
306 for (
unsigned N=0; N<CONTRACTED_DIM; N++)
308 for (
unsigned c=0; c<DIM3; c++)
310 for (
unsigned b=0; b<DIM2; b++)
312 for (
unsigned a=0; a<DIM1; a++)
322 *iter += rMatrix(d,N) * *other_tensor_iter;
330 if (N != CONTRACTED_DIM-1)
332 iter-= DIM1*DIM2*DIM3;
336 other_tensor_iter -= CONTRACTED_DIM*DIM1*DIM2*DIM3;
340template<
unsigned DIM1,
unsigned DIM2,
unsigned DIM3,
unsigned DIM4>
348 return mData[GetVectorIndex(M,N,P,Q)];
351template<
unsigned DIM1,
unsigned DIM2,
unsigned DIM3,
unsigned DIM4>
354 for (
unsigned i=0; i<mData.size(); i++)
double & operator()(unsigned M, unsigned N, unsigned P, unsigned Q)
std::vector< double > & rGetData()
std::vector< double > mData
unsigned GetVectorIndex(unsigned M, unsigned N, unsigned P, unsigned Q)
void SetAsContractionOnSecondDimension(const c_matrix< double, DIM2, CONTRACTED_DIM > &rMatrix, FourthOrderTensor< DIM1, CONTRACTED_DIM, DIM3, DIM4 > &rTensor)
void SetAsContractionOnFourthDimension(const c_matrix< double, DIM4, CONTRACTED_DIM > &rMatrix, FourthOrderTensor< DIM1, DIM2, DIM3, CONTRACTED_DIM > &rTensor)
void SetAsContractionOnFirstDimension(const c_matrix< double, DIM1, CONTRACTED_DIM > &rMatrix, FourthOrderTensor< CONTRACTED_DIM, DIM2, DIM3, DIM4 > &rTensor)
void SetAsContractionOnThirdDimension(const c_matrix< double, DIM3, CONTRACTED_DIM > &rMatrix, FourthOrderTensor< DIM1, DIM2, CONTRACTED_DIM, DIM4 > &rTensor)