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