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