escript  Revision_
DataVectorOps.h
Go to the documentation of this file.
1 
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2016 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  typename DataTypes::CplxVectorType::ElementType in00,in10,in20,in01,in11,in21,in02,in12,in22;
691  typename DataTypes::CplxVectorType::ElementType ev0,ev1,ev2;
692  int s=inShape[0];
693  if (s==1) {
694  in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
695  eigenvalues1(in00,&ev0);
696  ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
697 
698  } else if (s==2) {
699  in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
700  in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
701  in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
702  in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
703  eigenvalues2(in00,(in01+in10)/2.,in11,&ev0,&ev1);
704  ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
705  ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
706 
707  } else if (s==3) {
708  // this doesn't work yet
709 // in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
710 // in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
711 // in20=in[inOffset+DataTypes::getRelIndex(inShape,2,0)];
712 // in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
713 // in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
714 // in21=in[inOffset+DataTypes::getRelIndex(inShape,2,1)];
715 // in02=in[inOffset+DataTypes::getRelIndex(inShape,0,2)];
716 // in12=in[inOffset+DataTypes::getRelIndex(inShape,1,2)];
717 // in22=in[inOffset+DataTypes::getRelIndex(inShape,2,2)];
718 // eigenvalues3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
719 // &ev0,&ev1,&ev2);
720 // ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
721 // ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
722 // ev[evOffset+DataTypes::getRelIndex(evShape,2)]=ev2;
723 
724  }
725  }
726 
727 
744  inline
745  void
752  const double tol=1.e-13)
753  {
754  double in00,in10,in20,in01,in11,in21,in02,in12,in22;
755  double V00,V10,V20,V01,V11,V21,V02,V12,V22;
756  double ev0,ev1,ev2;
757  int s=inShape[0];
758  if (s==1) {
759  in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
760  eigenvalues_and_eigenvectors1(in00,&ev0,&V00,tol);
761  ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
762  V[inOffset+DataTypes::getRelIndex(VShape,0,0)]=V00;
763  } else if (s==2) {
764  in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
765  in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
766  in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
767  in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
768  eigenvalues_and_eigenvectors2(in00,(in01+in10)/2.,in11,
769  &ev0,&ev1,&V00,&V10,&V01,&V11,tol);
770  ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
771  ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
772  V[inOffset+DataTypes::getRelIndex(VShape,0,0)]=V00;
773  V[inOffset+DataTypes::getRelIndex(VShape,1,0)]=V10;
774  V[inOffset+DataTypes::getRelIndex(VShape,0,1)]=V01;
775  V[inOffset+DataTypes::getRelIndex(VShape,1,1)]=V11;
776  } else if (s==3) {
777  in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
778  in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
779  in20=in[inOffset+DataTypes::getRelIndex(inShape,2,0)];
780  in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
781  in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
782  in21=in[inOffset+DataTypes::getRelIndex(inShape,2,1)];
783  in02=in[inOffset+DataTypes::getRelIndex(inShape,0,2)];
784  in12=in[inOffset+DataTypes::getRelIndex(inShape,1,2)];
785  in22=in[inOffset+DataTypes::getRelIndex(inShape,2,2)];
786  eigenvalues_and_eigenvectors3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
787  &ev0,&ev1,&ev2,
788  &V00,&V10,&V20,&V01,&V11,&V21,&V02,&V12,&V22,tol);
789  ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
790  ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
791  ev[evOffset+DataTypes::getRelIndex(evShape,2)]=ev2;
792  V[inOffset+DataTypes::getRelIndex(VShape,0,0)]=V00;
793  V[inOffset+DataTypes::getRelIndex(VShape,1,0)]=V10;
794  V[inOffset+DataTypes::getRelIndex(VShape,2,0)]=V20;
795  V[inOffset+DataTypes::getRelIndex(VShape,0,1)]=V01;
796  V[inOffset+DataTypes::getRelIndex(VShape,1,1)]=V11;
797  V[inOffset+DataTypes::getRelIndex(VShape,2,1)]=V21;
798  V[inOffset+DataTypes::getRelIndex(VShape,0,2)]=V02;
799  V[inOffset+DataTypes::getRelIndex(VShape,1,2)]=V12;
800  V[inOffset+DataTypes::getRelIndex(VShape,2,2)]=V22;
801 
802  }
803  }
804 
805 
810 template <class VEC>
811 inline
812 bool
813 checkOffset(const VEC& data,
814  const DataTypes::ShapeType& shape,
815  typename VEC::size_type offset)
816 {
817  return (data.size() >= (offset+DataTypes::noValues(shape)));
818 }
819 
823 template <class ResVEC, class LVEC, class RSCALAR>
824 void
825 binaryOpVectorRightScalar(ResVEC& res, // where result is to be stored
826  typename ResVEC::size_type resOffset, // offset in the result vector to start storing results
827  const typename ResVEC::size_type samplesToProcess, // number of samples to be updated in the result
828  const typename ResVEC::size_type sampleSize, // number of values in each sample
829  const LVEC& left, // LHS of calculation
830  typename LVEC::size_type leftOffset, // where to start reading LHS values
831  const RSCALAR* right, // RHS of the calculation
832  const bool rightreset, // true if RHS is providing a single sample of 1 value only
833  escript::ES_optype operation, // operation to perform
834  bool singleleftsample) // set to false for normal operation
835 {
836  size_t substep=(rightreset?0:1);
837  switch (operation)
838  {
839  case ADD:
840  {
841 #pragma omp parallel for
842  for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
843  {
844  typename LVEC::size_type leftbase=leftOffset+(singleleftsample?0:i*sampleSize);
845  const RSCALAR* rpos=right+(rightreset?0:i*substep);
846 
847  for (typename ResVEC::size_type j=0;j<sampleSize;++j)
848  {
849  res[i*sampleSize+resOffset+j]=left[leftbase+j]+*rpos;
850  }
851  }
852  }
853  break;
854  case POW:
855  {
856 #pragma omp parallel for
857  for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
858  {
859  typename LVEC::size_type leftbase=leftOffset+(singleleftsample?0:i*sampleSize);
860  const RSCALAR* rpos=right+(rightreset?0:i*substep);
861 
862  for (typename ResVEC::size_type j=0;j<sampleSize;++j)
863  {
864  res[i*sampleSize+resOffset+j]=pow(left[leftbase+j],*rpos);
865  }
866  }
867  }
868  break;
869  case SUB:
870  {
871 #pragma omp parallel for
872  for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
873  {
874  typename LVEC::size_type leftbase=leftOffset+(singleleftsample?0:i*sampleSize);
875  const RSCALAR* rpos=right+(rightreset?0:i*substep);
876 
877  for (typename ResVEC::size_type j=0;j<sampleSize;++j)
878  {
879  res[i*sampleSize+resOffset+j]=left[leftbase+j]-*rpos;
880  }
881  }
882  }
883  break;
884  case MUL:
885  {
886 #pragma omp parallel for
887  for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
888  {
889  typename LVEC::size_type leftbase=leftOffset+(singleleftsample?0:i*sampleSize);
890  const RSCALAR* rpos=right+(rightreset?0:i*substep);
891 
892  for (typename ResVEC::size_type j=0;j<sampleSize;++j)
893  {
894  res[i*sampleSize+resOffset+j]=left[leftbase+j] * *rpos;
895  }
896  }
897  }
898  break;
899  case DIV:
900  {
901 #pragma omp parallel for
902  for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
903  {
904  typename LVEC::size_type leftbase=leftOffset+(singleleftsample?0:i*sampleSize);
905  const RSCALAR* rpos=right+(rightreset?0:i*substep);
906 
907  for (typename ResVEC::size_type j=0;j<sampleSize;++j)
908  {
909  res[i*sampleSize+resOffset+j]=left[leftbase+j]/ *rpos;
910  }
911  }
912  }
913  break;
914  default:
915  throw DataException("Unsupported binary operation");
916  }
917 }
918 
919 template<>
920 void
921 binaryOpVectorRightScalar(DataTypes::RealVectorType& res, // where result is to be stored
922  typename DataTypes::RealVectorType::size_type resOffset, // offset in the result vector to start storing results
923  const typename DataTypes::RealVectorType::size_type samplesToProcess, // number of samples to be updated in the result
924  const typename DataTypes::RealVectorType::size_type sampleSize, // number of values in each sample
925  const DataTypes::RealVectorType& left, // LHS of calculation
926  typename DataTypes::RealVectorType::size_type leftOffset, // where to start reading LHS values
927  const DataTypes::real_t* right, // RHS of the calculation
928  const bool rightreset, // true if RHS is providing a single sample of 1 value only
929  escript::ES_optype operation, // operation to perform
930  bool singleleftsample);
931 
935 template <class ResVEC, class LSCALAR, class RVEC>
936 void
937 binaryOpVectorLeftScalar(ResVEC& res, // where result is to be stored
938  typename ResVEC::size_type resOffset, // offset in the result vector to start storing results
939  const typename ResVEC::size_type samplesToProcess, // number of samples to be updated in the result
940  const typename ResVEC::size_type sampleSize, // number of values in each sample
941  const LSCALAR* left, // LHS of calculation
942  const bool leftreset, // true if LHS is providing a single sample of 1 value only
943  const RVEC& right, // RHS of the calculation
944  typename RVEC::size_type rightOffset, // where to start reading RHS values
945  escript::ES_optype operation, // operation to perform
946  bool singlerightsample) // right consists of a single sample
947 {
948  size_t substep=(leftreset?0:1);
949  switch (operation)
950  {
951  case ADD:
952  {
953 #pragma omp parallel for
954  for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
955  {
956  typename RVEC::size_type rightbase=rightOffset+(singlerightsample?0:i*sampleSize);
957  const LSCALAR* lpos=left+(leftreset?0:i*substep);
958  for (typename ResVEC::size_type j=0;j<sampleSize;++j)
959  {
960  res[i*sampleSize+resOffset+j]=*lpos+right[rightbase+j];
961  }
962  }
963  }
964  break;
965  case POW:
966  {
967 #pragma omp parallel for
968  for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
969  {
970  typename RVEC::size_type rightbase=rightOffset+(singlerightsample?0:i*sampleSize);
971  const LSCALAR* lpos=left+(leftreset?0:i*substep);
972  for (typename ResVEC::size_type j=0;j<sampleSize;++j)
973  {
974  res[i*sampleSize+resOffset+j]=pow(*lpos,right[rightbase+j]);
975  }
976  }
977  }
978  break;
979  case SUB:
980  {
981 #pragma omp parallel for
982  for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
983  {
984  typename RVEC::size_type rightbase=rightOffset+(singlerightsample?0:i*sampleSize);
985  const LSCALAR* lpos=left+(leftreset?0:i*substep);
986  for (typename ResVEC::size_type j=0;j<sampleSize;++j)
987  {
988  res[i*sampleSize+resOffset+j]=*lpos-right[rightbase+j];
989  }
990  }
991  }
992  break;
993  case MUL:
994  {
995 #pragma omp parallel for
996  for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
997  {
998  typename RVEC::size_type rightbase=rightOffset+(singlerightsample?0:i*sampleSize);
999  const LSCALAR* lpos=left+(leftreset?0:i*substep);
1000  for (typename ResVEC::size_type j=0;j<sampleSize;++j)
1001  {
1002  res[i*sampleSize+resOffset+j]=*lpos*right[rightbase+j];
1003  }
1004  }
1005  }
1006  break;
1007  case DIV:
1008  {
1009 #pragma omp parallel for
1010  for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
1011  {
1012  typename RVEC::size_type rightbase=rightOffset+(singlerightsample?0:i*sampleSize);
1013  const LSCALAR* lpos=left+(leftreset?0:i*substep);
1014  for (typename ResVEC::size_type j=0;j<sampleSize;++j)
1015  {
1016  res[i*sampleSize+resOffset+j]=*lpos/right[rightbase+j];
1017  }
1018  }
1019  }
1020  break;
1021  default:
1022  throw DataException("Unsupported binary operation");
1023  }
1024 }
1025 
1026 template <>
1027 void
1028 binaryOpVectorLeftScalar(DataTypes::RealVectorType& res, // where result is to be stored
1029  typename DataTypes::RealVectorType::size_type resOffset, // offset in the result vector to start storing results
1030  const typename DataTypes::RealVectorType::size_type samplesToProcess, // number of samples to be updated in the result
1031  const typename DataTypes::RealVectorType::size_type sampleSize, // number of values in each sample
1032  const DataTypes::real_t* left, // LHS of calculation
1033  const bool leftreset, // true if LHS is providing a single sample of 1 value only
1034  const DataTypes::RealVectorType& right, // RHS of the calculation
1035  typename DataTypes::RealVectorType::size_type rightOffset, // where to start reading RHS values
1036  escript::ES_optype operation, // operation to perform
1037  bool singlerightsample); // right consists of a single sample
1038 
1042 template <class ResVEC, class LVEC, class RVEC>
1043 void
1044 binaryOpVector(ResVEC& res, // where result is to be stored
1045  typename ResVEC::size_type resOffset, // offset in the result vector to start storing results
1046  const typename ResVEC::size_type samplesToProcess, // number of samples to be updated in the result
1047  const typename ResVEC::size_type sampleSize, // number of values in each sample
1048  const LVEC& left, // LHS of calculation
1049  typename LVEC::size_type leftOffset, // where to start reading LHS values
1050  const bool leftreset, // Is LHS only supplying a single sample instead of a bunch of them
1051  const RVEC& right, // RHS of the calculation
1052  typename RVEC::size_type rightOffset, // where to start reading RHS values
1053  const bool rightreset, // Is RHS only supplying a single sample instead of a bunch of them
1054  escript::ES_optype operation) // operation to perform
1055 {
1056  switch (operation)
1057  {
1058  case ADD:
1059  {
1060 #pragma omp parallel for
1061  for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
1062  {
1063  typename LVEC::size_type leftbase=leftOffset+(leftreset?0:i*sampleSize);
1064  typename RVEC::size_type rightbase=rightOffset+(rightreset?0:i*sampleSize);
1065  for (typename ResVEC::size_type j=0;j<sampleSize;++j)
1066  {
1067  res[i*sampleSize+resOffset+j]=left[leftbase+j]+right[rightbase+j];
1068  }
1069  }
1070  }
1071  break;
1072  case POW:
1073  {
1074 #pragma omp parallel for
1075  for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
1076  {
1077  typename LVEC::size_type leftbase=leftOffset+(leftreset?0:i*sampleSize);
1078  typename RVEC::size_type rightbase=rightOffset+(rightreset?0:i*sampleSize);
1079  for (typename ResVEC::size_type j=0;j<sampleSize;++j)
1080  {
1081  res[i*sampleSize+resOffset+j]=pow(left[leftbase+j],right[rightbase+j]);
1082  }
1083  }
1084  }
1085  break;
1086  case SUB:
1087  {
1088 #pragma omp parallel for
1089  for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
1090  {
1091  typename LVEC::size_type leftbase=leftOffset+(leftreset?0:i*sampleSize);
1092  typename RVEC::size_type rightbase=rightOffset+(rightreset?0:i*sampleSize);
1093  for (typename ResVEC::size_type j=0;j<sampleSize;++j)
1094  {
1095  res[i*sampleSize+resOffset+j]=left[leftbase+j]-right[rightbase+j];
1096  }
1097  }
1098  }
1099  break;
1100  case MUL:
1101  {
1102 #pragma omp parallel for
1103  for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
1104  {
1105  typename LVEC::size_type leftbase=leftOffset+(leftreset?0:i*sampleSize);
1106  typename RVEC::size_type rightbase=rightOffset+(rightreset?0:i*sampleSize);
1107  for (typename ResVEC::size_type j=0;j<sampleSize;++j)
1108  {
1109  res[i*sampleSize+resOffset+j]=left[leftbase+j]*right[rightbase+j];
1110  }
1111  }
1112  }
1113  break;
1114  case DIV:
1115  {
1116 #pragma omp parallel for
1117  for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
1118  {
1119  typename LVEC::size_type leftbase=leftOffset+(leftreset?0:i*sampleSize);
1120  typename RVEC::size_type rightbase=rightOffset+(rightreset?0:i*sampleSize);
1121  for (typename ResVEC::size_type j=0;j<sampleSize;++j)
1122  {
1123  res[i*sampleSize+resOffset+j]=left[leftbase+j]/right[rightbase+j];
1124  }
1125  }
1126  }
1127  break;
1128  default:
1129  throw DataException("Unsupported binary operation");
1130  }
1131 }
1132 
1133 template <>
1134 void
1135 binaryOpVector(DataTypes::RealVectorType& res, // where result is to be stored
1136  typename DataTypes::RealVectorType::size_type resOffset, // offset in the result vector to start storing results
1137  const typename DataTypes::RealVectorType::size_type samplesToProcess, // number of samples to be updated in the result
1138  const typename DataTypes::RealVectorType::size_type sampleSize, // number of values in each sample
1139  const DataTypes::RealVectorType& left, // LHS of calculation
1140  typename DataTypes::RealVectorType::size_type leftOffset, // where to start reading LHS values
1141  const bool leftreset, // Is LHS only supplying a single sample instead of a bunch of them
1142  const DataTypes::RealVectorType& right, // RHS of the calculation
1143  typename DataTypes::RealVectorType::size_type rightOffset, // where to start reading RHS values
1144  const bool rightreset, // Is RHS only supplying a single sample instead of a bunch of them
1145  escript::ES_optype operation); // operation to perform
1146 
1147 #define OPVECLAZYBODY(X) \
1148  for (size_t j=0;j<onumsteps;++j)\
1149  {\
1150  for (size_t i=0;i<numsteps;++i,res+=resultStep) \
1151  { \
1152  for (size_t s=0; s<chunksize; ++s)\
1153  {\
1154  res[s] = X;\
1155  }\
1156 /* tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X);*/ \
1157  lroffset+=leftstep; \
1158  rroffset+=rightstep; \
1159  }\
1160  lroffset+=oleftstep;\
1161  rroffset+=orightstep;\
1162  }
1163 
1170 template <class ResELT, class LELT, class RELT>
1171 void
1173  const LELT* left,
1174  const RELT* right,
1175  const size_t chunksize,
1176  const size_t onumsteps,
1177  const size_t numsteps,
1178  const size_t resultStep,
1179  const size_t leftstep,
1180  const size_t rightstep,
1181  const size_t oleftstep,
1182  const size_t orightstep,
1183  size_t lroffset,
1184  size_t rroffset,
1185  escript::ES_optype operation) // operation to perform
1186 {
1187  switch (operation)
1188  {
1189  case ADD:
1190  OPVECLAZYBODY((left[lroffset+s]+right[rroffset+s]));
1191  break;
1192  case POW:
1193  OPVECLAZYBODY(pow(left[lroffset+s],right[rroffset+s]))
1194  break;
1195  case SUB:
1196  OPVECLAZYBODY(left[lroffset+s]-right[rroffset+s])
1197  break;
1198  case MUL:
1199  OPVECLAZYBODY(left[lroffset+s]*right[rroffset+s])
1200  break;
1201  case DIV:
1202  OPVECLAZYBODY(left[lroffset+s]/right[rroffset+s])
1203  break;
1204  case LESS:
1205  OPVECLAZYBODY(left[lroffset+s]<right[rroffset+s])
1206  break;
1207  case GREATER:
1208  OPVECLAZYBODY(left[lroffset+s]>right[rroffset+s])
1209  break;
1210  case GREATER_EQUAL:
1211  OPVECLAZYBODY(left[lroffset+s]>=right[rroffset+s])
1212  break;
1213  case LESS_EQUAL:
1214  OPVECLAZYBODY(left[lroffset+s]<=right[rroffset+s])
1215  break;
1216  default:
1217  ESYS_ASSERT(false, "Invalid operation. This should never happen!");
1218  // I can't throw here because this will be called inside a parallel section
1219  }
1220 }
1221 
1225 /* trying to make a single version for all Tagged+Expanded interactions */
1226 template <class ResVEC, class LVEC, class RVEC>
1227 void
1228 binaryOpVectorTagged(ResVEC& res, // where result is to be stored
1229  const typename ResVEC::size_type samplesToProcess, // number of samples to be updated in the result
1230  const typename ResVEC::size_type DPPSample, // number of datapoints per sample
1231  const typename ResVEC::size_type DPSize, // datapoint size
1232  const LVEC& left, // LHS of calculation
1233  bool leftscalar,
1234  const RVEC& right, // RHS of the calculation
1235  bool rightscalar,
1236  bool lefttagged, // true if left object is the tagged one
1237  const DataTagged& tagsource, // where to get tag offsets from
1238  escript::ES_optype operation) // operation to perform
1239 {
1240  typename ResVEC::size_type lstep=leftscalar?1:DPSize;
1241  typename ResVEC::size_type rstep=rightscalar?1:DPSize;
1242  typename ResVEC::size_type limit=samplesToProcess*DPPSample;
1243  switch (operation)
1244  {
1245  case ADD:
1246  {
1247 #pragma omp parallel for
1248  for (typename ResVEC::size_type i=0;i<limit;++i)
1249  {
1250  typename LVEC::size_type leftbase=(lefttagged?tagsource.getPointOffset(i/DPPSample,0):i*lstep); // only one of these
1251  typename RVEC::size_type rightbase=(lefttagged?i*rstep:tagsource.getPointOffset(i/DPPSample,0)); // will apply
1252 
1253  for (typename ResVEC::size_type j=0;j<DPSize;++j)
1254  {
1255  res[i*DPSize+j]=left[leftbase+j*(!leftscalar)]+right[rightbase+j*(!rightscalar)];
1256  }
1257  }
1258  }
1259  break;
1260  case POW:
1261  {
1262 #pragma omp parallel for
1263  for (typename ResVEC::size_type i=0;i<limit;++i)
1264  {
1265  typename LVEC::size_type leftbase=(lefttagged?tagsource.getPointOffset(i/DPPSample,0):i*lstep); // only one of these
1266  typename RVEC::size_type rightbase=(lefttagged?i*rstep:tagsource.getPointOffset(i/DPPSample,0)); // will apply
1267 
1268  for (typename ResVEC::size_type j=0;j<DPSize;++j)
1269  {
1270  res[i*DPSize+j]=pow(left[leftbase+j*(!leftscalar)],right[rightbase+j*(!rightscalar)]);
1271  }
1272  }
1273  }
1274  break;
1275  case SUB:
1276  {
1277 #pragma omp parallel for
1278  for (typename ResVEC::size_type i=0;i<limit;++i)
1279  {
1280  typename LVEC::size_type leftbase=(lefttagged?tagsource.getPointOffset(i/DPPSample,0):i*lstep); // only one of these
1281  typename RVEC::size_type rightbase=(lefttagged?i*rstep:tagsource.getPointOffset(i/DPPSample,0)); // will apply
1282 
1283  for (typename ResVEC::size_type j=0;j<DPSize;++j)
1284  {
1285  res[i*DPSize+j]=left[leftbase+j*(!leftscalar)]-right[rightbase+j*(!rightscalar)];
1286  }
1287  }
1288  }
1289  break;
1290  case MUL:
1291  {
1292 #pragma omp parallel for
1293  for (typename ResVEC::size_type i=0;i<limit;++i)
1294  {
1295  typename LVEC::size_type leftbase=(lefttagged?tagsource.getPointOffset(i/DPPSample,0):i*lstep); // only one of these
1296  typename RVEC::size_type rightbase=(lefttagged?i*rstep:tagsource.getPointOffset(i/DPPSample,0)); // will apply
1297 
1298  for (typename ResVEC::size_type j=0;j<DPSize;++j)
1299  {
1300  res[i*DPSize+j]=left[leftbase+j*(!leftscalar)]*right[rightbase+j*(!rightscalar)];
1301  }
1302  }
1303  }
1304  break;
1305  case DIV:
1306  {
1307 #pragma omp parallel for
1308  for (typename ResVEC::size_type i=0;i<limit;++i)
1309  {
1310  typename LVEC::size_type leftbase=(lefttagged?tagsource.getPointOffset(i/DPPSample,0):i*lstep); // only one of these
1311  typename RVEC::size_type rightbase=(lefttagged?i*rstep:tagsource.getPointOffset(i/DPPSample,0)); // will apply
1312 
1313  for (typename ResVEC::size_type j=0;j<DPSize;++j)
1314  {
1315  res[i*DPSize+j]=left[leftbase+j*(!leftscalar)]/right[rightbase+j*(!rightscalar)];
1316  }
1317  }
1318  }
1319  break;
1320  default:
1321  throw DataException("Unsupported binary operation");
1322  }
1323 }
1324 
1325 template<>
1326 void
1327 binaryOpVectorTagged(DataTypes::RealVectorType& res, // where result is to be stored
1328  const typename DataTypes::RealVectorType::size_type samplesToProcess, // number of samples to be updated in the result
1329  const typename DataTypes::RealVectorType::size_type DPPSample, // number of datapoints per sample
1330  const typename DataTypes::RealVectorType::size_type DPSize, // datapoint size
1331  const DataTypes::RealVectorType& left, // LHS of calculation
1332  const bool leftscalar,
1333  const DataTypes::RealVectorType& right, // RHS of the calculation
1334  const bool rightscalar,
1335  const bool lefttagged, // true if left object is the tagged one
1336  const DataTagged& tagsource, // where to get tag offsets from
1337  escript::ES_optype operation);
1338 
1339 
1340 
1341 
1358 template <class BinaryFunction>
1359 inline
1362  const DataTypes::ShapeType& leftShape,
1364  BinaryFunction operation,
1365  DataTypes::real_t initial_value)
1366 {
1367  ESYS_ASSERT((left.size()>0)&&checkOffset(left,leftShape,offset),
1368  "Couldn't perform reductionOp due to insufficient storage.");
1369  DataTypes::real_t current_value=initial_value;
1370  for (DataTypes::RealVectorType::size_type i=0;i<DataTypes::noValues(leftShape);i++) {
1371  current_value=operation(current_value,left[offset+i]);
1372  }
1373  return current_value;
1374 }
1375 
1376 template <class BinaryFunction>
1377 inline
1380  const DataTypes::ShapeType& leftShape,
1382  BinaryFunction operation,
1383  DataTypes::real_t initial_value)
1384 {
1385  ESYS_ASSERT((left.size()>0)&&checkOffset(left,leftShape,offset),
1386  "Couldn't perform reductionOp due to insufficient storage.");
1387  DataTypes::real_t current_value=initial_value;
1388  for (DataTypes::RealVectorType::size_type i=0;i<DataTypes::noValues(leftShape);i++) {
1389  current_value=operation(current_value,left[offset+i]);
1390  }
1391  return current_value;
1392 }
1393 
1394 
1411 int
1413  const DataTypes::ShapeType& inShape,
1416  const DataTypes::ShapeType& outShape,
1418  int count,
1419  LapackInverseHelper& helper);
1420 
1428 void
1429 matrixInverseError(int err);
1430 
1435 inline
1436 bool
1438 {
1439  for (size_t z=inOffset;z<inOffset+count;++z)
1440  {
1441  if (nancheck(in[z]))
1442  {
1443  return true;
1444  }
1445  }
1446  return false;
1447 }
1448 
1449 } // end namespace escript
1450 
1451 #endif // __ESCRIPT_DATAMATHS_H__
1452 
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 binaryOpVectorLazyHelper(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:1172
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
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:148
virtual DataTypes::RealVectorType::size_type getPointOffset(int sampleNo, int dataPointNo) const
getPointOffset
Definition: DataTagged.cpp:840
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:438
bool vectorHasNaN(const DataTypes::RealVectorType &in, DataTypes::RealVectorType::size_type inOffset, size_t count)
returns true if the vector contains NaN
Definition: DataVectorOps.h:1437
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
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:198
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:170
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:1147
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:329
bool checkOffset(const VEC &data, const DataTypes::ShapeType &shape, typename VEC::size_type offset)
Definition: DataVectorOps.h:813
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:237
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:38
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:1361
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:746
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:194