ESScript  Revision_4488
DataMaths.h
Go to the documentation of this file.
1 
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2013 by University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development since 2012 by School of Earth Sciences
13 *
14 *****************************************************************************/
15 
16 
17 #ifndef escript_DataMaths_20080822_H
18 #define escript_DataMaths_20080822_H
19 #include "DataAbstract.h"
20 #include "DataException.h"
21 #include "LocalOps.h"
22 #include "LapackInverseHelper.h"
23 
34 namespace escript
35 {
36 namespace DataMaths
37 {
38 
62  template <class UnaryFunction>
63  void
66  UnaryFunction operation);
67 
82  template <class BinaryFunction>
83  void
85  const DataTypes::ShapeType& leftShape,
87  const DataTypes::ValueType& right,
88  const DataTypes::ShapeType& rightShape,
90  BinaryFunction operation);
91 
108  template <class BinaryFunction>
109  void
111  const DataTypes::ShapeType& shape,
113  double right,
114  BinaryFunction operation);
115 
132  template <class BinaryFunction>
133  double
134  reductionOp(const DataTypes::ValueType& left,
135  const DataTypes::ShapeType& shape,
137  BinaryFunction operation,
138  double initial_value);
139 
155  void
156  matMult(const DataTypes::ValueType& left,
157  const DataTypes::ShapeType& leftShape,
159  const DataTypes::ValueType& right,
160  const DataTypes::ShapeType& rightShape,
162  DataTypes::ValueType& result,
163  const DataTypes::ShapeType& resultShape);
164 // Hmmmm why is there no offset for the result??
165 
166 
167 
168 
180  const DataTypes::ShapeType& right);
181 
194  inline
195  void
197  const DataTypes::ShapeType& inShape,
200  const DataTypes::ShapeType& evShape,
202  {
203  if (DataTypes::getRank(inShape) == 2) {
204  int i0, i1;
205  int s0=inShape[0];
206  int s1=inShape[1];
207  for (i0=0; i0<s0; i0++) {
208  for (i1=0; i1<s1; i1++) {
209  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = (in[inOffset+DataTypes::getRelIndex(inShape,i0,i1)] + in[inOffset+DataTypes::getRelIndex(inShape,i1,i0)]) / 2.0;
210  }
211  }
212  }
213  else if (DataTypes::getRank(inShape) == 4) {
214  int i0, i1, i2, i3;
215  int s0=inShape[0];
216  int s1=inShape[1];
217  int s2=inShape[2];
218  int s3=inShape[3];
219  for (i0=0; i0<s0; i0++) {
220  for (i1=0; i1<s1; i1++) {
221  for (i2=0; i2<s2; i2++) {
222  for (i3=0; i3<s3; i3++) {
223  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;
224  }
225  }
226  }
227  }
228  }
229  }
230 
243  inline
244  void
246  const DataTypes::ShapeType& inShape,
249  const DataTypes::ShapeType& evShape,
251  {
252  if (DataTypes::getRank(inShape) == 2) {
253  int i0, i1;
254  int s0=inShape[0];
255  int s1=inShape[1];
256  for (i0=0; i0<s0; i0++) {
257  for (i1=0; i1<s1; i1++) {
258  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = (in[inOffset+DataTypes::getRelIndex(inShape,i0,i1)] - in[inOffset+DataTypes::getRelIndex(inShape,i1,i0)]) / 2.0;
259  }
260  }
261  }
262  else if (DataTypes::getRank(inShape) == 4) {
263  int i0, i1, i2, i3;
264  int s0=inShape[0];
265  int s1=inShape[1];
266  int s2=inShape[2];
267  int s3=inShape[3];
268  for (i0=0; i0<s0; i0++) {
269  for (i1=0; i1<s1; i1++) {
270  for (i2=0; i2<s2; i2++) {
271  for (i3=0; i3<s3; i3++) {
272  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;
273  }
274  }
275  }
276  }
277  }
278  }
279 
292  inline
293  void
295  const DataTypes::ShapeType& inShape,
298  const DataTypes::ShapeType& evShape,
300  int axis_offset)
301  {
302  for (int j=0;j<DataTypes::noValues(evShape);++j)
303  {
304  ev[evOffset+j]=0;
305  }
306  if (DataTypes::getRank(inShape) == 2) {
307  int s0=inShape[0]; // Python wrapper limits to square matrix
308  int i;
309  for (i=0; i<s0; i++) {
310  ev[evOffset/*+DataTypes::getRelIndex(evShape)*/] += in[inOffset+DataTypes::getRelIndex(inShape,i,i)];
311  }
312  }
313  else if (DataTypes::getRank(inShape) == 3) {
314  if (axis_offset==0) {
315  int s0=inShape[0];
316  int s2=inShape[2];
317  int i0, i2;
318  for (i0=0; i0<s0; i0++) {
319  for (i2=0; i2<s2; i2++) {
320  ev[evOffset+DataTypes::getRelIndex(evShape,i2)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i0,i2)];
321  }
322  }
323  }
324  else if (axis_offset==1) {
325  int s0=inShape[0];
326  int s1=inShape[1];
327  int i0, i1;
328  for (i0=0; i0<s0; i0++) {
329  for (i1=0; i1<s1; i1++) {
330  ev[evOffset+DataTypes::getRelIndex(evShape,i0)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i1)];
331  }
332  }
333  }
334  }
335  else if (DataTypes::getRank(inShape) == 4) {
336  if (axis_offset==0) {
337  int s0=inShape[0];
338  int s2=inShape[2];
339  int s3=inShape[3];
340  int i0, i2, i3;
341  for (i0=0; i0<s0; i0++) {
342  for (i2=0; i2<s2; i2++) {
343  for (i3=0; i3<s3; i3++) {
344  ev[evOffset+DataTypes::getRelIndex(evShape,i2,i3)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i0,i2,i3)];
345  }
346  }
347  }
348  }
349  else if (axis_offset==1) {
350  int s0=inShape[0];
351  int s1=inShape[1];
352  int s3=inShape[3];
353  int i0, i1, i3;
354  for (i0=0; i0<s0; i0++) {
355  for (i1=0; i1<s1; i1++) {
356  for (i3=0; i3<s3; i3++) {
357  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i3)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i1,i3)];
358  }
359  }
360  }
361  }
362  else if (axis_offset==2) {
363  int s0=inShape[0];
364  int s1=inShape[1];
365  int s2=inShape[2];
366  int i0, i1, i2;
367  for (i0=0; i0<s0; i0++) {
368  for (i1=0; i1<s1; i1++) {
369  for (i2=0; i2<s2; i2++) {
370  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i2,i2)];
371  }
372  }
373  }
374  }
375  }
376  }
377 
391  inline
392  void
394  const DataTypes::ShapeType& inShape,
397  const DataTypes::ShapeType& evShape,
399  int axis_offset)
400  {
401  int inRank=DataTypes::getRank(inShape);
402  if ( inRank== 4) {
403  int s0=evShape[0];
404  int s1=evShape[1];
405  int s2=evShape[2];
406  int s3=evShape[3];
407  int i0, i1, i2, i3;
408  if (axis_offset==1) {
409  for (i0=0; i0<s0; i0++) {
410  for (i1=0; i1<s1; i1++) {
411  for (i2=0; i2<s2; i2++) {
412  for (i3=0; i3<s3; i3++) {
413  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i3,i0,i1,i2)];
414  }
415  }
416  }
417  }
418  }
419  else if (axis_offset==2) {
420  for (i0=0; i0<s0; i0++) {
421  for (i1=0; i1<s1; i1++) {
422  for (i2=0; i2<s2; i2++) {
423  for (i3=0; i3<s3; i3++) {
424  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i3,i0,i1)];
425  }
426  }
427  }
428  }
429  }
430  else if (axis_offset==3) {
431  for (i0=0; i0<s0; i0++) {
432  for (i1=0; i1<s1; i1++) {
433  for (i2=0; i2<s2; i2++) {
434  for (i3=0; i3<s3; i3++) {
435  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i2,i3,i0)];
436  }
437  }
438  }
439  }
440  }
441  else {
442  for (i0=0; i0<s0; i0++) {
443  for (i1=0; i1<s1; i1++) {
444  for (i2=0; i2<s2; i2++) {
445  for (i3=0; i3<s3; i3++) {
446  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i2,i3)];
447  }
448  }
449  }
450  }
451  }
452  }
453  else if (inRank == 3) {
454  int s0=evShape[0];
455  int s1=evShape[1];
456  int s2=evShape[2];
457  int i0, i1, i2;
458  if (axis_offset==1) {
459  for (i0=0; i0<s0; i0++) {
460  for (i1=0; i1<s1; i1++) {
461  for (i2=0; i2<s2; i2++) {
462  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i0,i1)];
463  }
464  }
465  }
466  }
467  else if (axis_offset==2) {
468  for (i0=0; i0<s0; i0++) {
469  for (i1=0; i1<s1; i1++) {
470  for (i2=0; i2<s2; i2++) {
471  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i2,i0)];
472  }
473  }
474  }
475  }
476  else {
477  // Copy the matrix unchanged
478  for (i0=0; i0<s0; i0++) {
479  for (i1=0; i1<s1; i1++) {
480  for (i2=0; i2<s2; i2++) {
481  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i2)];
482  }
483  }
484  }
485  }
486  }
487  else if (inRank == 2) {
488  int s0=evShape[0];
489  int s1=evShape[1];
490  int i0, i1;
491  if (axis_offset==1) {
492  for (i0=0; i0<s0; i0++) {
493  for (i1=0; i1<s1; i1++) {
494  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0)];
495  }
496  }
497  }
498  else {
499  for (i0=0; i0<s0; i0++) {
500  for (i1=0; i1<s1; i1++) {
501  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1)];
502  }
503  }
504  }
505  }
506  else if (inRank == 1) {
507  int s0=evShape[0];
508  int i0;
509  for (i0=0; i0<s0; i0++) {
510  ev[evOffset+DataTypes::getRelIndex(evShape,i0)] = in[inOffset+DataTypes::getRelIndex(inShape,i0)];
511  }
512  }
513  else if (inRank == 0) {
514  ev[evOffset/*+DataTypes::getRelIndex(evShape,)*/] = in[inOffset/*+DataTypes::getRelIndex(inShape,)*/];
515  }
516  else {
517  throw DataException("Error - DataArrayView::transpose can only be calculated for rank 0, 1, 2, 3 or 4 objects.");
518  }
519  }
520 
535  inline
536  void
538  const DataTypes::ShapeType& inShape,
541  const DataTypes::ShapeType& evShape,
543  int axis0,
544  int axis1)
545  {
546  int inRank=DataTypes::getRank(inShape);
547  if (inRank == 4) {
548  int s0=evShape[0];
549  int s1=evShape[1];
550  int s2=evShape[2];
551  int s3=evShape[3];
552  int i0, i1, i2, i3;
553  if (axis0==0) {
554  if (axis1==1) {
555  for (i0=0; i0<s0; i0++) {
556  for (i1=0; i1<s1; i1++) {
557  for (i2=0; i2<s2; i2++) {
558  for (i3=0; i3<s3; i3++) {
559  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0,i2,i3)];
560  }
561  }
562  }
563  }
564  } else if (axis1==2) {
565  for (i0=0; i0<s0; i0++) {
566  for (i1=0; i1<s1; i1++) {
567  for (i2=0; i2<s2; i2++) {
568  for (i3=0; i3<s3; i3++) {
569  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i1,i0,i3)];
570  }
571  }
572  }
573  }
574 
575  } else if (axis1==3) {
576  for (i0=0; i0<s0; i0++) {
577  for (i1=0; i1<s1; i1++) {
578  for (i2=0; i2<s2; i2++) {
579  for (i3=0; i3<s3; i3++) {
580  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i3,i1,i2,i0)];
581  }
582  }
583  }
584  }
585  }
586  } else if (axis0==1) {
587  if (axis1==2) {
588  for (i0=0; i0<s0; i0++) {
589  for (i1=0; i1<s1; i1++) {
590  for (i2=0; i2<s2; i2++) {
591  for (i3=0; i3<s3; i3++) {
592  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i2,i1,i3)];
593  }
594  }
595  }
596  }
597  } else if (axis1==3) {
598  for (i0=0; i0<s0; i0++) {
599  for (i1=0; i1<s1; i1++) {
600  for (i2=0; i2<s2; i2++) {
601  for (i3=0; i3<s3; i3++) {
602  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i3,i2,i1)];
603  }
604  }
605  }
606  }
607  }
608  } else if (axis0==2) {
609  if (axis1==3) {
610  for (i0=0; i0<s0; i0++) {
611  for (i1=0; i1<s1; i1++) {
612  for (i2=0; i2<s2; i2++) {
613  for (i3=0; i3<s3; i3++) {
614  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i3,i2)];
615  }
616  }
617  }
618  }
619  }
620  }
621 
622  } else if ( inRank == 3) {
623  int s0=evShape[0];
624  int s1=evShape[1];
625  int s2=evShape[2];
626  int i0, i1, i2;
627  if (axis0==0) {
628  if (axis1==1) {
629  for (i0=0; i0<s0; i0++) {
630  for (i1=0; i1<s1; i1++) {
631  for (i2=0; i2<s2; i2++) {
632  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0,i2)];
633  }
634  }
635  }
636  } else if (axis1==2) {
637  for (i0=0; i0<s0; i0++) {
638  for (i1=0; i1<s1; i1++) {
639  for (i2=0; i2<s2; i2++) {
640  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i1,i0)];
641  }
642  }
643  }
644  }
645  } else if (axis0==1) {
646  if (axis1==2) {
647  for (i0=0; i0<s0; i0++) {
648  for (i1=0; i1<s1; i1++) {
649  for (i2=0; i2<s2; i2++) {
650  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i2,i1)];
651  }
652  }
653  }
654  }
655  }
656  } else if ( inRank == 2) {
657  int s0=evShape[0];
658  int s1=evShape[1];
659  int i0, i1;
660  if (axis0==0) {
661  if (axis1==1) {
662  for (i0=0; i0<s0; i0++) {
663  for (i1=0; i1<s1; i1++) {
664  ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0)];
665  }
666  }
667  }
668  }
669  } else {
670  throw DataException("Error - DataArrayView::swapaxes can only be calculated for rank 2, 3 or 4 objects.");
671  }
672  }
673 
686  inline
687  void
689  const DataTypes::ShapeType& inShape,
692  const DataTypes::ShapeType& evShape,
694  {
695  double in00,in10,in20,in01,in11,in21,in02,in12,in22;
696  double ev0,ev1,ev2;
697  int s=inShape[0];
698  if (s==1) {
699  in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
700  eigenvalues1(in00,&ev0);
701  ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
702 
703  } else if (s==2) {
704  in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
705  in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
706  in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
707  in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
708  eigenvalues2(in00,(in01+in10)/2.,in11,&ev0,&ev1);
709  ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
710  ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
711 
712  } else if (s==3) {
713  in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
714  in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
715  in20=in[inOffset+DataTypes::getRelIndex(inShape,2,0)];
716  in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
717  in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
718  in21=in[inOffset+DataTypes::getRelIndex(inShape,2,1)];
719  in02=in[inOffset+DataTypes::getRelIndex(inShape,0,2)];
720  in12=in[inOffset+DataTypes::getRelIndex(inShape,1,2)];
721  in22=in[inOffset+DataTypes::getRelIndex(inShape,2,2)];
722  eigenvalues3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
723  &ev0,&ev1,&ev2);
724  ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
725  ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
726  ev[evOffset+DataTypes::getRelIndex(evShape,2)]=ev2;
727 
728  }
729  }
730 
747  inline
748  void
751  DataTypes::ValueType& ev, const DataTypes::ShapeType& evShape,
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 inline
814 bool
816  const DataTypes::ShapeType& shape,
818 {
819  return (data.size() >= (offset+DataTypes::noValues(shape)));
820 }
821 
822 template <class UnaryFunction>
823 inline
824 void
827  UnaryFunction operation)
828 {
829  EsysAssert((data.size()>0)&&checkOffset(data,shape,offset),
830  "Error - Couldn't perform unaryOp due to insufficient storage.");
832  for (DataTypes::ValueType::size_type i=0;i<nVals;i++) {
833  data[offset+i]=operation(data[offset+i]);
834  }
835 }
836 
837 
838 template <class BinaryFunction>
839 inline
840 void
842  const DataTypes::ShapeType& leftShape,
844  const DataTypes::ValueType& right,
845  const DataTypes::ShapeType& rightShape,
847  BinaryFunction operation)
848 {
849  EsysAssert(leftShape==rightShape,
850  "Error - Couldn't perform binaryOp due to shape mismatch,");
851  EsysAssert(((left.size()>0)&&checkOffset(left,leftShape, leftOffset)),
852  "Error - Couldn't perform binaryOp due to insufficient storage in left object.");
853  EsysAssert(((right.size()>0)&&checkOffset(right,rightShape,rightOffset)),
854  "Error - Couldn't perform binaryOp due to insufficient storage in right object.");
855  for (DataTypes::ValueType::size_type i=0;i<DataTypes::noValues(leftShape);i++) {
856  left[leftOffset+i]=operation(left[leftOffset+i],right[rightOffset+i]);
857  }
858 }
859 
860 template <class BinaryFunction>
861 inline
862 void
864  const DataTypes::ShapeType& leftShape,
866  double right,
867  BinaryFunction operation)
868 {
869  EsysAssert(((left.size()>0)&&checkOffset(left,leftShape,offset)),
870  "Error - Couldn't perform binaryOp due to insufficient storage in left object.");
871  for (DataTypes::ValueType::size_type i=0;i<DataTypes::noValues(leftShape);i++) {
872  left[offset+i]=operation(left[offset+i],right);
873  }
874 }
875 
876 template <class BinaryFunction>
877 inline
878 double
880  const DataTypes::ShapeType& leftShape,
882  BinaryFunction operation,
883  double initial_value)
884 {
885  EsysAssert(((left.size()>0)&&checkOffset(left,leftShape,offset)),
886  "Error - Couldn't perform reductionOp due to insufficient storage.");
887  double current_value=initial_value;
888  for (DataTypes::ValueType::size_type i=0;i<DataTypes::noValues(leftShape);i++) {
889  current_value=operation(current_value,left[offset+i]);
890  }
891  return current_value;
892 }
893 
910 int
912  const DataTypes::ShapeType& inShape,
915  const DataTypes::ShapeType& outShape,
917  int count,
918  LapackInverseHelper& helper);
919 
927 void
928 matrixInverseError(int err);
929 
934 inline
935 bool
937 {
938  for (size_t z=inOffset;z<inOffset+count;++z)
939  {
940  if (nancheck(in[z]))
941  {
942  return true;
943  }
944  }
945  return false;
946 }
947 
948 } // end namespace DataMath
949 } // end namespace escript
950 #endif
951