escript  Revision_
DataMaths.h
Go to the documentation of this file.
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