00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015 #ifndef escript_DataMaths_20080822_H
00016 #define escript_DataMaths_20080822_H
00017 #include "DataAbstract.h"
00018 #include "DataException.h"
00019 #include "LocalOps.h"
00020
00031 namespace escript
00032 {
00033 namespace DataMaths
00034 {
00035
00059 template <class UnaryFunction>
00060 void
00061 unaryOp(DataTypes::ValueType& data, const DataTypes::ShapeType& shape,
00062 DataTypes::ValueType::size_type offset,
00063 UnaryFunction operation);
00064
00079 template <class BinaryFunction>
00080 void
00081 binaryOp(DataTypes::ValueType& left,
00082 const DataTypes::ShapeType& leftShape,
00083 DataTypes::ValueType::size_type leftOffset,
00084 const DataTypes::ValueType& right,
00085 const DataTypes::ShapeType& rightShape,
00086 DataTypes::ValueType::size_type rightOffset,
00087 BinaryFunction operation);
00088
00105 template <class BinaryFunction>
00106 void
00107 binaryOp(DataTypes::ValueType& left,
00108 const DataTypes::ShapeType& shape,
00109 DataTypes::ValueType::size_type offset,
00110 double right,
00111 BinaryFunction operation);
00112
00129 template <class BinaryFunction>
00130 double
00131 reductionOp(const DataTypes::ValueType& left,
00132 const DataTypes::ShapeType& shape,
00133 DataTypes::ValueType::size_type offset,
00134 BinaryFunction operation,
00135 double initial_value);
00136
00151 ESCRIPT_DLL_API
00152 void
00153 matMult(const DataTypes::ValueType& left,
00154 const DataTypes::ShapeType& leftShape,
00155 DataTypes::ValueType::size_type leftOffset,
00156 const DataTypes::ValueType& right,
00157 const DataTypes::ShapeType& rightShape,
00158 DataTypes::ValueType::size_type rightOffset,
00159 DataTypes::ValueType& result,
00160 const DataTypes::ShapeType& resultShape);
00161
00162
00163
00164
00165
00174 ESCRIPT_DLL_API
00175 DataTypes::ShapeType
00176 determineResultShape(const DataTypes::ShapeType& left,
00177 const DataTypes::ShapeType& right);
00178
00190 ESCRIPT_DLL_API
00191 inline
00192 void
00193 symmetric(const DataTypes::ValueType& in,
00194 const DataTypes::ShapeType& inShape,
00195 DataTypes::ValueType::size_type inOffset,
00196 DataTypes::ValueType& ev,
00197 const DataTypes::ShapeType& evShape,
00198 DataTypes::ValueType::size_type evOffset)
00199 {
00200 if (DataTypes::getRank(inShape) == 2) {
00201 int i0, i1;
00202 int s0=inShape[0];
00203 int s1=inShape[1];
00204 for (i0=0; i0<s0; i0++) {
00205 for (i1=0; i1<s1; i1++) {
00206 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = (in[inOffset+DataTypes::getRelIndex(inShape,i0,i1)] + in[inOffset+DataTypes::getRelIndex(inShape,i1,i0)]) / 2.0;
00207 }
00208 }
00209 }
00210 else if (DataTypes::getRank(inShape) == 4) {
00211 int i0, i1, i2, i3;
00212 int s0=inShape[0];
00213 int s1=inShape[1];
00214 int s2=inShape[2];
00215 int s3=inShape[3];
00216 for (i0=0; i0<s0; i0++) {
00217 for (i1=0; i1<s1; i1++) {
00218 for (i2=0; i2<s2; i2++) {
00219 for (i3=0; i3<s3; i3++) {
00220 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;
00221 }
00222 }
00223 }
00224 }
00225 }
00226 }
00227
00239 ESCRIPT_DLL_API
00240 inline
00241 void
00242 nonsymmetric(const DataTypes::ValueType& in,
00243 const DataTypes::ShapeType& inShape,
00244 DataTypes::ValueType::size_type inOffset,
00245 DataTypes::ValueType& ev,
00246 const DataTypes::ShapeType& evShape,
00247 DataTypes::ValueType::size_type evOffset)
00248 {
00249 if (DataTypes::getRank(inShape) == 2) {
00250 int i0, i1;
00251 int s0=inShape[0];
00252 int s1=inShape[1];
00253 for (i0=0; i0<s0; i0++) {
00254 for (i1=0; i1<s1; i1++) {
00255 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = (in[inOffset+DataTypes::getRelIndex(inShape,i0,i1)] - in[inOffset+DataTypes::getRelIndex(inShape,i1,i0)]) / 2.0;
00256 }
00257 }
00258 }
00259 else if (DataTypes::getRank(inShape) == 4) {
00260 int i0, i1, i2, i3;
00261 int s0=inShape[0];
00262 int s1=inShape[1];
00263 int s2=inShape[2];
00264 int s3=inShape[3];
00265 for (i0=0; i0<s0; i0++) {
00266 for (i1=0; i1<s1; i1++) {
00267 for (i2=0; i2<s2; i2++) {
00268 for (i3=0; i3<s3; i3++) {
00269 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;
00270 }
00271 }
00272 }
00273 }
00274 }
00275 }
00276
00289 inline
00290 void
00291 trace(const DataTypes::ValueType& in,
00292 const DataTypes::ShapeType& inShape,
00293 DataTypes::ValueType::size_type inOffset,
00294 DataTypes::ValueType& ev,
00295 const DataTypes::ShapeType& evShape,
00296 DataTypes::ValueType::size_type evOffset,
00297 int axis_offset)
00298 {
00299 for (int j=0;j<DataTypes::noValues(evShape);++j)
00300 {
00301 ev[evOffset+j]=0;
00302 }
00303 if (DataTypes::getRank(inShape) == 2) {
00304 int s0=inShape[0];
00305 int i;
00306 for (i=0; i<s0; i++) {
00307 ev[evOffset] += in[inOffset+DataTypes::getRelIndex(inShape,i,i)];
00308 }
00309 }
00310 else if (DataTypes::getRank(inShape) == 3) {
00311 if (axis_offset==0) {
00312 int s0=inShape[0];
00313 int s2=inShape[2];
00314 int i0, i2;
00315 for (i0=0; i0<s0; i0++) {
00316 for (i2=0; i2<s2; i2++) {
00317 ev[evOffset+DataTypes::getRelIndex(evShape,i2)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i0,i2)];
00318 }
00319 }
00320 }
00321 else if (axis_offset==1) {
00322 int s0=inShape[0];
00323 int s1=inShape[1];
00324 int i0, i1;
00325 for (i0=0; i0<s0; i0++) {
00326 for (i1=0; i1<s1; i1++) {
00327 ev[evOffset+DataTypes::getRelIndex(evShape,i0)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i1)];
00328 }
00329 }
00330 }
00331 }
00332 else if (DataTypes::getRank(inShape) == 4) {
00333 if (axis_offset==0) {
00334 int s0=inShape[0];
00335 int s2=inShape[2];
00336 int s3=inShape[3];
00337 int i0, i2, i3;
00338 for (i0=0; i0<s0; i0++) {
00339 for (i2=0; i2<s2; i2++) {
00340 for (i3=0; i3<s3; i3++) {
00341 ev[evOffset+DataTypes::getRelIndex(evShape,i2,i3)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i0,i2,i3)];
00342 }
00343 }
00344 }
00345 }
00346 else if (axis_offset==1) {
00347 int s0=inShape[0];
00348 int s1=inShape[1];
00349 int s3=inShape[3];
00350 int i0, i1, i3;
00351 for (i0=0; i0<s0; i0++) {
00352 for (i1=0; i1<s1; i1++) {
00353 for (i3=0; i3<s3; i3++) {
00354 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i3)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i1,i3)];
00355 }
00356 }
00357 }
00358 }
00359 else if (axis_offset==2) {
00360 int s0=inShape[0];
00361 int s1=inShape[1];
00362 int s2=inShape[2];
00363 int i0, i1, i2;
00364 for (i0=0; i0<s0; i0++) {
00365 for (i1=0; i1<s1; i1++) {
00366 for (i2=0; i2<s2; i2++) {
00367 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i2,i2)];
00368 }
00369 }
00370 }
00371 }
00372 }
00373 }
00374
00387 ESCRIPT_DLL_API
00388 inline
00389 void
00390 transpose(const DataTypes::ValueType& in,
00391 const DataTypes::ShapeType& inShape,
00392 DataTypes::ValueType::size_type inOffset,
00393 DataTypes::ValueType& ev,
00394 const DataTypes::ShapeType& evShape,
00395 DataTypes::ValueType::size_type evOffset,
00396 int axis_offset)
00397 {
00398 int inRank=DataTypes::getRank(inShape);
00399 if ( inRank== 4) {
00400 int s0=evShape[0];
00401 int s1=evShape[1];
00402 int s2=evShape[2];
00403 int s3=evShape[3];
00404 int i0, i1, i2, i3;
00405 if (axis_offset==1) {
00406 for (i0=0; i0<s0; i0++) {
00407 for (i1=0; i1<s1; i1++) {
00408 for (i2=0; i2<s2; i2++) {
00409 for (i3=0; i3<s3; i3++) {
00410 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i3,i0,i1,i2)];
00411 }
00412 }
00413 }
00414 }
00415 }
00416 else if (axis_offset==2) {
00417 for (i0=0; i0<s0; i0++) {
00418 for (i1=0; i1<s1; i1++) {
00419 for (i2=0; i2<s2; i2++) {
00420 for (i3=0; i3<s3; i3++) {
00421 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i3,i0,i1)];
00422 }
00423 }
00424 }
00425 }
00426 }
00427 else if (axis_offset==3) {
00428 for (i0=0; i0<s0; i0++) {
00429 for (i1=0; i1<s1; i1++) {
00430 for (i2=0; i2<s2; i2++) {
00431 for (i3=0; i3<s3; i3++) {
00432 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i2,i3,i0)];
00433 }
00434 }
00435 }
00436 }
00437 }
00438 else {
00439 for (i0=0; i0<s0; i0++) {
00440 for (i1=0; i1<s1; i1++) {
00441 for (i2=0; i2<s2; i2++) {
00442 for (i3=0; i3<s3; i3++) {
00443 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i2,i3)];
00444 }
00445 }
00446 }
00447 }
00448 }
00449 }
00450 else if (inRank == 3) {
00451 int s0=evShape[0];
00452 int s1=evShape[1];
00453 int s2=evShape[2];
00454 int i0, i1, i2;
00455 if (axis_offset==1) {
00456 for (i0=0; i0<s0; i0++) {
00457 for (i1=0; i1<s1; i1++) {
00458 for (i2=0; i2<s2; i2++) {
00459 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i0,i1)];
00460 }
00461 }
00462 }
00463 }
00464 else if (axis_offset==2) {
00465 for (i0=0; i0<s0; i0++) {
00466 for (i1=0; i1<s1; i1++) {
00467 for (i2=0; i2<s2; i2++) {
00468 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i2,i0)];
00469 }
00470 }
00471 }
00472 }
00473 else {
00474
00475 for (i0=0; i0<s0; i0++) {
00476 for (i1=0; i1<s1; i1++) {
00477 for (i2=0; i2<s2; i2++) {
00478 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i2)];
00479 }
00480 }
00481 }
00482 }
00483 }
00484 else if (inRank == 2) {
00485 int s0=evShape[0];
00486 int s1=evShape[1];
00487 int i0, i1;
00488 if (axis_offset==1) {
00489 for (i0=0; i0<s0; i0++) {
00490 for (i1=0; i1<s1; i1++) {
00491 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0)];
00492 }
00493 }
00494 }
00495 else {
00496 for (i0=0; i0<s0; i0++) {
00497 for (i1=0; i1<s1; i1++) {
00498 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1)];
00499 }
00500 }
00501 }
00502 }
00503 else if (inRank == 1) {
00504 int s0=evShape[0];
00505 int i0;
00506 for (i0=0; i0<s0; i0++) {
00507 ev[evOffset+DataTypes::getRelIndex(evShape,i0)] = in[inOffset+DataTypes::getRelIndex(inShape,i0)];
00508 }
00509 }
00510 else if (inRank == 0) {
00511 ev[evOffset] = in[inOffset];
00512 }
00513 else {
00514 throw DataException("Error - DataArrayView::transpose can only be calculated for rank 0, 1, 2, 3 or 4 objects.");
00515 }
00516 }
00517
00531 ESCRIPT_DLL_API
00532 inline
00533 void
00534 swapaxes(const DataTypes::ValueType& in,
00535 const DataTypes::ShapeType& inShape,
00536 DataTypes::ValueType::size_type inOffset,
00537 DataTypes::ValueType& ev,
00538 const DataTypes::ShapeType& evShape,
00539 DataTypes::ValueType::size_type evOffset,
00540 int axis0,
00541 int axis1)
00542 {
00543 int inRank=DataTypes::getRank(inShape);
00544 if (inRank == 4) {
00545 int s0=evShape[0];
00546 int s1=evShape[1];
00547 int s2=evShape[2];
00548 int s3=evShape[3];
00549 int i0, i1, i2, i3;
00550 if (axis0==0) {
00551 if (axis1==1) {
00552 for (i0=0; i0<s0; i0++) {
00553 for (i1=0; i1<s1; i1++) {
00554 for (i2=0; i2<s2; i2++) {
00555 for (i3=0; i3<s3; i3++) {
00556 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0,i2,i3)];
00557 }
00558 }
00559 }
00560 }
00561 } else if (axis1==2) {
00562 for (i0=0; i0<s0; i0++) {
00563 for (i1=0; i1<s1; i1++) {
00564 for (i2=0; i2<s2; i2++) {
00565 for (i3=0; i3<s3; i3++) {
00566 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i1,i0,i3)];
00567 }
00568 }
00569 }
00570 }
00571
00572 } else if (axis1==3) {
00573 for (i0=0; i0<s0; i0++) {
00574 for (i1=0; i1<s1; i1++) {
00575 for (i2=0; i2<s2; i2++) {
00576 for (i3=0; i3<s3; i3++) {
00577 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i3,i1,i2,i0)];
00578 }
00579 }
00580 }
00581 }
00582 }
00583 } else if (axis0==1) {
00584 if (axis1==2) {
00585 for (i0=0; i0<s0; i0++) {
00586 for (i1=0; i1<s1; i1++) {
00587 for (i2=0; i2<s2; i2++) {
00588 for (i3=0; i3<s3; i3++) {
00589 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i2,i1,i3)];
00590 }
00591 }
00592 }
00593 }
00594 } else if (axis1==3) {
00595 for (i0=0; i0<s0; i0++) {
00596 for (i1=0; i1<s1; i1++) {
00597 for (i2=0; i2<s2; i2++) {
00598 for (i3=0; i3<s3; i3++) {
00599 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i3,i2,i1)];
00600 }
00601 }
00602 }
00603 }
00604 }
00605 } else if (axis0==2) {
00606 if (axis1==3) {
00607 for (i0=0; i0<s0; i0++) {
00608 for (i1=0; i1<s1; i1++) {
00609 for (i2=0; i2<s2; i2++) {
00610 for (i3=0; i3<s3; i3++) {
00611 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i3,i2)];
00612 }
00613 }
00614 }
00615 }
00616 }
00617 }
00618
00619 } else if ( inRank == 3) {
00620 int s0=evShape[0];
00621 int s1=evShape[1];
00622 int s2=evShape[2];
00623 int i0, i1, i2;
00624 if (axis0==0) {
00625 if (axis1==1) {
00626 for (i0=0; i0<s0; i0++) {
00627 for (i1=0; i1<s1; i1++) {
00628 for (i2=0; i2<s2; i2++) {
00629 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0,i2)];
00630 }
00631 }
00632 }
00633 } else if (axis1==2) {
00634 for (i0=0; i0<s0; i0++) {
00635 for (i1=0; i1<s1; i1++) {
00636 for (i2=0; i2<s2; i2++) {
00637 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i1,i0)];
00638 }
00639 }
00640 }
00641 }
00642 } else if (axis0==1) {
00643 if (axis1==2) {
00644 for (i0=0; i0<s0; i0++) {
00645 for (i1=0; i1<s1; i1++) {
00646 for (i2=0; i2<s2; i2++) {
00647 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i2,i1)];
00648 }
00649 }
00650 }
00651 }
00652 }
00653 } else if ( inRank == 2) {
00654 int s0=evShape[0];
00655 int s1=evShape[1];
00656 int i0, i1;
00657 if (axis0==0) {
00658 if (axis1==1) {
00659 for (i0=0; i0<s0; i0++) {
00660 for (i1=0; i1<s1; i1++) {
00661 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0)];
00662 }
00663 }
00664 }
00665 }
00666 } else {
00667 throw DataException("Error - DataArrayView::swapaxes can only be calculated for rank 2, 3 or 4 objects.");
00668 }
00669 }
00670
00682 ESCRIPT_DLL_API
00683 inline
00684 void
00685 eigenvalues(const DataTypes::ValueType& in,
00686 const DataTypes::ShapeType& inShape,
00687 DataTypes::ValueType::size_type inOffset,
00688 DataTypes::ValueType& ev,
00689 const DataTypes::ShapeType& evShape,
00690 DataTypes::ValueType::size_type evOffset)
00691 {
00692 double in00,in10,in20,in01,in11,in21,in02,in12,in22;
00693 double ev0,ev1,ev2;
00694 int s=inShape[0];
00695 if (s==1) {
00696 in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
00697 eigenvalues1(in00,&ev0);
00698 ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
00699
00700 } else if (s==2) {
00701 in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
00702 in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
00703 in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
00704 in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
00705 eigenvalues2(in00,(in01+in10)/2.,in11,&ev0,&ev1);
00706 ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
00707 ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
00708
00709 } else if (s==3) {
00710 in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
00711 in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
00712 in20=in[inOffset+DataTypes::getRelIndex(inShape,2,0)];
00713 in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
00714 in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
00715 in21=in[inOffset+DataTypes::getRelIndex(inShape,2,1)];
00716 in02=in[inOffset+DataTypes::getRelIndex(inShape,0,2)];
00717 in12=in[inOffset+DataTypes::getRelIndex(inShape,1,2)];
00718 in22=in[inOffset+DataTypes::getRelIndex(inShape,2,2)];
00719 eigenvalues3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
00720 &ev0,&ev1,&ev2);
00721 ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
00722 ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
00723 ev[evOffset+DataTypes::getRelIndex(evShape,2)]=ev2;
00724
00725 }
00726 }
00727
00743 ESCRIPT_DLL_API
00744 inline
00745 void
00746 eigenvalues_and_eigenvectors(const DataTypes::ValueType& in, const DataTypes::ShapeType& inShape,
00747 DataTypes::ValueType::size_type inOffset,
00748 DataTypes::ValueType& ev, const DataTypes::ShapeType& evShape,
00749 DataTypes::ValueType::size_type evOffset,
00750 DataTypes::ValueType& V, const DataTypes::ShapeType& VShape,
00751 DataTypes::ValueType::size_type VOffset,
00752 const double tol=1.e-13)
00753 {
00754 double in00,in10,in20,in01,in11,in21,in02,in12,in22;
00755 double V00,V10,V20,V01,V11,V21,V02,V12,V22;
00756 double ev0,ev1,ev2;
00757 int s=inShape[0];
00758 if (s==1) {
00759 in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
00760 eigenvalues_and_eigenvectors1(in00,&ev0,&V00,tol);
00761 ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
00762 V[inOffset+DataTypes::getRelIndex(VShape,0,0)]=V00;
00763 } else if (s==2) {
00764 in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
00765 in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
00766 in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
00767 in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
00768 eigenvalues_and_eigenvectors2(in00,(in01+in10)/2.,in11,
00769 &ev0,&ev1,&V00,&V10,&V01,&V11,tol);
00770 ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
00771 ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
00772 V[inOffset+DataTypes::getRelIndex(VShape,0,0)]=V00;
00773 V[inOffset+DataTypes::getRelIndex(VShape,1,0)]=V10;
00774 V[inOffset+DataTypes::getRelIndex(VShape,0,1)]=V01;
00775 V[inOffset+DataTypes::getRelIndex(VShape,1,1)]=V11;
00776 } else if (s==3) {
00777 in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
00778 in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
00779 in20=in[inOffset+DataTypes::getRelIndex(inShape,2,0)];
00780 in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
00781 in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
00782 in21=in[inOffset+DataTypes::getRelIndex(inShape,2,1)];
00783 in02=in[inOffset+DataTypes::getRelIndex(inShape,0,2)];
00784 in12=in[inOffset+DataTypes::getRelIndex(inShape,1,2)];
00785 in22=in[inOffset+DataTypes::getRelIndex(inShape,2,2)];
00786 eigenvalues_and_eigenvectors3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
00787 &ev0,&ev1,&ev2,
00788 &V00,&V10,&V20,&V01,&V11,&V21,&V02,&V12,&V22,tol);
00789 ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
00790 ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
00791 ev[evOffset+DataTypes::getRelIndex(evShape,2)]=ev2;
00792 V[inOffset+DataTypes::getRelIndex(VShape,0,0)]=V00;
00793 V[inOffset+DataTypes::getRelIndex(VShape,1,0)]=V10;
00794 V[inOffset+DataTypes::getRelIndex(VShape,2,0)]=V20;
00795 V[inOffset+DataTypes::getRelIndex(VShape,0,1)]=V01;
00796 V[inOffset+DataTypes::getRelIndex(VShape,1,1)]=V11;
00797 V[inOffset+DataTypes::getRelIndex(VShape,2,1)]=V21;
00798 V[inOffset+DataTypes::getRelIndex(VShape,0,2)]=V02;
00799 V[inOffset+DataTypes::getRelIndex(VShape,1,2)]=V12;
00800 V[inOffset+DataTypes::getRelIndex(VShape,2,2)]=V22;
00801
00802 }
00803 }
00804
00805
00810 inline
00811 bool
00812 checkOffset(const DataTypes::ValueType& data,
00813 const DataTypes::ShapeType& shape,
00814 DataTypes::ValueType::size_type offset)
00815 {
00816 return (data.size() >= (offset+DataTypes::noValues(shape)));
00817 }
00818
00819 template <class UnaryFunction>
00820 inline
00821 void
00822 unaryOp(DataTypes::ValueType& data, const DataTypes::ShapeType& shape,
00823 DataTypes::ValueType::size_type offset,
00824 UnaryFunction operation)
00825 {
00826 EsysAssert((data.size()>0)&&checkOffset(data,shape,offset),
00827 "Error - Couldn't perform unaryOp due to insufficient storage.");
00828 DataTypes::ValueType::size_type nVals=DataTypes::noValues(shape);
00829 for (DataTypes::ValueType::size_type i=0;i<nVals;i++) {
00830 data[offset+i]=operation(data[offset+i]);
00831 }
00832 }
00833
00834
00835 template <class BinaryFunction>
00836 inline
00837 void
00838 binaryOp(DataTypes::ValueType& left,
00839 const DataTypes::ShapeType& leftShape,
00840 DataTypes::ValueType::size_type leftOffset,
00841 const DataTypes::ValueType& right,
00842 const DataTypes::ShapeType& rightShape,
00843 DataTypes::ValueType::size_type rightOffset,
00844 BinaryFunction operation)
00845 {
00846 EsysAssert(leftShape==rightShape,
00847 "Error - Couldn't perform binaryOp due to shape mismatch,");
00848 EsysAssert(((left.size()>0)&&checkOffset(left,leftShape, leftOffset)),
00849 "Error - Couldn't perform binaryOp due to insufficient storage in left object.");
00850 EsysAssert(((right.size()>0)&&checkOffset(right,rightShape,rightOffset)),
00851 "Error - Couldn't perform binaryOp due to insufficient storage in right object.");
00852 for (DataTypes::ValueType::size_type i=0;i<DataTypes::noValues(leftShape);i++) {
00853 left[leftOffset+i]=operation(left[leftOffset+i],right[rightOffset+i]);
00854 }
00855 }
00856
00857 template <class BinaryFunction>
00858 inline
00859 void
00860 binaryOp(DataTypes::ValueType& left,
00861 const DataTypes::ShapeType& leftShape,
00862 DataTypes::ValueType::size_type offset,
00863 double right,
00864 BinaryFunction operation)
00865 {
00866 EsysAssert(((left.size()>0)&&checkOffset(left,leftShape,offset)),
00867 "Error - Couldn't perform binaryOp due to insufficient storage in left object.");
00868 for (DataTypes::ValueType::size_type i=0;i<DataTypes::noValues(leftShape);i++) {
00869 left[offset+i]=operation(left[offset+i],right);
00870 }
00871 }
00872
00873 template <class BinaryFunction>
00874 inline
00875 double
00876 reductionOp(const DataTypes::ValueType& left,
00877 const DataTypes::ShapeType& leftShape,
00878 DataTypes::ValueType::size_type offset,
00879 BinaryFunction operation,
00880 double initial_value)
00881 {
00882 EsysAssert(((left.size()>0)&&checkOffset(left,leftShape,offset)),
00883 "Error - Couldn't perform reductionOp due to insufficient storage.");
00884 double current_value=initial_value;
00885 for (DataTypes::ValueType::size_type i=0;i<DataTypes::noValues(leftShape);i++) {
00886 current_value=operation(current_value,left[offset+i]);
00887 }
00888 return current_value;
00889 }
00890
00891
00892
00893 }
00894 }
00895 #endif
00896