18 #ifndef __ESCRIPT_DATAMATHS_H__
19 #define __ESCRIPT_DATAMATHS_H__
100 template<
typename VEC>
105 typename VEC::size_type inOffset,
108 typename VEC::size_type evOffset)
114 for (i0=0; i0<s0; i0++) {
115 for (i1=0; i1<s1; i1++) {
126 for (i0=0; i0<s0; i0++) {
127 for (i1=0; i1<s1; i1++) {
128 for (i2=0; i2<s2; i2++) {
129 for (i3=0; i3<s3; i3++) {
130 ev[evOffset+
DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = (in[inOffset+
DataTypes::getRelIndex(inShape,i0,i1,i2,i3)] + in[inOffset+
DataTypes::getRelIndex(inShape,i2,i3,i0,i1)]) / 2.0;
149 template<
typename VEC>
154 typename VEC::size_type inOffset,
157 typename VEC::size_type evOffset)
163 for (i0=0; i0<s0; i0++) {
164 for (i1=0; i1<s1; i1++) {
175 for (i0=0; i0<s0; i0++) {
176 for (i1=0; i1<s1; i1++) {
177 for (i2=0; i2<s2; i2++) {
178 for (i3=0; i3<s3; i3++) {
179 ev[evOffset+
DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = (in[inOffset+
DataTypes::getRelIndex(inShape,i0,i1,i2,i3)] - in[inOffset+
DataTypes::getRelIndex(inShape,i2,i3,i0,i1)]) / 2.0;
244 typename VEC::size_type inOffset,
247 typename VEC::size_type evOffset,
257 for (i=0; i<s0; i++) {
262 if (axis_offset==0) {
266 for (i0=0; i0<s0; i0++) {
267 for (i2=0; i2<s2; i2++) {
272 else if (axis_offset==1) {
276 for (i0=0; i0<s0; i0++) {
277 for (i1=0; i1<s1; i1++) {
284 if (axis_offset==0) {
289 for (i0=0; i0<s0; i0++) {
290 for (i2=0; i2<s2; i2++) {
291 for (i3=0; i3<s3; i3++) {
297 else if (axis_offset==1) {
302 for (i0=0; i0<s0; i0++) {
303 for (i1=0; i1<s1; i1++) {
304 for (i3=0; i3<s3; i3++) {
310 else if (axis_offset==2) {
315 for (i0=0; i0<s0; i0++) {
316 for (i1=0; i1<s1; i1++) {
317 for (i2=0; i2<s2; i2++) {
345 typename VEC::size_type inOffset,
348 typename VEC::size_type evOffset,
358 if (axis_offset==1) {
359 for (i0=0; i0<s0; i0++) {
360 for (i1=0; i1<s1; i1++) {
361 for (i2=0; i2<s2; i2++) {
362 for (i3=0; i3<s3; i3++) {
369 else if (axis_offset==2) {
370 for (i0=0; i0<s0; i0++) {
371 for (i1=0; i1<s1; i1++) {
372 for (i2=0; i2<s2; i2++) {
373 for (i3=0; i3<s3; i3++) {
380 else if (axis_offset==3) {
381 for (i0=0; i0<s0; i0++) {
382 for (i1=0; i1<s1; i1++) {
383 for (i2=0; i2<s2; i2++) {
384 for (i3=0; i3<s3; i3++) {
392 for (i0=0; i0<s0; i0++) {
393 for (i1=0; i1<s1; i1++) {
394 for (i2=0; i2<s2; i2++) {
395 for (i3=0; i3<s3; i3++) {
403 else if (inRank == 3) {
408 if (axis_offset==1) {
409 for (i0=0; i0<s0; i0++) {
410 for (i1=0; i1<s1; i1++) {
411 for (i2=0; i2<s2; i2++) {
417 else if (axis_offset==2) {
418 for (i0=0; i0<s0; i0++) {
419 for (i1=0; i1<s1; i1++) {
420 for (i2=0; i2<s2; i2++) {
428 for (i0=0; i0<s0; i0++) {
429 for (i1=0; i1<s1; i1++) {
430 for (i2=0; i2<s2; i2++) {
437 else if (inRank == 2) {
441 if (axis_offset==1) {
442 for (i0=0; i0<s0; i0++) {
443 for (i1=0; i1<s1; i1++) {
449 for (i0=0; i0<s0; i0++) {
450 for (i1=0; i1<s1; i1++) {
456 else if (inRank == 1) {
459 for (i0=0; i0<s0; i0++) {
463 else if (inRank == 0) {
464 ev[evOffset] = in[inOffset];
467 throw DataException(
"Error - DataArrayView::transpose can only be calculated for rank 0, 1, 2, 3 or 4 objects.");
490 typename VEC::size_type inOffset,
493 typename VEC::size_type evOffset,
506 for (i0=0; i0<s0; i0++) {
507 for (i1=0; i1<s1; i1++) {
508 for (i2=0; i2<s2; i2++) {
509 for (i3=0; i3<s3; i3++) {
515 }
else if (axis1==2) {
516 for (i0=0; i0<s0; i0++) {
517 for (i1=0; i1<s1; i1++) {
518 for (i2=0; i2<s2; i2++) {
519 for (i3=0; i3<s3; i3++) {
526 }
else if (axis1==3) {
527 for (i0=0; i0<s0; i0++) {
528 for (i1=0; i1<s1; i1++) {
529 for (i2=0; i2<s2; i2++) {
530 for (i3=0; i3<s3; i3++) {
537 }
else if (axis0==1) {
539 for (i0=0; i0<s0; i0++) {
540 for (i1=0; i1<s1; i1++) {
541 for (i2=0; i2<s2; i2++) {
542 for (i3=0; i3<s3; i3++) {
548 }
else if (axis1==3) {
549 for (i0=0; i0<s0; i0++) {
550 for (i1=0; i1<s1; i1++) {
551 for (i2=0; i2<s2; i2++) {
552 for (i3=0; i3<s3; i3++) {
559 }
else if (axis0==2) {
561 for (i0=0; i0<s0; i0++) {
562 for (i1=0; i1<s1; i1++) {
563 for (i2=0; i2<s2; i2++) {
564 for (i3=0; i3<s3; i3++) {
573 }
else if ( inRank == 3) {
580 for (i0=0; i0<s0; i0++) {
581 for (i1=0; i1<s1; i1++) {
582 for (i2=0; i2<s2; i2++) {
587 }
else if (axis1==2) {
588 for (i0=0; i0<s0; i0++) {
589 for (i1=0; i1<s1; i1++) {
590 for (i2=0; i2<s2; i2++) {
596 }
else if (axis0==1) {
598 for (i0=0; i0<s0; i0++) {
599 for (i1=0; i1<s1; i1++) {
600 for (i2=0; i2<s2; i2++) {
607 }
else if ( inRank == 2) {
613 for (i0=0; i0<s0; i0++) {
614 for (i1=0; i1<s1; i1++) {
621 throw DataException(
"Error - DataArrayView::swapaxes can only be calculated for rank 2, 3 or 4 objects.");
673 eigenvalues3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
691 #pragma clang diagnostic push
692 #pragma clang diagnostic ignored "-Wunused-variable"
695 #pragma clang diagnostic pop
756 const double tol=1.e-13)
758 double in00,in10,in20,in01,in11,in21,in02,in12,in22;
759 double V00,V10,V20,V01,V11,V21,V02,V12,V22;
773 &ev0,&ev1,&V00,&V10,&V01,&V11,tol);
792 &V00,&V10,&V20,&V01,&V11,&V21,&V02,&V12,&V22,tol);
819 typename VEC::size_type offset)
827 template <
class ResVEC,
class LVEC,
class RSCALAR>
830 typename ResVEC::size_type resOffset,
831 const typename ResVEC::size_type samplesToProcess,
832 const typename ResVEC::size_type sampleSize,
834 typename LVEC::size_type leftOffset,
835 const RSCALAR* right,
836 const bool rightreset,
838 bool singleleftsample)
840 size_t substep=(rightreset?0:1);
845 #pragma omp parallel for
846 for (
typename ResVEC::size_type i=0;i<samplesToProcess;++i)
848 typename LVEC::size_type leftbase=leftOffset+(singleleftsample?0:i*sampleSize);
849 const RSCALAR* rpos=right+(rightreset?0:i*substep);
851 for (
typename ResVEC::size_type j=0;j<sampleSize;++j)
853 res[i*sampleSize+resOffset+j]=left[leftbase+j]+*rpos;
860 #pragma omp parallel for
861 for (
typename ResVEC::size_type i=0;i<samplesToProcess;++i)
863 typename LVEC::size_type leftbase=leftOffset+(singleleftsample?0:i*sampleSize);
864 const RSCALAR* rpos=right+(rightreset?0:i*substep);
866 for (
typename ResVEC::size_type j=0;j<sampleSize;++j)
868 res[i*sampleSize+resOffset+j]=pow(left[leftbase+j],*rpos);
875 #pragma omp parallel for
876 for (
typename ResVEC::size_type i=0;i<samplesToProcess;++i)
878 typename LVEC::size_type leftbase=leftOffset+(singleleftsample?0:i*sampleSize);
879 const RSCALAR* rpos=right+(rightreset?0:i*substep);
881 for (
typename ResVEC::size_type j=0;j<sampleSize;++j)
883 res[i*sampleSize+resOffset+j]=left[leftbase+j]-*rpos;
890 #pragma omp parallel for
891 for (
typename ResVEC::size_type i=0;i<samplesToProcess;++i)
893 typename LVEC::size_type leftbase=leftOffset+(singleleftsample?0:i*sampleSize);
894 const RSCALAR* rpos=right+(rightreset?0:i*substep);
896 for (
typename ResVEC::size_type j=0;j<sampleSize;++j)
898 res[i*sampleSize+resOffset+j]=left[leftbase+j] * *rpos;
905 #pragma omp parallel for
906 for (
typename ResVEC::size_type i=0;i<samplesToProcess;++i)
908 typename LVEC::size_type leftbase=leftOffset+(singleleftsample?0:i*sampleSize);
909 const RSCALAR* rpos=right+(rightreset?0:i*substep);
911 for (
typename ResVEC::size_type j=0;j<sampleSize;++j)
913 res[i*sampleSize+resOffset+j]=left[leftbase+j]/ *rpos;
932 const bool rightreset,
934 bool singleleftsample);
939 template <
class ResVEC,
class LSCALAR,
class RVEC>
942 typename ResVEC::size_type resOffset,
943 const typename ResVEC::size_type samplesToProcess,
944 const typename ResVEC::size_type sampleSize,
946 const bool leftreset,
948 typename RVEC::size_type rightOffset,
950 bool singlerightsample)
952 size_t substep=(leftreset?0:1);
957 #pragma omp parallel for
958 for (
typename ResVEC::size_type i=0;i<samplesToProcess;++i)
960 typename RVEC::size_type rightbase=rightOffset+(singlerightsample?0:i*sampleSize);
961 const LSCALAR* lpos=left+(leftreset?0:i*substep);
962 for (
typename ResVEC::size_type j=0;j<sampleSize;++j)
964 res[i*sampleSize+resOffset+j]=*lpos+right[rightbase+j];
971 #pragma omp parallel for
972 for (
typename ResVEC::size_type i=0;i<samplesToProcess;++i)
974 typename RVEC::size_type rightbase=rightOffset+(singlerightsample?0:i*sampleSize);
975 const LSCALAR* lpos=left+(leftreset?0:i*substep);
976 for (
typename ResVEC::size_type j=0;j<sampleSize;++j)
978 res[i*sampleSize+resOffset+j]=pow(*lpos,right[rightbase+j]);
985 #pragma omp parallel for
986 for (
typename ResVEC::size_type i=0;i<samplesToProcess;++i)
988 typename RVEC::size_type rightbase=rightOffset+(singlerightsample?0:i*sampleSize);
989 const LSCALAR* lpos=left+(leftreset?0:i*substep);
990 for (
typename ResVEC::size_type j=0;j<sampleSize;++j)
992 res[i*sampleSize+resOffset+j]=*lpos-right[rightbase+j];
999 #pragma omp parallel for
1000 for (
typename ResVEC::size_type i=0;i<samplesToProcess;++i)
1002 typename RVEC::size_type rightbase=rightOffset+(singlerightsample?0:i*sampleSize);
1003 const LSCALAR* lpos=left+(leftreset?0:i*substep);
1004 for (
typename ResVEC::size_type j=0;j<sampleSize;++j)
1006 res[i*sampleSize+resOffset+j]=*lpos*right[rightbase+j];
1013 #pragma omp parallel for
1014 for (
typename ResVEC::size_type i=0;i<samplesToProcess;++i)
1016 typename RVEC::size_type rightbase=rightOffset+(singlerightsample?0:i*sampleSize);
1017 const LSCALAR* lpos=left+(leftreset?0:i*substep);
1018 for (
typename ResVEC::size_type j=0;j<sampleSize;++j)
1020 res[i*sampleSize+resOffset+j]=*lpos/right[rightbase+j];
1037 const bool leftreset,
1041 bool singlerightsample);
1046 template <
class ResVEC,
class LVEC,
class RVEC>
1049 typename ResVEC::size_type resOffset,
1050 const typename ResVEC::size_type samplesToProcess,
1051 const typename ResVEC::size_type sampleSize,
1053 typename LVEC::size_type leftOffset,
1054 const bool leftreset,
1056 typename RVEC::size_type rightOffset,
1057 const bool rightreset,
1064 #pragma omp parallel for
1065 for (
typename ResVEC::size_type i=0;i<samplesToProcess;++i)
1067 typename LVEC::size_type leftbase=leftOffset+(leftreset?0:i*sampleSize);
1068 typename RVEC::size_type rightbase=rightOffset+(rightreset?0:i*sampleSize);
1069 for (
typename ResVEC::size_type j=0;j<sampleSize;++j)
1071 res[i*sampleSize+resOffset+j]=left[leftbase+j]+right[rightbase+j];
1078 #pragma omp parallel for
1079 for (
typename ResVEC::size_type i=0;i<samplesToProcess;++i)
1081 typename LVEC::size_type leftbase=leftOffset+(leftreset?0:i*sampleSize);
1082 typename RVEC::size_type rightbase=rightOffset+(rightreset?0:i*sampleSize);
1083 for (
typename ResVEC::size_type j=0;j<sampleSize;++j)
1085 res[i*sampleSize+resOffset+j]=pow(left[leftbase+j],right[rightbase+j]);
1092 #pragma omp parallel for
1093 for (
typename ResVEC::size_type i=0;i<samplesToProcess;++i)
1095 typename LVEC::size_type leftbase=leftOffset+(leftreset?0:i*sampleSize);
1096 typename RVEC::size_type rightbase=rightOffset+(rightreset?0:i*sampleSize);
1097 for (
typename ResVEC::size_type j=0;j<sampleSize;++j)
1099 res[i*sampleSize+resOffset+j]=left[leftbase+j]-right[rightbase+j];
1106 #pragma omp parallel for
1107 for (
typename ResVEC::size_type i=0;i<samplesToProcess;++i)
1109 typename LVEC::size_type leftbase=leftOffset+(leftreset?0:i*sampleSize);
1110 typename RVEC::size_type rightbase=rightOffset+(rightreset?0:i*sampleSize);
1111 for (
typename ResVEC::size_type j=0;j<sampleSize;++j)
1113 res[i*sampleSize+resOffset+j]=left[leftbase+j]*right[rightbase+j];
1120 #pragma omp parallel for
1121 for (
typename ResVEC::size_type i=0;i<samplesToProcess;++i)
1123 typename LVEC::size_type leftbase=leftOffset+(leftreset?0:i*sampleSize);
1124 typename RVEC::size_type rightbase=rightOffset+(rightreset?0:i*sampleSize);
1125 for (
typename ResVEC::size_type j=0;j<sampleSize;++j)
1127 res[i*sampleSize+resOffset+j]=left[leftbase+j]/right[rightbase+j];
1145 const bool leftreset,
1148 const bool rightreset,
1151 #define OPVECLAZYBODY(X) \
1152 for (size_t j=0;j<onumsteps;++j)\
1154 for (size_t i=0;i<numsteps;++i,res+=resultStep) \
1156 for (size_t s=0; s<chunksize; ++s)\
1161 lroffset+=leftstep; \
1162 rroffset+=rightstep; \
1164 lroffset+=oleftstep;\
1165 rroffset+=orightstep;\
1174 template <
class ResELT,
class LELT,
class RELT>
1179 const size_t chunksize,
1180 const size_t onumsteps,
1181 const size_t numsteps,
1182 const size_t resultStep,
1183 const size_t leftstep,
1184 const size_t rightstep,
1185 const size_t oleftstep,
1186 const size_t orightstep,
1209 ESYS_ASSERT(
false,
"Invalid operation. This should never happen!");
1221 template <
class ResELT,
class LELT,
class RELT>
1226 const size_t chunksize,
1227 const size_t onumsteps,
1228 const size_t numsteps,
1229 const size_t resultStep,
1230 const size_t leftstep,
1231 const size_t rightstep,
1232 const size_t oleftstep,
1233 const size_t orightstep,
1253 ESYS_ASSERT(
false,
"Invalid operation. This should never happen!");
1262 template <
class ResVEC,
class LVEC,
class RVEC>
1265 const typename ResVEC::size_type samplesToProcess,
1266 const typename ResVEC::size_type DPPSample,
1267 const typename ResVEC::size_type DPSize,
1276 typename ResVEC::size_type lstep=leftscalar?1:DPSize;
1277 typename ResVEC::size_type rstep=rightscalar?1:DPSize;
1278 typename ResVEC::size_type limit=samplesToProcess*DPPSample;
1283 #pragma omp parallel for
1284 for (
typename ResVEC::size_type i=0;i<limit;++i)
1286 typename LVEC::size_type leftbase=(lefttagged?tagsource.
getPointOffset(i/DPPSample,0):i*lstep);
1287 typename RVEC::size_type rightbase=(lefttagged?i*rstep:tagsource.
getPointOffset(i/DPPSample,0));
1289 for (
typename ResVEC::size_type j=0;j<DPSize;++j)
1291 res[i*DPSize+j]=left[leftbase+j*(!leftscalar)]+right[rightbase+j*(!rightscalar)];
1298 #pragma omp parallel for
1299 for (
typename ResVEC::size_type i=0;i<limit;++i)
1301 typename LVEC::size_type leftbase=(lefttagged?tagsource.
getPointOffset(i/DPPSample,0):i*lstep);
1302 typename RVEC::size_type rightbase=(lefttagged?i*rstep:tagsource.
getPointOffset(i/DPPSample,0));
1304 for (
typename ResVEC::size_type j=0;j<DPSize;++j)
1306 res[i*DPSize+j]=pow(left[leftbase+j*(!leftscalar)],right[rightbase+j*(!rightscalar)]);
1313 #pragma omp parallel for
1314 for (
typename ResVEC::size_type i=0;i<limit;++i)
1316 typename LVEC::size_type leftbase=(lefttagged?tagsource.
getPointOffset(i/DPPSample,0):i*lstep);
1317 typename RVEC::size_type rightbase=(lefttagged?i*rstep:tagsource.
getPointOffset(i/DPPSample,0));
1319 for (
typename ResVEC::size_type j=0;j<DPSize;++j)
1321 res[i*DPSize+j]=left[leftbase+j*(!leftscalar)]-right[rightbase+j*(!rightscalar)];
1328 #pragma omp parallel for
1329 for (
typename ResVEC::size_type i=0;i<limit;++i)
1331 typename LVEC::size_type leftbase=(lefttagged?tagsource.
getPointOffset(i/DPPSample,0):i*lstep);
1332 typename RVEC::size_type rightbase=(lefttagged?i*rstep:tagsource.
getPointOffset(i/DPPSample,0));
1334 for (
typename ResVEC::size_type j=0;j<DPSize;++j)
1336 res[i*DPSize+j]=left[leftbase+j*(!leftscalar)]*right[rightbase+j*(!rightscalar)];
1343 #pragma omp parallel for
1344 for (
typename ResVEC::size_type i=0;i<limit;++i)
1346 typename LVEC::size_type leftbase=(lefttagged?tagsource.
getPointOffset(i/DPPSample,0):i*lstep);
1347 typename RVEC::size_type rightbase=(lefttagged?i*rstep:tagsource.
getPointOffset(i/DPPSample,0));
1349 for (
typename ResVEC::size_type j=0;j<DPSize;++j)
1351 res[i*DPSize+j]=left[leftbase+j*(!leftscalar)]/right[rightbase+j*(!rightscalar)];
1368 const bool leftscalar,
1370 const bool rightscalar,
1371 const bool lefttagged,
1372 const DataTagged& tagsource,
1394 template <
class BinaryFunction>
1400 BinaryFunction operation,
1404 "Couldn't perform reductionOp due to insufficient storage.");
1407 current_value=operation(current_value,left[offset+i]);
1409 return current_value;
1412 template <
class BinaryFunction>
1418 BinaryFunction operation,
1422 "Couldn't perform reductionOp due to insufficient storage.");
1425 current_value=operation(current_value,left[offset+i]);
1427 return current_value;
1455 LapackInverseHelper& helper);
1475 for (
size_t z=inOffset;z<inOffset+count;++z)
1489 for (
size_t z=inOffset;z<inOffset+count;++z)
1502 #endif // __ESCRIPT_DATAMATHS_H__