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