ESYS13: DataMaths.h Source File

Go to the documentation of this file.
00001 
00002 /*******************************************************
00003 *
00004 * Copyright (c) 2003-2010 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