escript  Revision_
DataVectorOps.h
Go to the documentation of this file.
1 
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2020 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Apache License, version 2.0
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014-2017 by Centre for Geoscience Computing (GeoComp)
14 * Development from 2019 by School of Earth and Environmental Sciences
15 **
16 *****************************************************************************/
17 
18 #ifndef __ESCRIPT_DATAMATHS_H__
19 #define __ESCRIPT_DATAMATHS_H__
20 
21 #include "DataAbstract.h"
22 #include "DataException.h"
23 #include "ArrayOps.h"
24 #include "LapackInverseHelper.h"
25 #include "DataTagged.h"
26 #include <complex>
37 namespace escript
38 {
39 
61  void
63  const DataTypes::ShapeType& leftShape,
65  const DataTypes::RealVectorType& right,
66  const DataTypes::ShapeType& rightShape,
69  const DataTypes::ShapeType& resultShape);
70 // Hmmmm why is there no offset for the result??
71 
72 
73 
74 
86  const DataTypes::ShapeType& right);
87 
88 
100  template<typename VEC>
101  inline
102  void
103  symmetric(const VEC& in,
104  const DataTypes::ShapeType& inShape,
105  typename VEC::size_type inOffset,
106  VEC& ev,
107  const DataTypes::ShapeType& evShape,
108  typename VEC::size_type evOffset)
109  {
110  if (DataTypes::getRank(inShape) == 2) {
111  int i0, i1;
112  int s0=inShape[0];
113  int s1=inShape[1];
114  for (i0=0; i0<s0; i0++) {
115  for (i1=0; i1<s1; i1++) {
116  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = (in[inOffset+DataTypes::getRelIndex(inShape,i0,i1)] + in[inOffset+DataTypes::getRelIndex(inShape,i1,i0)]) / 2.0;
117  }
118  }
119  }
120  else if (DataTypes::getRank(inShape) == 4) {
121  int i0, i1, i2, i3;
122  int s0=inShape[0];
123  int s1=inShape[1];
124  int s2=inShape[2];
125  int s3=inShape[3];
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;
131  }
132  }
133  }
134  }
135  }
136  }
137 
149  template<typename VEC>
150  inline
151  void
152  antisymmetric(const VEC& in,
153  const DataTypes::ShapeType& inShape,
154  typename VEC::size_type inOffset,
155  VEC& ev,
156  const DataTypes::ShapeType& evShape,
157  typename VEC::size_type evOffset)
158  {
159  if (DataTypes::getRank(inShape) == 2) {
160  int i0, i1;
161  int s0=inShape[0];
162  int s1=inShape[1];
163  for (i0=0; i0<s0; i0++) {
164  for (i1=0; i1<s1; i1++) {
165  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = (in[inOffset+DataTypes::getRelIndex(inShape,i0,i1)] - in[inOffset+DataTypes::getRelIndex(inShape,i1,i0)]) / 2.0;
166  }
167  }
168  }
169  else if (DataTypes::getRank(inShape) == 4) {
170  int i0, i1, i2, i3;
171  int s0=inShape[0];
172  int s1=inShape[1];
173  int s2=inShape[2];
174  int s3=inShape[3];
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;
180  }
181  }
182  }
183  }
184  }
185  }
186 
187 
188 
200  void
202  const DataTypes::ShapeType& inShape,
205  const DataTypes::ShapeType& evShape,
207 
219  void
221  const DataTypes::ShapeType& inShape,
222  typename DataTypes::CplxVectorType::size_type inOffset,
224  const DataTypes::ShapeType& evShape,
225  typename DataTypes::CplxVectorType::size_type evOffset);
226 
239  template <class VEC>
240  inline
241  void
242  trace(const VEC& in,
243  const DataTypes::ShapeType& inShape,
244  typename VEC::size_type inOffset,
245  VEC& ev,
246  const DataTypes::ShapeType& evShape,
247  typename VEC::size_type evOffset,
248  int axis_offset)
249  {
250  for (int j=0;j<DataTypes::noValues(evShape);++j)
251  {
252  ev[evOffset+j]=0;
253  }
254  if (DataTypes::getRank(inShape) == 2) {
255  int s0=inShape[0]; // Python wrapper limits to square matrix
256  int i;
257  for (i=0; i<s0; i++) {
258  ev[evOffset/*+DataTypes::getRelIndex(evShape)*/] += in[inOffset+DataTypes::getRelIndex(inShape,i,i)];
259  }
260  }
261  else if (DataTypes::getRank(inShape) == 3) {
262  if (axis_offset==0) {
263  int s0=inShape[0];
264  int s2=inShape[2];
265  int i0, i2;
266  for (i0=0; i0<s0; i0++) {
267  for (i2=0; i2<s2; i2++) {
268  ev[evOffset+DataTypes::getRelIndex(evShape,i2)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i0,i2)];
269  }
270  }
271  }
272  else if (axis_offset==1) {
273  int s0=inShape[0];
274  int s1=inShape[1];
275  int i0, i1;
276  for (i0=0; i0<s0; i0++) {
277  for (i1=0; i1<s1; i1++) {
278  ev[evOffset+DataTypes::getRelIndex(evShape,i0)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i1)];
279  }
280  }
281  }
282  }
283  else if (DataTypes::getRank(inShape) == 4) {
284  if (axis_offset==0) {
285  int s0=inShape[0];
286  int s2=inShape[2];
287  int s3=inShape[3];
288  int i0, i2, i3;
289  for (i0=0; i0<s0; i0++) {
290  for (i2=0; i2<s2; i2++) {
291  for (i3=0; i3<s3; i3++) {
292  ev[evOffset+DataTypes::getRelIndex(evShape,i2,i3)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i0,i2,i3)];
293  }
294  }
295  }
296  }
297  else if (axis_offset==1) {
298  int s0=inShape[0];
299  int s1=inShape[1];
300  int s3=inShape[3];
301  int i0, i1, i3;
302  for (i0=0; i0<s0; i0++) {
303  for (i1=0; i1<s1; i1++) {
304  for (i3=0; i3<s3; i3++) {
305  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i3)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i1,i3)];
306  }
307  }
308  }
309  }
310  else if (axis_offset==2) {
311  int s0=inShape[0];
312  int s1=inShape[1];
313  int s2=inShape[2];
314  int i0, i1, i2;
315  for (i0=0; i0<s0; i0++) {
316  for (i1=0; i1<s1; i1++) {
317  for (i2=0; i2<s2; i2++) {
318  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i2,i2)];
319  }
320  }
321  }
322  }
323  }
324  }
325 
326 
339  // ESCRIPT_DLL_API can only export template specialisation!
340  template <class VEC>
341  inline
342  void
343  transpose(const VEC& in,
344  const DataTypes::ShapeType& inShape,
345  typename VEC::size_type inOffset,
346  VEC& ev,
347  const DataTypes::ShapeType& evShape,
348  typename VEC::size_type evOffset,
349  int axis_offset)
350  {
351  int inRank=DataTypes::getRank(inShape);
352  if ( inRank== 4) {
353  int s0=evShape[0];
354  int s1=evShape[1];
355  int s2=evShape[2];
356  int s3=evShape[3];
357  int i0, i1, i2, i3;
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++) {
363  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i3,i0,i1,i2)];
364  }
365  }
366  }
367  }
368  }
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++) {
374  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i3,i0,i1)];
375  }
376  }
377  }
378  }
379  }
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++) {
385  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i2,i3,i0)];
386  }
387  }
388  }
389  }
390  }
391  else {
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++) {
396  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i2,i3)];
397  }
398  }
399  }
400  }
401  }
402  }
403  else if (inRank == 3) {
404  int s0=evShape[0];
405  int s1=evShape[1];
406  int s2=evShape[2];
407  int i0, i1, i2;
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++) {
412  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i0,i1)];
413  }
414  }
415  }
416  }
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++) {
421  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i2,i0)];
422  }
423  }
424  }
425  }
426  else {
427  // Copy the matrix unchanged
428  for (i0=0; i0<s0; i0++) {
429  for (i1=0; i1<s1; i1++) {
430  for (i2=0; i2<s2; i2++) {
431  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i2)];
432  }
433  }
434  }
435  }
436  }
437  else if (inRank == 2) {
438  int s0=evShape[0];
439  int s1=evShape[1];
440  int i0, i1;
441  if (axis_offset==1) {
442  for (i0=0; i0<s0; i0++) {
443  for (i1=0; i1<s1; i1++) {
444  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0)];
445  }
446  }
447  }
448  else {
449  for (i0=0; i0<s0; i0++) {
450  for (i1=0; i1<s1; i1++) {
451  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1)];
452  }
453  }
454  }
455  }
456  else if (inRank == 1) {
457  int s0=evShape[0];
458  int i0;
459  for (i0=0; i0<s0; i0++) {
460  ev[evOffset+DataTypes::getRelIndex(evShape,i0)] = in[inOffset+DataTypes::getRelIndex(inShape,i0)];
461  }
462  }
463  else if (inRank == 0) {
464  ev[evOffset/*+DataTypes::getRelIndex(evShape,)*/] = in[inOffset/*+DataTypes::getRelIndex(inShape,)*/];
465  }
466  else {
467  throw DataException("Error - DataArrayView::transpose can only be calculated for rank 0, 1, 2, 3 or 4 objects.");
468  }
469  }
470 
484  // ESCRIPT_DLL_API can only export template specialisation!
485  template <class VEC>
486  inline
487  void
488  swapaxes(const VEC& in,
489  const DataTypes::ShapeType& inShape,
490  typename VEC::size_type inOffset,
491  VEC& ev,
492  const DataTypes::ShapeType& evShape,
493  typename VEC::size_type evOffset,
494  int axis0,
495  int axis1)
496  {
497  int inRank=DataTypes::getRank(inShape);
498  if (inRank == 4) {
499  int s0=evShape[0];
500  int s1=evShape[1];
501  int s2=evShape[2];
502  int s3=evShape[3];
503  int i0, i1, i2, i3;
504  if (axis0==0) {
505  if (axis1==1) {
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++) {
510  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0,i2,i3)];
511  }
512  }
513  }
514  }
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++) {
520  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i1,i0,i3)];
521  }
522  }
523  }
524  }
525 
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++) {
531  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i3,i1,i2,i0)];
532  }
533  }
534  }
535  }
536  }
537  } else if (axis0==1) {
538  if (axis1==2) {
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++) {
543  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i2,i1,i3)];
544  }
545  }
546  }
547  }
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++) {
553  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i3,i2,i1)];
554  }
555  }
556  }
557  }
558  }
559  } else if (axis0==2) {
560  if (axis1==3) {
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++) {
565  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i3,i2)];
566  }
567  }
568  }
569  }
570  }
571  }
572 
573  } else if ( inRank == 3) {
574  int s0=evShape[0];
575  int s1=evShape[1];
576  int s2=evShape[2];
577  int i0, i1, i2;
578  if (axis0==0) {
579  if (axis1==1) {
580  for (i0=0; i0<s0; i0++) {
581  for (i1=0; i1<s1; i1++) {
582  for (i2=0; i2<s2; i2++) {
583  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0,i2)];
584  }
585  }
586  }
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++) {
591  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i1,i0)];
592  }
593  }
594  }
595  }
596  } else if (axis0==1) {
597  if (axis1==2) {
598  for (i0=0; i0<s0; i0++) {
599  for (i1=0; i1<s1; i1++) {
600  for (i2=0; i2<s2; i2++) {
601  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i2,i1)];
602  }
603  }
604  }
605  }
606  }
607  } else if ( inRank == 2) {
608  int s0=evShape[0];
609  int s1=evShape[1];
610  int i0, i1;
611  if (axis0==0) {
612  if (axis1==1) {
613  for (i0=0; i0<s0; i0++) {
614  for (i1=0; i1<s1; i1++) {
615  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0)];
616  }
617  }
618  }
619  }
620  } else {
621  throw DataException("Error - DataArrayView::swapaxes can only be calculated for rank 2, 3 or 4 objects.");
622  }
623  }
624 
637  inline
638  void
640  const DataTypes::ShapeType& inShape,
641  typename DataTypes::RealVectorType::size_type inOffset,
643  const DataTypes::ShapeType& evShape,
644  typename DataTypes::RealVectorType::size_type evOffset)
645  {
646  typename DataTypes::RealVectorType::ElementType in00,in10,in20,in01,in11,in21,in02,in12,in22;
647  typename DataTypes::RealVectorType::ElementType ev0,ev1,ev2;
648  int s=inShape[0];
649  if (s==1) {
650  in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
651  eigenvalues1(in00,&ev0);
652  ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
653 
654  } else if (s==2) {
655  in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
656  in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
657  in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
658  in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
659  eigenvalues2(in00,(in01+in10)/2.,in11,&ev0,&ev1);
660  ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
661  ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
662 
663  } else if (s==3) {
664  in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
665  in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
666  in20=in[inOffset+DataTypes::getRelIndex(inShape,2,0)];
667  in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
668  in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
669  in21=in[inOffset+DataTypes::getRelIndex(inShape,2,1)];
670  in02=in[inOffset+DataTypes::getRelIndex(inShape,0,2)];
671  in12=in[inOffset+DataTypes::getRelIndex(inShape,1,2)];
672  in22=in[inOffset+DataTypes::getRelIndex(inShape,2,2)];
673  eigenvalues3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
674  &ev0,&ev1,&ev2);
675  ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
676  ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
677  ev[evOffset+DataTypes::getRelIndex(evShape,2)]=ev2;
678 
679  }
680  }
681 
682  inline
683  void
685  const DataTypes::ShapeType& inShape,
686  typename DataTypes::CplxVectorType::size_type inOffset,
688  const DataTypes::ShapeType& evShape,
689  typename DataTypes::CplxVectorType::size_type evOffset)
690  {
691 #pragma clang diagnostic push
692 #pragma clang diagnostic ignored "-Wunused-variable"
693  typename DataTypes::CplxVectorType::ElementType in00,in10,in20,in01,in11,in21,in02,in12,in22;
694  typename DataTypes::CplxVectorType::ElementType ev0,ev1,ev2;
695 #pragma clang diagnostic pop
696  int s=inShape[0];
697  if (s==1) {
698  in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
699  eigenvalues1(in00,&ev0);
700  ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
701 
702  } else if (s==2) {
703  in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
704  in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
705  in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
706  in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
707  eigenvalues2(in00,(in01+in10)/2.,in11,&ev0,&ev1);
708  ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
709  ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
710 
711  } else if (s==3) {
712  // this doesn't work yet
713 // in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
714 // in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
715 // in20=in[inOffset+DataTypes::getRelIndex(inShape,2,0)];
716 // in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
717 // in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
718 // in21=in[inOffset+DataTypes::getRelIndex(inShape,2,1)];
719 // in02=in[inOffset+DataTypes::getRelIndex(inShape,0,2)];
720 // in12=in[inOffset+DataTypes::getRelIndex(inShape,1,2)];
721 // in22=in[inOffset+DataTypes::getRelIndex(inShape,2,2)];
722 // eigenvalues3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
723 // &ev0,&ev1,&ev2);
724 // ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
725 // ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
726 // ev[evOffset+DataTypes::getRelIndex(evShape,2)]=ev2;
727 
728  }
729  }
730 
731 
748  inline
749  void
756  const double tol=1.e-13)
757  {
758  double in00,in10,in20,in01,in11,in21,in02,in12,in22;
759  double V00,V10,V20,V01,V11,V21,V02,V12,V22;
760  double ev0,ev1,ev2;
761  int s=inShape[0];
762  if (s==1) {
763  in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
764  eigenvalues_and_eigenvectors1(in00,&ev0,&V00,tol);
765  ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
766  V[inOffset+DataTypes::getRelIndex(VShape,0,0)]=V00;
767  } else if (s==2) {
768  in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
769  in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
770  in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
771  in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
772  eigenvalues_and_eigenvectors2(in00,(in01+in10)/2.,in11,
773  &ev0,&ev1,&V00,&V10,&V01,&V11,tol);
774  ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
775  ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
776  V[inOffset+DataTypes::getRelIndex(VShape,0,0)]=V00;
777  V[inOffset+DataTypes::getRelIndex(VShape,1,0)]=V10;
778  V[inOffset+DataTypes::getRelIndex(VShape,0,1)]=V01;
779  V[inOffset+DataTypes::getRelIndex(VShape,1,1)]=V11;
780  } else if (s==3) {
781  in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
782  in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
783  in20=in[inOffset+DataTypes::getRelIndex(inShape,2,0)];
784  in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
785  in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
786  in21=in[inOffset+DataTypes::getRelIndex(inShape,2,1)];
787  in02=in[inOffset+DataTypes::getRelIndex(inShape,0,2)];
788  in12=in[inOffset+DataTypes::getRelIndex(inShape,1,2)];
789  in22=in[inOffset+DataTypes::getRelIndex(inShape,2,2)];
790  eigenvalues_and_eigenvectors3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
791  &ev0,&ev1,&ev2,
792  &V00,&V10,&V20,&V01,&V11,&V21,&V02,&V12,&V22,tol);
793  ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
794  ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
795  ev[evOffset+DataTypes::getRelIndex(evShape,2)]=ev2;
796  V[inOffset+DataTypes::getRelIndex(VShape,0,0)]=V00;
797  V[inOffset+DataTypes::getRelIndex(VShape,1,0)]=V10;
798  V[inOffset+DataTypes::getRelIndex(VShape,2,0)]=V20;
799  V[inOffset+DataTypes::getRelIndex(VShape,0,1)]=V01;
800  V[inOffset+DataTypes::getRelIndex(VShape,1,1)]=V11;
801  V[inOffset+DataTypes::getRelIndex(VShape,2,1)]=V21;
802  V[inOffset+DataTypes::getRelIndex(VShape,0,2)]=V02;
803  V[inOffset+DataTypes::getRelIndex(VShape,1,2)]=V12;
804  V[inOffset+DataTypes::getRelIndex(VShape,2,2)]=V22;
805 
806  }
807  }
808 
809 
814 template <class VEC>
815 inline
816 bool
817 checkOffset(const VEC& data,
818  const DataTypes::ShapeType& shape,
819  typename VEC::size_type offset)
820 {
821  return (data.size() >= (offset+DataTypes::noValues(shape)));
822 }
823 
827 template <class ResVEC, class LVEC, class RSCALAR>
828 void
829 binaryOpVectorRightScalar(ResVEC& res, // where result is to be stored
830  typename ResVEC::size_type resOffset, // offset in the result vector to start storing results
831  const typename ResVEC::size_type samplesToProcess, // number of samples to be updated in the result
832  const typename ResVEC::size_type sampleSize, // number of values in each sample
833  const LVEC& left, // LHS of calculation
834  typename LVEC::size_type leftOffset, // where to start reading LHS values
835  const RSCALAR* right, // RHS of the calculation
836  const bool rightreset, // true if RHS is providing a single sample of 1 value only
837  escript::ES_optype operation, // operation to perform
838  bool singleleftsample) // set to false for normal operation
839 {
840  size_t substep=(rightreset?0:1);
841  switch (operation)
842  {
843  case ADD:
844  {
845 #pragma omp parallel for
846  for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
847  {
848  typename LVEC::size_type leftbase=leftOffset+(singleleftsample?0:i*sampleSize);
849  const RSCALAR* rpos=right+(rightreset?0:i*substep);
850 
851  for (typename ResVEC::size_type j=0;j<sampleSize;++j)
852  {
853  res[i*sampleSize+resOffset+j]=left[leftbase+j]+*rpos;
854  }
855  }
856  }
857  break;
858  case POW:
859  {
860 #pragma omp parallel for
861  for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
862  {
863  typename LVEC::size_type leftbase=leftOffset+(singleleftsample?0:i*sampleSize);
864  const RSCALAR* rpos=right+(rightreset?0:i*substep);
865 
866  for (typename ResVEC::size_type j=0;j<sampleSize;++j)
867  {
868  res[i*sampleSize+resOffset+j]=pow(left[leftbase+j],*rpos);
869  }
870  }
871  }
872  break;
873  case SUB:
874  {
875 #pragma omp parallel for
876  for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
877  {
878  typename LVEC::size_type leftbase=leftOffset+(singleleftsample?0:i*sampleSize);
879  const RSCALAR* rpos=right+(rightreset?0:i*substep);
880 
881  for (typename ResVEC::size_type j=0;j<sampleSize;++j)
882  {
883  res[i*sampleSize+resOffset+j]=left[leftbase+j]-*rpos;
884  }
885  }
886  }
887  break;
888  case MUL:
889  {
890 #pragma omp parallel for
891  for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
892  {
893  typename LVEC::size_type leftbase=leftOffset+(singleleftsample?0:i*sampleSize);
894  const RSCALAR* rpos=right+(rightreset?0:i*substep);
895 
896  for (typename ResVEC::size_type j=0;j<sampleSize;++j)
897  {
898  res[i*sampleSize+resOffset+j]=left[leftbase+j] * *rpos;
899  }
900  }
901  }
902  break;
903  case DIV:
904  {
905 #pragma omp parallel for
906  for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
907  {
908  typename LVEC::size_type leftbase=leftOffset+(singleleftsample?0:i*sampleSize);
909  const RSCALAR* rpos=right+(rightreset?0:i*substep);
910 
911  for (typename ResVEC::size_type j=0;j<sampleSize;++j)
912  {
913  res[i*sampleSize+resOffset+j]=left[leftbase+j]/ *rpos;
914  }
915  }
916  }
917  break;
918  default:
919  throw DataException("Unsupported binary operation");
920  }
921 }
922 
923 template<>
924 void
925 binaryOpVectorRightScalar(DataTypes::RealVectorType& res, // where result is to be stored
926  typename DataTypes::RealVectorType::size_type resOffset, // offset in the result vector to start storing results
927  const typename DataTypes::RealVectorType::size_type samplesToProcess, // number of samples to be updated in the result
928  const typename DataTypes::RealVectorType::size_type sampleSize, // number of values in each sample
929  const DataTypes::RealVectorType& left, // LHS of calculation
930  typename DataTypes::RealVectorType::size_type leftOffset, // where to start reading LHS values
931  const DataTypes::real_t* right, // RHS of the calculation
932  const bool rightreset, // true if RHS is providing a single sample of 1 value only
933  escript::ES_optype operation, // operation to perform
934  bool singleleftsample);
935 
939 template <class ResVEC, class LSCALAR, class RVEC>
940 void
941 binaryOpVectorLeftScalar(ResVEC& res, // where result is to be stored
942  typename ResVEC::size_type resOffset, // offset in the result vector to start storing results
943  const typename ResVEC::size_type samplesToProcess, // number of samples to be updated in the result
944  const typename ResVEC::size_type sampleSize, // number of values in each sample
945  const LSCALAR* left, // LHS of calculation
946  const bool leftreset, // true if LHS is providing a single sample of 1 value only
947  const RVEC& right, // RHS of the calculation
948  typename RVEC::size_type rightOffset, // where to start reading RHS values
949  escript::ES_optype operation, // operation to perform
950  bool singlerightsample) // right consists of a single sample
951 {
952  size_t substep=(leftreset?0:1);
953  switch (operation)
954  {
955  case ADD:
956  {
957 #pragma omp parallel for
958  for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
959  {
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)
963  {
964  res[i*sampleSize+resOffset+j]=*lpos+right[rightbase+j];
965  }
966  }
967  }
968  break;
969  case POW:
970  {
971 #pragma omp parallel for
972  for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
973  {
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)
977  {
978  res[i*sampleSize+resOffset+j]=pow(*lpos,right[rightbase+j]);
979  }
980  }
981  }
982  break;
983  case SUB:
984  {
985 #pragma omp parallel for
986  for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
987  {
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)
991  {
992  res[i*sampleSize+resOffset+j]=*lpos-right[rightbase+j];
993  }
994  }
995  }
996  break;
997  case MUL:
998  {
999 #pragma omp parallel for
1000  for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
1001  {
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)
1005  {
1006  res[i*sampleSize+resOffset+j]=*lpos*right[rightbase+j];
1007  }
1008  }
1009  }
1010  break;
1011  case DIV:
1012  {
1013 #pragma omp parallel for
1014  for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
1015  {
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)
1019  {
1020  res[i*sampleSize+resOffset+j]=*lpos/right[rightbase+j];
1021  }
1022  }
1023  }
1024  break;
1025  default:
1026  throw DataException("Unsupported binary operation");
1027  }
1028 }
1029 
1030 template <>
1031 void
1032 binaryOpVectorLeftScalar(DataTypes::RealVectorType& res, // where result is to be stored
1033  typename DataTypes::RealVectorType::size_type resOffset, // offset in the result vector to start storing results
1034  const typename DataTypes::RealVectorType::size_type samplesToProcess, // number of samples to be updated in the result
1035  const typename DataTypes::RealVectorType::size_type sampleSize, // number of values in each sample
1036  const DataTypes::real_t* left, // LHS of calculation
1037  const bool leftreset, // true if LHS is providing a single sample of 1 value only
1038  const DataTypes::RealVectorType& right, // RHS of the calculation
1039  typename DataTypes::RealVectorType::size_type rightOffset, // where to start reading RHS values
1040  escript::ES_optype operation, // operation to perform
1041  bool singlerightsample); // right consists of a single sample
1042 
1046 template <class ResVEC, class LVEC, class RVEC>
1047 void
1048 binaryOpVector(ResVEC& res, // where result is to be stored
1049  typename ResVEC::size_type resOffset, // offset in the result vector to start storing results
1050  const typename ResVEC::size_type samplesToProcess, // number of samples to be updated in the result
1051  const typename ResVEC::size_type sampleSize, // number of values in each sample
1052  const LVEC& left, // LHS of calculation
1053  typename LVEC::size_type leftOffset, // where to start reading LHS values
1054  const bool leftreset, // Is LHS only supplying a single sample instead of a bunch of them
1055  const RVEC& right, // RHS of the calculation
1056  typename RVEC::size_type rightOffset, // where to start reading RHS values
1057  const bool rightreset, // Is RHS only supplying a single sample instead of a bunch of them
1058  escript::ES_optype operation) // operation to perform
1059 {
1060  switch (operation)
1061  {
1062  case ADD:
1063  {
1064 #pragma omp parallel for
1065  for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
1066  {
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)
1070  {
1071  res[i*sampleSize+resOffset+j]=left[leftbase+j]+right[rightbase+j];
1072  }
1073  }
1074  }
1075  break;
1076  case POW:
1077  {
1078 #pragma omp parallel for
1079  for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
1080  {
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)
1084  {
1085  res[i*sampleSize+resOffset+j]=pow(left[leftbase+j],right[rightbase+j]);
1086  }
1087  }
1088  }
1089  break;
1090  case SUB:
1091  {
1092 #pragma omp parallel for
1093  for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
1094  {
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)
1098  {
1099  res[i*sampleSize+resOffset+j]=left[leftbase+j]-right[rightbase+j];
1100  }
1101  }
1102  }
1103  break;
1104  case MUL:
1105  {
1106 #pragma omp parallel for
1107  for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
1108  {
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)
1112  {
1113  res[i*sampleSize+resOffset+j]=left[leftbase+j]*right[rightbase+j];
1114  }
1115  }
1116  }
1117  break;
1118  case DIV:
1119  {
1120 #pragma omp parallel for
1121  for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
1122  {
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)
1126  {
1127  res[i*sampleSize+resOffset+j]=left[leftbase+j]/right[rightbase+j];
1128  }
1129  }
1130  }
1131  break;
1132  default:
1133  throw DataException("Unsupported binary operation");
1134  }
1135 }
1136 
1137 template <>
1138 void
1139 binaryOpVector(DataTypes::RealVectorType& res, // where result is to be stored
1140  typename DataTypes::RealVectorType::size_type resOffset, // offset in the result vector to start storing results
1141  const typename DataTypes::RealVectorType::size_type samplesToProcess, // number of samples to be updated in the result
1142  const typename DataTypes::RealVectorType::size_type sampleSize, // number of values in each sample
1143  const DataTypes::RealVectorType& left, // LHS of calculation
1144  typename DataTypes::RealVectorType::size_type leftOffset, // where to start reading LHS values
1145  const bool leftreset, // Is LHS only supplying a single sample instead of a bunch of them
1146  const DataTypes::RealVectorType& right, // RHS of the calculation
1147  typename DataTypes::RealVectorType::size_type rightOffset, // where to start reading RHS values
1148  const bool rightreset, // Is RHS only supplying a single sample instead of a bunch of them
1149  escript::ES_optype operation); // operation to perform
1150 
1151 #define OPVECLAZYBODY(X) \
1152  for (size_t j=0;j<onumsteps;++j)\
1153  {\
1154  for (size_t i=0;i<numsteps;++i,res+=resultStep) \
1155  { \
1156  for (size_t s=0; s<chunksize; ++s)\
1157  {\
1158  res[s] = X;\
1159  }\
1160 /* tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X);*/ \
1161  lroffset+=leftstep; \
1162  rroffset+=rightstep; \
1163  }\
1164  lroffset+=oleftstep;\
1165  rroffset+=orightstep;\
1166  }
1167 
1174 template <class ResELT, class LELT, class RELT>
1175 void
1177  const LELT* left,
1178  const RELT* right,
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,
1187  size_t lroffset,
1188  size_t rroffset,
1189  escript::ES_optype operation) // operation to perform
1190 {
1191  switch (operation)
1192  {
1193  case ADD:
1194  OPVECLAZYBODY((left[lroffset+s]+right[rroffset+s]));
1195  break;
1196  case POW:
1197  OPVECLAZYBODY(pow(left[lroffset+s],right[rroffset+s]))
1198  break;
1199  case SUB:
1200  OPVECLAZYBODY(left[lroffset+s]-right[rroffset+s])
1201  break;
1202  case MUL:
1203  OPVECLAZYBODY(left[lroffset+s]*right[rroffset+s])
1204  break;
1205  case DIV:
1206  OPVECLAZYBODY(left[lroffset+s]/right[rroffset+s])
1207  break;
1208  default:
1209  ESYS_ASSERT(false, "Invalid operation. This should never happen!");
1210  // I can't throw here because this will be called inside a parallel section
1211  }
1212 }
1213 
1214 
1221 template <class ResELT, class LELT, class RELT>
1222 void
1224  const LELT* left,
1225  const RELT* right,
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,
1234  size_t lroffset,
1235  size_t rroffset,
1236  escript::ES_optype operation) // operation to perform
1237 {
1238  switch (operation)
1239  {
1240  case LESS:
1241  OPVECLAZYBODY(left[lroffset+s]<right[rroffset+s])
1242  break;
1243  case GREATER:
1244  OPVECLAZYBODY(left[lroffset+s]>right[rroffset+s])
1245  break;
1246  case GREATER_EQUAL:
1247  OPVECLAZYBODY(left[lroffset+s]>=right[rroffset+s])
1248  break;
1249  case LESS_EQUAL:
1250  OPVECLAZYBODY(left[lroffset+s]<=right[rroffset+s])
1251  break;
1252  default:
1253  ESYS_ASSERT(false, "Invalid operation. This should never happen!");
1254  // I can't throw here because this will be called inside a parallel section
1255  }
1256 }
1257 
1261 /* trying to make a single version for all Tagged+Expanded interactions */
1262 template <class ResVEC, class LVEC, class RVEC>
1263 void
1264 binaryOpVectorTagged(ResVEC& res, // where result is to be stored
1265  const typename ResVEC::size_type samplesToProcess, // number of samples to be updated in the result
1266  const typename ResVEC::size_type DPPSample, // number of datapoints per sample
1267  const typename ResVEC::size_type DPSize, // datapoint size
1268  const LVEC& left, // LHS of calculation
1269  bool leftscalar,
1270  const RVEC& right, // RHS of the calculation
1271  bool rightscalar,
1272  bool lefttagged, // true if left object is the tagged one
1273  const DataTagged& tagsource, // where to get tag offsets from
1274  escript::ES_optype operation) // operation to perform
1275 {
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;
1279  switch (operation)
1280  {
1281  case ADD:
1282  {
1283 #pragma omp parallel for
1284  for (typename ResVEC::size_type i=0;i<limit;++i)
1285  {
1286  typename LVEC::size_type leftbase=(lefttagged?tagsource.getPointOffset(i/DPPSample,0):i*lstep); // only one of these
1287  typename RVEC::size_type rightbase=(lefttagged?i*rstep:tagsource.getPointOffset(i/DPPSample,0)); // will apply
1288 
1289  for (typename ResVEC::size_type j=0;j<DPSize;++j)
1290  {
1291  res[i*DPSize+j]=left[leftbase+j*(!leftscalar)]+right[rightbase+j*(!rightscalar)];
1292  }
1293  }
1294  }
1295  break;
1296  case POW:
1297  {
1298 #pragma omp parallel for
1299  for (typename ResVEC::size_type i=0;i<limit;++i)
1300  {
1301  typename LVEC::size_type leftbase=(lefttagged?tagsource.getPointOffset(i/DPPSample,0):i*lstep); // only one of these
1302  typename RVEC::size_type rightbase=(lefttagged?i*rstep:tagsource.getPointOffset(i/DPPSample,0)); // will apply
1303 
1304  for (typename ResVEC::size_type j=0;j<DPSize;++j)
1305  {
1306  res[i*DPSize+j]=pow(left[leftbase+j*(!leftscalar)],right[rightbase+j*(!rightscalar)]);
1307  }
1308  }
1309  }
1310  break;
1311  case SUB:
1312  {
1313 #pragma omp parallel for
1314  for (typename ResVEC::size_type i=0;i<limit;++i)
1315  {
1316  typename LVEC::size_type leftbase=(lefttagged?tagsource.getPointOffset(i/DPPSample,0):i*lstep); // only one of these
1317  typename RVEC::size_type rightbase=(lefttagged?i*rstep:tagsource.getPointOffset(i/DPPSample,0)); // will apply
1318 
1319  for (typename ResVEC::size_type j=0;j<DPSize;++j)
1320  {
1321  res[i*DPSize+j]=left[leftbase+j*(!leftscalar)]-right[rightbase+j*(!rightscalar)];
1322  }
1323  }
1324  }
1325  break;
1326  case MUL:
1327  {
1328 #pragma omp parallel for
1329  for (typename ResVEC::size_type i=0;i<limit;++i)
1330  {
1331  typename LVEC::size_type leftbase=(lefttagged?tagsource.getPointOffset(i/DPPSample,0):i*lstep); // only one of these
1332  typename RVEC::size_type rightbase=(lefttagged?i*rstep:tagsource.getPointOffset(i/DPPSample,0)); // will apply
1333 
1334  for (typename ResVEC::size_type j=0;j<DPSize;++j)
1335  {
1336  res[i*DPSize+j]=left[leftbase+j*(!leftscalar)]*right[rightbase+j*(!rightscalar)];
1337  }
1338  }
1339  }
1340  break;
1341  case DIV:
1342  {
1343 #pragma omp parallel for
1344  for (typename ResVEC::size_type i=0;i<limit;++i)
1345  {
1346  typename LVEC::size_type leftbase=(lefttagged?tagsource.getPointOffset(i/DPPSample,0):i*lstep); // only one of these
1347  typename RVEC::size_type rightbase=(lefttagged?i*rstep:tagsource.getPointOffset(i/DPPSample,0)); // will apply
1348 
1349  for (typename ResVEC::size_type j=0;j<DPSize;++j)
1350  {
1351  res[i*DPSize+j]=left[leftbase+j*(!leftscalar)]/right[rightbase+j*(!rightscalar)];
1352  }
1353  }
1354  }
1355  break;
1356  default:
1357  throw DataException("Unsupported binary operation");
1358  }
1359 }
1360 
1361 template<>
1362 void
1363 binaryOpVectorTagged(DataTypes::RealVectorType& res, // where result is to be stored
1364  const typename DataTypes::RealVectorType::size_type samplesToProcess, // number of samples to be updated in the result
1365  const typename DataTypes::RealVectorType::size_type DPPSample, // number of datapoints per sample
1366  const typename DataTypes::RealVectorType::size_type DPSize, // datapoint size
1367  const DataTypes::RealVectorType& left, // LHS of calculation
1368  const bool leftscalar,
1369  const DataTypes::RealVectorType& right, // RHS of the calculation
1370  const bool rightscalar,
1371  const bool lefttagged, // true if left object is the tagged one
1372  const DataTagged& tagsource, // where to get tag offsets from
1373  escript::ES_optype operation);
1374 
1375 
1376 
1377 
1394 template <class BinaryFunction>
1395 inline
1398  const DataTypes::ShapeType& leftShape,
1400  BinaryFunction operation,
1401  DataTypes::real_t initial_value)
1402 {
1403  ESYS_ASSERT((left.size()>0)&&checkOffset(left,leftShape,offset),
1404  "Couldn't perform reductionOp due to insufficient storage.");
1405  DataTypes::real_t current_value=initial_value;
1406  for (DataTypes::RealVectorType::size_type i=0;i<DataTypes::noValues(leftShape);i++) {
1407  current_value=operation(current_value,left[offset+i]);
1408  }
1409  return current_value;
1410 }
1411 
1412 template <class BinaryFunction>
1413 inline
1416  const DataTypes::ShapeType& leftShape,
1418  BinaryFunction operation,
1419  DataTypes::real_t initial_value)
1420 {
1421  ESYS_ASSERT((left.size()>0)&&checkOffset(left,leftShape,offset),
1422  "Couldn't perform reductionOp due to insufficient storage.");
1423  DataTypes::real_t current_value=initial_value;
1424  for (DataTypes::RealVectorType::size_type i=0;i<DataTypes::noValues(leftShape);i++) {
1425  current_value=operation(current_value,left[offset+i]);
1426  }
1427  return current_value;
1428 }
1429 
1430 
1447 int
1449  const DataTypes::ShapeType& inShape,
1452  const DataTypes::ShapeType& outShape,
1454  int count,
1455  LapackInverseHelper& helper);
1456 
1464 void
1465 matrixInverseError(int err);
1466 
1471 inline
1472 bool
1474 {
1475  for (size_t z=inOffset;z<inOffset+count;++z)
1476  {
1477  if (nancheck(in[z]))
1478  {
1479  return true;
1480  }
1481  }
1482  return false;
1483 }
1484 
1485 inline
1486 bool
1488 {
1489  for (size_t z=inOffset;z<inOffset+count;++z)
1490  {
1491  if (nancheck(in[z]))
1492  {
1493  return true;
1494  }
1495  }
1496  return false;
1497 }
1498 
1499 
1500 } // end namespace escript
1501 
1502 #endif // __ESCRIPT_DATAMATHS_H__
1503 
#define ESYS_ASSERT(a, b)
EsysAssert is a MACRO that will throw an exception if the boolean condition specified is false.
Definition: Assert.h:79
#define OPVECLAZYBODY(X)
Definition: DataVectorOps.h:1151
#define V(_K_, _I_)
Definition: ShapeFunctions.cpp:121
Definition: DataException.h:28
Simulates a full dataset accessible via sampleNo and dataPointNo.
Definition: DataTagged.h:46
virtual DataTypes::RealVectorType::size_type getPointOffset(int sampleNo, int dataPointNo) const
getPointOffset
Definition: DataTagged.cpp:982
size_type size() const
Return the number of elements in this DataVectorAlt.
Definition: DataVectorAlt.h:215
real_t ElementType
Definition: DataVectorAlt.h:43
DataTypes::vec_size_type size_type
Definition: DataVectorAlt.h:50
#define ESCRIPT_DLL_API
Definition: escriptcore/src/system_dep.h:30
double real_t
type of all real-valued scalars in escript
Definition: DataTypes.h:52
escript::DataTypes::DataVectorAlt< real_t > RealVectorType
Vector to store underlying data.
Definition: DataVector.h:44
int noValues(const ShapeType &shape)
Calculate the number of values in a datapoint with the given shape.
Definition: DataTypes.cpp:91
int getRank(const DataTypes::ShapeType &shape)
Return the rank (number of dimensions) of the given shape.
Definition: DataTypes.h:225
std::vector< int > ShapeType
The shape of a single datapoint.
Definition: DataTypes.h:44
escript::DataTypes::DataVectorAlt< cplx_t > CplxVectorType
Definition: DataVector.h:45
vec_size_type getRelIndex(const DataTypes::ShapeType &shape, vec_size_type i)
Compute the offset (in 1D vector) of a given subscript with a shape.
Definition: DataTypes.h:240
Definition: AbstractContinuousDomain.cpp:23
DataTypes::real_t reductionOpVector(const DataTypes::RealVectorType &left, const DataTypes::ShapeType &leftShape, DataTypes::RealVectorType::size_type offset, BinaryFunction operation, DataTypes::real_t initial_value)
Perform the given data point reduction operation on the data point specified by the given offset into...
Definition: DataVectorOps.h:1397
void eigenvalues1(const DataTypes::real_t A00, DataTypes::real_t *ev0)
solves a 1x1 eigenvalue A*V=ev*V problem
Definition: ArrayOps.h:164
void eigenvalues_and_eigenvectors(const DataTypes::RealVectorType &in, const DataTypes::ShapeType &inShape, DataTypes::RealVectorType::size_type inOffset, DataTypes::RealVectorType &ev, const DataTypes::ShapeType &evShape, DataTypes::RealVectorType::size_type evOffset, DataTypes::RealVectorType &V, const DataTypes::ShapeType &VShape, DataTypes::RealVectorType::size_type VOffset, const double tol=1.e-13)
solves a local eigenvalue problem
Definition: DataVectorOps.h:750
void matMult(const DataTypes::RealVectorType &left, const DataTypes::ShapeType &leftShape, DataTypes::RealVectorType::size_type leftOffset, const DataTypes::RealVectorType &right, const DataTypes::ShapeType &rightShape, DataTypes::RealVectorType::size_type rightOffset, DataTypes::RealVectorType &result, const DataTypes::ShapeType &resultShape)
Perform a matrix multiply of the given views.
Definition: DataVectorOps.cpp:39
void eigenvalues(const DataTypes::RealVectorType &in, const DataTypes::ShapeType &inShape, typename DataTypes::RealVectorType::size_type inOffset, DataTypes::RealVectorType &ev, const DataTypes::ShapeType &evShape, typename DataTypes::RealVectorType::size_type evOffset)
solves a local eigenvalue problem
Definition: DataVectorOps.h:639
void eigenvalues3(const DataTypes::real_t A00, const DataTypes::real_t A01, const DataTypes::real_t A02, const DataTypes::real_t A11, const DataTypes::real_t A12, const DataTypes::real_t A22, DataTypes::real_t *ev0, DataTypes::real_t *ev1, DataTypes::real_t *ev2)
solves a 3x3 eigenvalue A*V=ev*V problem for symmetric A
Definition: ArrayOps.h:210
void binaryOpVectorLeftScalar(DataTypes::RealVectorType &res, typename DataTypes::RealVectorType::size_type resOffset, const typename DataTypes::RealVectorType::size_type samplesToProcess, const typename DataTypes::RealVectorType::size_type sampleSize, const DataTypes::real_t *left, const bool leftreset, const DataTypes::RealVectorType &right, typename DataTypes::RealVectorType::size_type rightOffset, escript::ES_optype operation, bool singlerightsample)
Definition: DataVectorOps.cpp:639
void eigenvalues_and_eigenvectors3(const DataTypes::real_t A00, const DataTypes::real_t A01, const DataTypes::real_t A02, const DataTypes::real_t A11, const DataTypes::real_t A12, const DataTypes::real_t A22, DataTypes::real_t *ev0, DataTypes::real_t *ev1, DataTypes::real_t *ev2, DataTypes::real_t *V00, DataTypes::real_t *V10, DataTypes::real_t *V20, DataTypes::real_t *V01, DataTypes::real_t *V11, DataTypes::real_t *V21, DataTypes::real_t *V02, DataTypes::real_t *V12, DataTypes::real_t *V22, const DataTypes::real_t tol)
solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A. Eigenvectors are ordered by increasing valu...
Definition: ArrayOps.h:454
void symmetric(const VEC &in, const DataTypes::ShapeType &inShape, typename VEC::size_type inOffset, VEC &ev, const DataTypes::ShapeType &evShape, typename VEC::size_type evOffset)
computes a symmetric matrix from your square matrix A: (A + transpose(A)) / 2
Definition: DataVectorOps.h:103
void binaryOpVector(DataTypes::RealVectorType &res, typename DataTypes::RealVectorType::size_type resOffset, const typename DataTypes::RealVectorType::size_type samplesToProcess, const typename DataTypes::RealVectorType::size_type sampleSize, const DataTypes::RealVectorType &left, typename DataTypes::RealVectorType::size_type leftOffset, const bool leftreset, const DataTypes::RealVectorType &right, typename DataTypes::RealVectorType::size_type rightOffset, const bool rightreset, escript::ES_optype operation)
Definition: DataVectorOps.cpp:768
void swapaxes(const VEC &in, const DataTypes::ShapeType &inShape, typename VEC::size_type inOffset, VEC &ev, const DataTypes::ShapeType &evShape, typename VEC::size_type evOffset, int axis0, int axis1)
swaps the components axis0 and axis1.
Definition: DataVectorOps.h:488
void binaryOpVectorTagged(DataTypes::RealVectorType &res, const typename DataTypes::RealVectorType::size_type samplesToProcess, const typename DataTypes::RealVectorType::size_type DPPSample, const typename DataTypes::RealVectorType::size_type DPSize, const DataTypes::RealVectorType &left, const bool leftscalar, const DataTypes::RealVectorType &right, const bool rightscalar, const bool lefttagged, const DataTagged &tagsource, escript::ES_optype operation)
Definition: DataVectorOps.cpp:340
void eigenvalues_and_eigenvectors2(const DataTypes::real_t A00, const DataTypes::real_t A01, const DataTypes::real_t A11, DataTypes::real_t *ev0, DataTypes::real_t *ev1, DataTypes::real_t *V00, DataTypes::real_t *V10, DataTypes::real_t *V01, DataTypes::real_t *V11, const DataTypes::real_t tol)
solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A. Eigenvectors are ordered by increasing valu...
Definition: ArrayOps.h:345
void transpose(const VEC &in, const DataTypes::ShapeType &inShape, typename VEC::size_type inOffset, VEC &ev, const DataTypes::ShapeType &evShape, typename VEC::size_type evOffset, int axis_offset)
Transpose each data point of this Data object around the given axis.
Definition: DataVectorOps.h:343
bool vectorHasNaN(const DataTypes::RealVectorType &in, DataTypes::RealVectorType::size_type inOffset, size_t count)
returns true if the vector contains NaN
Definition: DataVectorOps.h:1473
bool checkOffset(const VEC &data, const DataTypes::ShapeType &shape, typename VEC::size_type offset)
Definition: DataVectorOps.h:817
void antisymmetric(const VEC &in, const DataTypes::ShapeType &inShape, typename VEC::size_type inOffset, VEC &ev, const DataTypes::ShapeType &evShape, typename VEC::size_type evOffset)
computes a antisymmetric matrix from your square matrix A: (A - transpose(A)) / 2
Definition: DataVectorOps.h:152
void eigenvalues2(const T A00, const T A01, const T A11, T *ev0, T *ev1)
solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A
Definition: ArrayOps.h:186
void eigenvalues_and_eigenvectors1(const DataTypes::real_t A00, DataTypes::real_t *ev0, DataTypes::real_t *V00, const DataTypes::real_t tol)
solves a 1x1 eigenvalue A*V=ev*V problem for symmetric A
Definition: ArrayOps.h:253
int matrix_inverse(const DataTypes::RealVectorType &in, const DataTypes::ShapeType &inShape, DataTypes::RealVectorType::size_type inOffset, DataTypes::RealVectorType &out, const DataTypes::ShapeType &outShape, DataTypes::RealVectorType::size_type outOffset, int count, LapackInverseHelper &helper)
computes the inverses of square (up to 3x3) matricies
Definition: DataVectorOps.cpp:207
DataTypes::ShapeType determineResultShape(const DataTypes::ShapeType &left, const DataTypes::ShapeType &right)
Determine the shape of the result array for a matrix multiplication of the given views.
Definition: DataVectorOps.cpp:170
void antihermitian(const DataTypes::CplxVectorType &in, const DataTypes::ShapeType &inShape, typename DataTypes::CplxVectorType::size_type inOffset, DataTypes::CplxVectorType &ev, const DataTypes::ShapeType &evShape, typename DataTypes::CplxVectorType::size_type evOffset)
computes a antihermitian matrix from your square matrix A: (A - adjoint(A)) / 2
Definition: DataVectorOps.cpp:963
bool nancheck(DataTypes::real_t d)
acts as a wrapper to isnan.
Definition: ArrayOps.h:120
void matrixInverseError(int err)
throws an appropriate exception based on failure of matrix_inverse.
Definition: DataVectorOps.cpp:186
ES_optype
Definition: ES_optype.h:30
@ GREATER
Definition: ES_optype.h:81
@ POW
Definition: ES_optype.h:37
@ LESS
Definition: ES_optype.h:80
@ MUL
Definition: ES_optype.h:35
@ LESS_EQUAL
Definition: ES_optype.h:83
@ SUB
Definition: ES_optype.h:34
@ DIV
Definition: ES_optype.h:36
@ ADD
Definition: ES_optype.h:33
@ GREATER_EQUAL
Definition: ES_optype.h:82
void binaryOpVectorLazyArithmeticHelper(ResELT *res, const LELT *left, const RELT *right, const size_t chunksize, const size_t onumsteps, const size_t numsteps, const size_t resultStep, const size_t leftstep, const size_t rightstep, const size_t oleftstep, const size_t orightstep, size_t lroffset, size_t rroffset, escript::ES_optype operation)
Definition: DataVectorOps.h:1176
void trace(const VEC &in, const DataTypes::ShapeType &inShape, typename VEC::size_type inOffset, VEC &ev, const DataTypes::ShapeType &evShape, typename VEC::size_type evOffset, int axis_offset)
computes the trace of a matrix
Definition: DataVectorOps.h:242
void binaryOpVectorRightScalar(DataTypes::RealVectorType &res, typename DataTypes::RealVectorType::size_type resOffset, const typename DataTypes::RealVectorType::size_type samplesToProcess, const typename DataTypes::RealVectorType::size_type sampleSize, const DataTypes::RealVectorType &left, typename DataTypes::RealVectorType::size_type leftOffset, const DataTypes::real_t *right, const bool rightreset, escript::ES_optype operation, bool singleleftsample)
Definition: DataVectorOps.cpp:500
void hermitian(const DataTypes::CplxVectorType &in, const DataTypes::ShapeType &inShape, DataTypes::CplxVectorType::size_type inOffset, DataTypes::CplxVectorType &ev, const DataTypes::ShapeType &evShape, DataTypes::CplxVectorType::size_type evOffset)
computes an hermitian matrix from your square matrix A: (A + adjoint(A)) / 2
Definition: DataVectorOps.cpp:916
void binaryOpVectorLazyRelationalHelper(ResELT *res, const LELT *left, const RELT *right, const size_t chunksize, const size_t onumsteps, const size_t numsteps, const size_t resultStep, const size_t leftstep, const size_t rightstep, const size_t oleftstep, const size_t orightstep, size_t lroffset, size_t rroffset, escript::ES_optype operation)
Definition: DataVectorOps.h:1223