escript  Revision_
DataMaths.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 
18 #ifndef escript_DataMaths_20080822_H
19 #define escript_DataMaths_20080822_H
20 #include "DataAbstract.h"
21 #include "DataException.h"
22 #include "LocalOps.h"
23 #include "LapackInverseHelper.h"
24 
35 namespace escript
36 {
37 namespace DataMaths
38 {
39 
63  template <class UnaryFunction>
64  void
67  UnaryFunction operation);
68 
83  template <class BinaryFunction>
84  void
86  const DataTypes::ShapeType& leftShape,
88  const DataTypes::ValueType& right,
89  const DataTypes::ShapeType& rightShape,
91  BinaryFunction operation);
92 
109  template <class BinaryFunction>
110  void
112  const DataTypes::ShapeType& shape,
114  double right,
115  BinaryFunction operation);
116 
133  template <class BinaryFunction>
134  double
135  reductionOp(const DataTypes::ValueType& left,
136  const DataTypes::ShapeType& shape,
138  BinaryFunction operation,
139  double initial_value);
140 
156  void
157  matMult(const DataTypes::ValueType& left,
158  const DataTypes::ShapeType& leftShape,
160  const DataTypes::ValueType& right,
161  const DataTypes::ShapeType& rightShape,
163  DataTypes::ValueType& result,
164  const DataTypes::ShapeType& resultShape);
165 // Hmmmm why is there no offset for the result??
166 
167 
168 
169 
181  const DataTypes::ShapeType& right);
182 
195  inline
196  void
198  const DataTypes::ShapeType& inShape,
201  const DataTypes::ShapeType& evShape,
203  {
204  if (DataTypes::getRank(inShape) == 2) {
205  int i0, i1;
206  int s0=inShape[0];
207  int s1=inShape[1];
208  for (i0=0; i0<s0; i0++) {
209  for (i1=0; i1<s1; i1++) {
210  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = (in[inOffset+DataTypes::getRelIndex(inShape,i0,i1)] + in[inOffset+DataTypes::getRelIndex(inShape,i1,i0)]) / 2.0;
211  }
212  }
213  }
214  else if (DataTypes::getRank(inShape) == 4) {
215  int i0, i1, i2, i3;
216  int s0=inShape[0];
217  int s1=inShape[1];
218  int s2=inShape[2];
219  int s3=inShape[3];
220  for (i0=0; i0<s0; i0++) {
221  for (i1=0; i1<s1; i1++) {
222  for (i2=0; i2<s2; i2++) {
223  for (i3=0; i3<s3; i3++) {
224  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;
225  }
226  }
227  }
228  }
229  }
230  }
231 
244  inline
245  void
247  const DataTypes::ShapeType& inShape,
250  const DataTypes::ShapeType& evShape,
252  {
253  if (DataTypes::getRank(inShape) == 2) {
254  int i0, i1;
255  int s0=inShape[0];
256  int s1=inShape[1];
257  for (i0=0; i0<s0; i0++) {
258  for (i1=0; i1<s1; i1++) {
259  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = (in[inOffset+DataTypes::getRelIndex(inShape,i0,i1)] - in[inOffset+DataTypes::getRelIndex(inShape,i1,i0)]) / 2.0;
260  }
261  }
262  }
263  else if (DataTypes::getRank(inShape) == 4) {
264  int i0, i1, i2, i3;
265  int s0=inShape[0];
266  int s1=inShape[1];
267  int s2=inShape[2];
268  int s3=inShape[3];
269  for (i0=0; i0<s0; i0++) {
270  for (i1=0; i1<s1; i1++) {
271  for (i2=0; i2<s2; i2++) {
272  for (i3=0; i3<s3; i3++) {
273  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;
274  }
275  }
276  }
277  }
278  }
279  }
280 
293  inline
294  void
296  const DataTypes::ShapeType& inShape,
299  const DataTypes::ShapeType& evShape,
301  int axis_offset)
302  {
303  for (int j=0;j<DataTypes::noValues(evShape);++j)
304  {
305  ev[evOffset+j]=0;
306  }
307  if (DataTypes::getRank(inShape) == 2) {
308  int s0=inShape[0]; // Python wrapper limits to square matrix
309  int i;
310  for (i=0; i<s0; i++) {
311  ev[evOffset/*+DataTypes::getRelIndex(evShape)*/] += in[inOffset+DataTypes::getRelIndex(inShape,i,i)];
312  }
313  }
314  else if (DataTypes::getRank(inShape) == 3) {
315  if (axis_offset==0) {
316  int s0=inShape[0];
317  int s2=inShape[2];
318  int i0, i2;
319  for (i0=0; i0<s0; i0++) {
320  for (i2=0; i2<s2; i2++) {
321  ev[evOffset+DataTypes::getRelIndex(evShape,i2)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i0,i2)];
322  }
323  }
324  }
325  else if (axis_offset==1) {
326  int s0=inShape[0];
327  int s1=inShape[1];
328  int i0, i1;
329  for (i0=0; i0<s0; i0++) {
330  for (i1=0; i1<s1; i1++) {
331  ev[evOffset+DataTypes::getRelIndex(evShape,i0)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i1)];
332  }
333  }
334  }
335  }
336  else if (DataTypes::getRank(inShape) == 4) {
337  if (axis_offset==0) {
338  int s0=inShape[0];
339  int s2=inShape[2];
340  int s3=inShape[3];
341  int i0, i2, i3;
342  for (i0=0; i0<s0; i0++) {
343  for (i2=0; i2<s2; i2++) {
344  for (i3=0; i3<s3; i3++) {
345  ev[evOffset+DataTypes::getRelIndex(evShape,i2,i3)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i0,i2,i3)];
346  }
347  }
348  }
349  }
350  else if (axis_offset==1) {
351  int s0=inShape[0];
352  int s1=inShape[1];
353  int s3=inShape[3];
354  int i0, i1, i3;
355  for (i0=0; i0<s0; i0++) {
356  for (i1=0; i1<s1; i1++) {
357  for (i3=0; i3<s3; i3++) {
358  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i3)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i1,i3)];
359  }
360  }
361  }
362  }
363  else if (axis_offset==2) {
364  int s0=inShape[0];
365  int s1=inShape[1];
366  int s2=inShape[2];
367  int i0, i1, i2;
368  for (i0=0; i0<s0; i0++) {
369  for (i1=0; i1<s1; i1++) {
370  for (i2=0; i2<s2; i2++) {
371  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i2,i2)];
372  }
373  }
374  }
375  }
376  }
377  }
378 
392  inline
393  void
395  const DataTypes::ShapeType& inShape,
398  const DataTypes::ShapeType& evShape,
400  int axis_offset)
401  {
402  int inRank=DataTypes::getRank(inShape);
403  if ( inRank== 4) {
404  int s0=evShape[0];
405  int s1=evShape[1];
406  int s2=evShape[2];
407  int s3=evShape[3];
408  int i0, i1, i2, i3;
409  if (axis_offset==1) {
410  for (i0=0; i0<s0; i0++) {
411  for (i1=0; i1<s1; i1++) {
412  for (i2=0; i2<s2; i2++) {
413  for (i3=0; i3<s3; i3++) {
414  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i3,i0,i1,i2)];
415  }
416  }
417  }
418  }
419  }
420  else if (axis_offset==2) {
421  for (i0=0; i0<s0; i0++) {
422  for (i1=0; i1<s1; i1++) {
423  for (i2=0; i2<s2; i2++) {
424  for (i3=0; i3<s3; i3++) {
425  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i3,i0,i1)];
426  }
427  }
428  }
429  }
430  }
431  else if (axis_offset==3) {
432  for (i0=0; i0<s0; i0++) {
433  for (i1=0; i1<s1; i1++) {
434  for (i2=0; i2<s2; i2++) {
435  for (i3=0; i3<s3; i3++) {
436  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i2,i3,i0)];
437  }
438  }
439  }
440  }
441  }
442  else {
443  for (i0=0; i0<s0; i0++) {
444  for (i1=0; i1<s1; i1++) {
445  for (i2=0; i2<s2; i2++) {
446  for (i3=0; i3<s3; i3++) {
447  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i2,i3)];
448  }
449  }
450  }
451  }
452  }
453  }
454  else if (inRank == 3) {
455  int s0=evShape[0];
456  int s1=evShape[1];
457  int s2=evShape[2];
458  int i0, i1, i2;
459  if (axis_offset==1) {
460  for (i0=0; i0<s0; i0++) {
461  for (i1=0; i1<s1; i1++) {
462  for (i2=0; i2<s2; i2++) {
463  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i0,i1)];
464  }
465  }
466  }
467  }
468  else if (axis_offset==2) {
469  for (i0=0; i0<s0; i0++) {
470  for (i1=0; i1<s1; i1++) {
471  for (i2=0; i2<s2; i2++) {
472  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i2,i0)];
473  }
474  }
475  }
476  }
477  else {
478  // Copy the matrix unchanged
479  for (i0=0; i0<s0; i0++) {
480  for (i1=0; i1<s1; i1++) {
481  for (i2=0; i2<s2; i2++) {
482  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i2)];
483  }
484  }
485  }
486  }
487  }
488  else if (inRank == 2) {
489  int s0=evShape[0];
490  int s1=evShape[1];
491  int i0, i1;
492  if (axis_offset==1) {
493  for (i0=0; i0<s0; i0++) {
494  for (i1=0; i1<s1; i1++) {
495  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0)];
496  }
497  }
498  }
499  else {
500  for (i0=0; i0<s0; i0++) {
501  for (i1=0; i1<s1; i1++) {
502  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1)];
503  }
504  }
505  }
506  }
507  else if (inRank == 1) {
508  int s0=evShape[0];
509  int i0;
510  for (i0=0; i0<s0; i0++) {
511  ev[evOffset+DataTypes::getRelIndex(evShape,i0)] = in[inOffset+DataTypes::getRelIndex(inShape,i0)];
512  }
513  }
514  else if (inRank == 0) {
515  ev[evOffset/*+DataTypes::getRelIndex(evShape,)*/] = in[inOffset/*+DataTypes::getRelIndex(inShape,)*/];
516  }
517  else {
518  throw DataException("Error - DataArrayView::transpose can only be calculated for rank 0, 1, 2, 3 or 4 objects.");
519  }
520  }
521 
536  inline
537  void
539  const DataTypes::ShapeType& inShape,
542  const DataTypes::ShapeType& evShape,
544  int axis0,
545  int axis1)
546  {
547  int inRank=DataTypes::getRank(inShape);
548  if (inRank == 4) {
549  int s0=evShape[0];
550  int s1=evShape[1];
551  int s2=evShape[2];
552  int s3=evShape[3];
553  int i0, i1, i2, i3;
554  if (axis0==0) {
555  if (axis1==1) {
556  for (i0=0; i0<s0; i0++) {
557  for (i1=0; i1<s1; i1++) {
558  for (i2=0; i2<s2; i2++) {
559  for (i3=0; i3<s3; i3++) {
560  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0,i2,i3)];
561  }
562  }
563  }
564  }
565  } else if (axis1==2) {
566  for (i0=0; i0<s0; i0++) {
567  for (i1=0; i1<s1; i1++) {
568  for (i2=0; i2<s2; i2++) {
569  for (i3=0; i3<s3; i3++) {
570  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i1,i0,i3)];
571  }
572  }
573  }
574  }
575 
576  } else if (axis1==3) {
577  for (i0=0; i0<s0; i0++) {
578  for (i1=0; i1<s1; i1++) {
579  for (i2=0; i2<s2; i2++) {
580  for (i3=0; i3<s3; i3++) {
581  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i3,i1,i2,i0)];
582  }
583  }
584  }
585  }
586  }
587  } else if (axis0==1) {
588  if (axis1==2) {
589  for (i0=0; i0<s0; i0++) {
590  for (i1=0; i1<s1; i1++) {
591  for (i2=0; i2<s2; i2++) {
592  for (i3=0; i3<s3; i3++) {
593  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i2,i1,i3)];
594  }
595  }
596  }
597  }
598  } else if (axis1==3) {
599  for (i0=0; i0<s0; i0++) {
600  for (i1=0; i1<s1; i1++) {
601  for (i2=0; i2<s2; i2++) {
602  for (i3=0; i3<s3; i3++) {
603  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i3,i2,i1)];
604  }
605  }
606  }
607  }
608  }
609  } else if (axis0==2) {
610  if (axis1==3) {
611  for (i0=0; i0<s0; i0++) {
612  for (i1=0; i1<s1; i1++) {
613  for (i2=0; i2<s2; i2++) {
614  for (i3=0; i3<s3; i3++) {
615  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i3,i2)];
616  }
617  }
618  }
619  }
620  }
621  }
622 
623  } else if ( inRank == 3) {
624  int s0=evShape[0];
625  int s1=evShape[1];
626  int s2=evShape[2];
627  int i0, i1, i2;
628  if (axis0==0) {
629  if (axis1==1) {
630  for (i0=0; i0<s0; i0++) {
631  for (i1=0; i1<s1; i1++) {
632  for (i2=0; i2<s2; i2++) {
633  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0,i2)];
634  }
635  }
636  }
637  } else if (axis1==2) {
638  for (i0=0; i0<s0; i0++) {
639  for (i1=0; i1<s1; i1++) {
640  for (i2=0; i2<s2; i2++) {
641  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i1,i0)];
642  }
643  }
644  }
645  }
646  } else if (axis0==1) {
647  if (axis1==2) {
648  for (i0=0; i0<s0; i0++) {
649  for (i1=0; i1<s1; i1++) {
650  for (i2=0; i2<s2; i2++) {
651  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i2,i1)];
652  }
653  }
654  }
655  }
656  }
657  } else if ( inRank == 2) {
658  int s0=evShape[0];
659  int s1=evShape[1];
660  int i0, i1;
661  if (axis0==0) {
662  if (axis1==1) {
663  for (i0=0; i0<s0; i0++) {
664  for (i1=0; i1<s1; i1++) {
665  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0)];
666  }
667  }
668  }
669  }
670  } else {
671  throw DataException("Error - DataArrayView::swapaxes can only be calculated for rank 2, 3 or 4 objects.");
672  }
673  }
674 
687  inline
688  void
690  const DataTypes::ShapeType& inShape,
693  const DataTypes::ShapeType& evShape,
695  {
696  double in00,in10,in20,in01,in11,in21,in02,in12,in22;
697  double ev0,ev1,ev2;
698  int s=inShape[0];
699  if (s==1) {
700  in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
701  eigenvalues1(in00,&ev0);
702  ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
703 
704  } else if (s==2) {
705  in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
706  in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
707  in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
708  in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
709  eigenvalues2(in00,(in01+in10)/2.,in11,&ev0,&ev1);
710  ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
711  ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
712 
713  } else if (s==3) {
714  in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
715  in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
716  in20=in[inOffset+DataTypes::getRelIndex(inShape,2,0)];
717  in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
718  in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
719  in21=in[inOffset+DataTypes::getRelIndex(inShape,2,1)];
720  in02=in[inOffset+DataTypes::getRelIndex(inShape,0,2)];
721  in12=in[inOffset+DataTypes::getRelIndex(inShape,1,2)];
722  in22=in[inOffset+DataTypes::getRelIndex(inShape,2,2)];
723  eigenvalues3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
724  &ev0,&ev1,&ev2);
725  ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
726  ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
727  ev[evOffset+DataTypes::getRelIndex(evShape,2)]=ev2;
728 
729  }
730  }
731 
748  inline
749  void
752  DataTypes::ValueType& ev, const DataTypes::ShapeType& evShape,
756  const double tol=1.e-13)
757  {
758  double in00,in10,in20,in01,in11,in21,in02,in12,in22;
759  double V00,V10,V20,V01,V11,V21,V02,V12,V22;
760  double ev0,ev1,ev2;
761  int s=inShape[0];
762  if (s==1) {
763  in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
764  eigenvalues_and_eigenvectors1(in00,&ev0,&V00,tol);
765  ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
766  V[inOffset+DataTypes::getRelIndex(VShape,0,0)]=V00;
767  } else if (s==2) {
768  in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
769  in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
770  in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
771  in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
772  eigenvalues_and_eigenvectors2(in00,(in01+in10)/2.,in11,
773  &ev0,&ev1,&V00,&V10,&V01,&V11,tol);
774  ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
775  ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
776  V[inOffset+DataTypes::getRelIndex(VShape,0,0)]=V00;
777  V[inOffset+DataTypes::getRelIndex(VShape,1,0)]=V10;
778  V[inOffset+DataTypes::getRelIndex(VShape,0,1)]=V01;
779  V[inOffset+DataTypes::getRelIndex(VShape,1,1)]=V11;
780  } else if (s==3) {
781  in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
782  in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
783  in20=in[inOffset+DataTypes::getRelIndex(inShape,2,0)];
784  in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
785  in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
786  in21=in[inOffset+DataTypes::getRelIndex(inShape,2,1)];
787  in02=in[inOffset+DataTypes::getRelIndex(inShape,0,2)];
788  in12=in[inOffset+DataTypes::getRelIndex(inShape,1,2)];
789  in22=in[inOffset+DataTypes::getRelIndex(inShape,2,2)];
790  eigenvalues_and_eigenvectors3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
791  &ev0,&ev1,&ev2,
792  &V00,&V10,&V20,&V01,&V11,&V21,&V02,&V12,&V22,tol);
793  ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
794  ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
795  ev[evOffset+DataTypes::getRelIndex(evShape,2)]=ev2;
796  V[inOffset+DataTypes::getRelIndex(VShape,0,0)]=V00;
797  V[inOffset+DataTypes::getRelIndex(VShape,1,0)]=V10;
798  V[inOffset+DataTypes::getRelIndex(VShape,2,0)]=V20;
799  V[inOffset+DataTypes::getRelIndex(VShape,0,1)]=V01;
800  V[inOffset+DataTypes::getRelIndex(VShape,1,1)]=V11;
801  V[inOffset+DataTypes::getRelIndex(VShape,2,1)]=V21;
802  V[inOffset+DataTypes::getRelIndex(VShape,0,2)]=V02;
803  V[inOffset+DataTypes::getRelIndex(VShape,1,2)]=V12;
804  V[inOffset+DataTypes::getRelIndex(VShape,2,2)]=V22;
805 
806  }
807  }
808 
809 
814 inline
815 bool
817  const DataTypes::ShapeType& shape,
819 {
820  return (data.size() >= (offset+DataTypes::noValues(shape)));
821 }
822 
823 template <class UnaryFunction>
824 inline
825 void
828  UnaryFunction operation)
829 {
830  EsysAssert((data.size()>0)&&checkOffset(data,shape,offset),
831  "Error - Couldn't perform unaryOp due to insufficient storage.");
833  for (DataTypes::ValueType::size_type i=0;i<nVals;i++) {
834  data[offset+i]=operation(data[offset+i]);
835  }
836 }
837 
838 
839 template <class BinaryFunction>
840 inline
841 void
843  const DataTypes::ShapeType& leftShape,
845  const DataTypes::ValueType& right,
846  const DataTypes::ShapeType& rightShape,
848  BinaryFunction operation)
849 {
850  EsysAssert(leftShape==rightShape,
851  "Error - Couldn't perform binaryOp due to shape mismatch,");
852  EsysAssert(((left.size()>0)&&checkOffset(left,leftShape, leftOffset)),
853  "Error - Couldn't perform binaryOp due to insufficient storage in left object.");
854  EsysAssert(((right.size()>0)&&checkOffset(right,rightShape,rightOffset)),
855  "Error - Couldn't perform binaryOp due to insufficient storage in right object.");
856  for (DataTypes::ValueType::size_type i=0;i<DataTypes::noValues(leftShape);i++) {
857  left[leftOffset+i]=operation(left[leftOffset+i],right[rightOffset+i]);
858  }
859 }
860 
861 template <class BinaryFunction>
862 inline
863 void
865  const DataTypes::ShapeType& leftShape,
867  double right,
868  BinaryFunction operation)
869 {
870  EsysAssert(((left.size()>0)&&checkOffset(left,leftShape,offset)),
871  "Error - Couldn't perform binaryOp due to insufficient storage in left object.");
872  for (DataTypes::ValueType::size_type i=0;i<DataTypes::noValues(leftShape);i++) {
873  left[offset+i]=operation(left[offset+i],right);
874  }
875 }
876 
877 template <class BinaryFunction>
878 inline
879 double
881  const DataTypes::ShapeType& leftShape,
883  BinaryFunction operation,
884  double initial_value)
885 {
886  EsysAssert(((left.size()>0)&&checkOffset(left,leftShape,offset)),
887  "Error - Couldn't perform reductionOp due to insufficient storage.");
888  double current_value=initial_value;
889  for (DataTypes::ValueType::size_type i=0;i<DataTypes::noValues(leftShape);i++) {
890  current_value=operation(current_value,left[offset+i]);
891  }
892  return current_value;
893 }
894 
911 int
913  const DataTypes::ShapeType& inShape,
916  const DataTypes::ShapeType& outShape,
918  int count,
919  LapackInverseHelper& helper);
920 
928 void
929 matrixInverseError(int err);
930 
935 inline
936 bool
938 {
939  for (size_t z=inOffset;z<inOffset+count;++z)
940  {
941  if (nancheck(in[z]))
942  {
943  return true;
944  }
945  }
946  return false;
947 }
948 
949 } // end namespace DataMath
950 } // end namespace escript
951 #endif
952 
DataVector implements an arbitrarily long vector of data values. DataVector is the underlying data co...
Definition: DataVector.h:44
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: DataMaths.cpp:173
void matMult(const DataTypes::ValueType &left, const DataTypes::ShapeType &leftShape, DataTypes::ValueType::size_type leftOffset, const DataTypes::ValueType &right, const DataTypes::ShapeType &rightShape, DataTypes::ValueType::size_type rightOffset, DataTypes::ValueType &result, const DataTypes::ShapeType &resultShape)
Perform a matrix multiply of the given views.
Definition: DataMaths.cpp:42
escript::DataVector ValueType
Vector to store underlying data.
Definition: DataTypes.h:37
void matrixInverseError(int err)
throws an appropriate exception based on failure of matrix_inverse.
Definition: DataMaths.cpp:189
Definition: AbstractContinuousDomain.cpp:24
void transpose(const DataTypes::ValueType &in, const DataTypes::ShapeType &inShape, DataTypes::ValueType::size_type inOffset, DataTypes::ValueType &ev, const DataTypes::ShapeType &evShape, DataTypes::ValueType::size_type evOffset, int axis_offset)
Transpose each data point of this Data object around the given axis.
Definition: DataMaths.h:394
void binaryOp(DataTypes::ValueType &left, const DataTypes::ShapeType &leftShape, DataTypes::ValueType::size_type leftOffset, const DataTypes::ValueType &right, const DataTypes::ShapeType &rightShape, DataTypes::ValueType::size_type rightOffset, BinaryFunction operation)
Perform the binary operation on the data points specified by the given offsets in the "left" and "rig...
Definition: DataMaths.h:842
size_type size() const
Return the number of elements in this DataVector.
Definition: DataVector.h:215
Definition: LapackInverseHelper.h:26
bool vectorHasNaN(const DataTypes::ValueType &in, DataTypes::ValueType::size_type inOffset, size_t count)
returns true if the vector contains NaN
Definition: DataMaths.h:937
std::vector< int > ShapeType
The shape of a single datapoint.
Definition: DataTypes.h:38
void unaryOp(DataTypes::ValueType &data, const DataTypes::ShapeType &shape, DataTypes::ValueType::size_type offset, UnaryFunction operation)
Perform the unary operation on the data point specified by the given offset. Applies the specified op...
Definition: DataMaths.h:826
void swapaxes(const DataTypes::ValueType &in, const DataTypes::ShapeType &inShape, DataTypes::ValueType::size_type inOffset, DataTypes::ValueType &ev, const DataTypes::ShapeType &evShape, DataTypes::ValueType::size_type evOffset, int axis0, int axis1)
swaps the components axis0 and axis1.
Definition: DataMaths.h:538
void eigenvalues(const DataTypes::ValueType &in, const DataTypes::ShapeType &inShape, DataTypes::ValueType::size_type inOffset, DataTypes::ValueType &ev, const DataTypes::ShapeType &evShape, DataTypes::ValueType::size_type evOffset)
solves a local eigenvalue problem
Definition: DataMaths.h:689
void symmetric(const DataTypes::ValueType &in, const DataTypes::ShapeType &inShape, DataTypes::ValueType::size_type inOffset, DataTypes::ValueType &ev, const DataTypes::ShapeType &evShape, DataTypes::ValueType::size_type evOffset)
computes a symmetric matrix from your square matrix A: (A + transpose(A)) / 2
Definition: DataMaths.h:197
void eigenvalues2(const double A00, const double A01, const double A11, double *ev0, double *ev1)
solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A
Definition: LocalOps.h:94
double reductionOp(const DataTypes::ValueType &left, const DataTypes::ShapeType &shape, DataTypes::ValueType::size_type offset, BinaryFunction operation, double initial_value)
Perform the given data point reduction operation on the data point specified by the given offset into...
Definition: DataMaths.h:880
Describes binary operations performed on double*.
bool nancheck(double d)
acts as a wrapper to isnan.
Definition: LocalOps.h:41
DataTypes::ValueType::size_type getRelIndex(const DataTypes::ShapeType &shape, DataTypes::ValueType::size_type i)
Compute the offset (in 1D vector) of a given subscript with a shape.
Definition: DataTypes.h:183
#define EsysAssert(AssertTest, AssertMessage)
EsysAssert is a MACRO that will throw an exception if the boolean condition specified is false...
Definition: EsysAssert.h:96
void eigenvalues_and_eigenvectors(const DataTypes::ValueType &in, const DataTypes::ShapeType &inShape, DataTypes::ValueType::size_type inOffset, DataTypes::ValueType &ev, const DataTypes::ShapeType &evShape, DataTypes::ValueType::size_type evOffset, DataTypes::ValueType &V, const DataTypes::ShapeType &VShape, DataTypes::ValueType::size_type VOffset, const double tol=1.e-13)
solves a local eigenvalue problem
Definition: DataMaths.h:750
int matrix_inverse(const DataTypes::ValueType &in, const DataTypes::ShapeType &inShape, DataTypes::ValueType::size_type inOffset, DataTypes::ValueType &out, const DataTypes::ShapeType &outShape, DataTypes::ValueType::size_type outOffset, int count, LapackInverseHelper &helper)
computes the inverses of square (up to 3x3) matricies
Definition: DataMaths.cpp:210
void eigenvalues_and_eigenvectors1(const double A00, double *ev0, double *V00, const double tol)
solves a 1x1 eigenvalue A*V=ev*V problem for symmetric A
Definition: LocalOps.h:161
#define V(_K_, _I_)
Definition: ShapeFunctions.cpp:120
void eigenvalues1(const double A00, double *ev0)
solves a 1x1 eigenvalue A*V=ev*V problem
Definition: LocalOps.h:78
int noValues(const ShapeType &shape)
Calculate the number of values in a datapoint with the given shape.
Definition: DataTypes.cpp:94
DataException exception class.
Definition: DataException.h:35
#define ESCRIPT_DLL_API
Definition: escriptcore/src/system_dep.h:54
void eigenvalues_and_eigenvectors2(const double A00, const double A01, const double A11, double *ev0, double *ev1, double *V00, double *V10, double *V01, double *V11, const double tol)
solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A. Eigenvectors are ordered by increasing valu...
Definition: LocalOps.h:253
void nonsymmetric(const DataTypes::ValueType &in, const DataTypes::ShapeType &inShape, DataTypes::ValueType::size_type inOffset, DataTypes::ValueType &ev, const DataTypes::ShapeType &evShape, DataTypes::ValueType::size_type evOffset)
computes a nonsymmetric matrix from your square matrix A: (A - transpose(A)) / 2
Definition: DataMaths.h:246
int getRank(const DataTypes::ShapeType &shape)
Return the rank (number of dimensions) of the given shape.
Definition: DataTypes.h:167
void eigenvalues3(const double A00, const double A01, const double A02, const double A11, const double A12, const double A22, double *ev0, double *ev1, double *ev2)
solves a 3x3 eigenvalue A*V=ev*V problem for symmetric A
Definition: LocalOps.h:118
void eigenvalues_and_eigenvectors3(const double A00, const double A01, const double A02, const double A11, const double A12, const double A22, double *ev0, double *ev1, double *ev2, double *V00, double *V10, double *V20, double *V01, double *V11, double *V21, double *V02, double *V12, double *V22, const double tol)
solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A. Eigenvectors are ordered by increasing valu...
Definition: LocalOps.h:362
long size_type
Definition: DataVector.h:60
void trace(const DataTypes::ValueType &in, const DataTypes::ShapeType &inShape, DataTypes::ValueType::size_type inOffset, DataTypes::ValueType &ev, const DataTypes::ShapeType &evShape, DataTypes::ValueType::size_type evOffset, int axis_offset)
computes the trace of a matrix
Definition: DataMaths.h:295
bool checkOffset(const DataTypes::ValueType &data, const DataTypes::ShapeType &shape, DataTypes::ValueType::size_type offset)
Definition: DataMaths.h:816