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