# Contents of /branches/trilinos_from_5897/escriptcore/src/DataMaths.h

Revision 5963 - (show annotations)
Mon Feb 22 06:59:27 2016 UTC (3 years, 4 months ago) by caltinay
File MIME type: text/plain
File size: 40454 byte(s)
```sync and fix.

```
 1 2 /***************************************************************************** 3 * 4 * Copyright (c) 2003-2016 by The University of Queensland 5 * http://www.uq.edu.au 6 * 7 * Primary Business: Queensland, Australia 8 * Licensed under the Open Software License version 3.0 9 * http://www.opensource.org/licenses/osl-3.0.php 10 * 11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC) 12 * Development 2012-2013 by School of Earth Sciences 13 * Development from 2014 by Centre for Geoscience Computing (GeoComp) 14 * 15 *****************************************************************************/ 16 17 18 #ifndef escript_DataMaths_20080822_H 19 #define escript_DataMaths_20080822_H 20 #include "DataAbstract.h" 21 #include "DataException.h" 22 #include "LocalOps.h" 23 #include "LapackInverseHelper.h" 24 25 /** 26 \file DataMaths.h 27 \brief Describes binary operations performed on DataVector. 28 29 30 For operations on DataAbstract see BinaryOp.h. 31 For operations on double* see LocalOps.h. 32 */ 33 34 35 namespace escript 36 { 37 namespace DataMaths 38 { 39 40 /** 41 \namespace escript::DataMaths 42 \brief Contains maths operations performed on data vectors. 43 44 In order to properly identify the datapoints, in most cases, the vector, shape and offset of the point must all be supplied. 45 Note that vector in this context refers to a data vector storing datapoints not a mathematical vector. (However, datapoints within the data vector could represent scalars, vectors, matricies, ...). 46 */ 47 48 49 /** 50 \brief 51 Perform the unary operation on the data point specified by the given 52 offset. Applies the specified operation to each value in the data 53 point. Operation must be a pointer to a function. 54 55 Called by escript::unaryOp. 56 57 \param data - vector containing the datapoint 58 \param shape - shape of the point 59 \param offset - offset of the point within data 60 \param operation - Input - 61 Operation to apply. Must be a pointer to a function. 62 */ 63 template 64 void 65 unaryOp(DataTypes::RealVectorType& data, const DataTypes::ShapeType& shape, 66 DataTypes::RealVectorType::size_type offset, 67 UnaryFunction operation); 68 69 /** 70 \brief 71 Perform the binary operation on the data points specified by the given 72 offsets in the "left" and "right" vectors. Applies the specified operation 73 to corresponding values in both data points. Operation must be a pointer 74 to a function. 75 76 Called by escript::binaryOp. 77 \param left,right - vectors containing the datapoints 78 \param leftShape,rightShape - shapes of datapoints in the vectors 79 \param leftOffset,rightOffset - beginnings of datapoints in the vectors 80 \param operation - Input - 81 Operation to apply. Must be a pointer to a function. 82 */ 83 template 84 void 85 binaryOp(DataTypes::RealVectorType& left, 86 const DataTypes::ShapeType& leftShape, 87 DataTypes::RealVectorType::size_type leftOffset, 88 const DataTypes::RealVectorType& right, 89 const DataTypes::ShapeType& rightShape, 90 DataTypes::RealVectorType::size_type rightOffset, 91 BinaryFunction operation); 92 93 /** 94 \brief 95 Perform the binary operation on the data point specified by the given 96 offset in the vector using the scalar value "right". Applies the specified 97 operation to values in the data point. Operation must be a pointer 98 to a function. 99 100 Called by escript::binaryOp. 101 102 \param left - vector containing the datapoints 103 \param shape - shape of datapoint in the vector 104 \param offset - beginning of datapoint in the vector 105 \param right - scalar value for the right hand side of the operation 106 \param operation - Input - 107 Operation to apply. Must be a pointer to a function. 108 */ 109 template 110 void 111 binaryOp(DataTypes::RealVectorType& left, 112 const DataTypes::ShapeType& shape, 113 DataTypes::RealVectorType::size_type offset, 114 double right, 115 BinaryFunction operation); 116 117 // ------------------------ 118 119 /** 120 \brief 121 Perform the binary operation on the data points specified by the given 122 offsets in the "left" and "right" vectors. Applies the specified operation 123 to corresponding values in both data points. Operation must be a pointer 124 to a function. 125 126 Called by escript::binaryOp. 127 \param left,right - vectors containing the datapoints 128 \param leftShape,rightShape - shapes of datapoints in the vectors 129 \param leftOffset,rightOffset - beginnings of datapoints in the vectors 130 \param operation - Input - 131 Operation to apply. Must be a pointer to a function. 132 */ 133 template 134 void 135 binaryOpVector(LVEC& left, 136 const DataTypes::ShapeType& leftShape, 137 typename LVEC::size_type leftOffset, 138 const RVEC& right, 139 const DataTypes::ShapeType& rightShape, 140 typename RVEC::size_type rightOffset, 141 escript::ESFunction operation); 142 143 /** 144 \brief 145 Perform the binary operation on the data point specified by the given 146 offset in the vector using the scalar value "right". Applies the specified 147 operation to values in the data point. Operation must be a pointer 148 to a function. 149 150 Called by escript::binaryOp. 151 152 \param left - vector containing the datapoints 153 \param shape - shape of datapoint in the vector 154 \param offset - beginning of datapoint in the vector 155 \param right - scalar value for the right hand side of the operation 156 \param operation - Input - 157 Operation to apply. Must be a pointer to a function. 158 */ 159 template 160 void 161 binaryOpVector(LVEC& left, 162 const DataTypes::ShapeType& shape, 163 typename LVEC::size_type offset, 164 SCALAR right, 165 escript::ESFunction operation); 166 167 // ------------------ 168 169 170 /** 171 \brief 172 Perform the given data point reduction operation on the data point 173 specified by the given offset into the view. Reduces all elements of 174 the data point using the given operation, returning the result as a 175 scalar. Operation must be a pointer to a function. 176 177 Called by escript::algorithm. 178 179 \param left - vector containing the datapoint 180 \param shape - shape of datapoints in the vector 181 \param offset - beginning of datapoint in the vector 182 \param operation - Input - 183 Operation to apply. Must be a pointer to a function. 184 \param initial_value 185 */ 186 template 187 double 188 reductionOp(const DataTypes::RealVectorType& left, 189 const DataTypes::ShapeType& shape, 190 DataTypes::RealVectorType::size_type offset, 191 BinaryFunction operation, 192 double initial_value); 193 194 /** 195 \brief 196 Perform a matrix multiply of the given views. 197 198 NB: Only multiplies together the two given datapoints, 199 would need to call this over all data-points to multiply the entire 200 Data objects involved. 201 202 \param left,right - vectors containing the datapoints 203 \param leftShape,rightShape - shapes of datapoints in the vectors 204 \param leftOffset,rightOffset - beginnings of datapoints in the vectors 205 \param result - Vector to store the resulting datapoint in 206 \param resultShape - expected shape of the resulting datapoint 207 */ 208 ESCRIPT_DLL_API 209 void 210 matMult(const DataTypes::RealVectorType& left, 211 const DataTypes::ShapeType& leftShape, 212 DataTypes::RealVectorType::size_type leftOffset, 213 const DataTypes::RealVectorType& right, 214 const DataTypes::ShapeType& rightShape, 215 DataTypes::RealVectorType::size_type rightOffset, 216 DataTypes::RealVectorType& result, 217 const DataTypes::ShapeType& resultShape); 218 // Hmmmm why is there no offset for the result?? 219 220 221 222 223 /** 224 \brief 225 Determine the shape of the result array for a matrix multiplication 226 of the given views. 227 228 \param left,right - shapes of the left and right matricies 229 \return the shape of the matrix which would result from multiplying left and right 230 */ 231 ESCRIPT_DLL_API 232 DataTypes::ShapeType 233 determineResultShape(const DataTypes::ShapeType& left, 234 const DataTypes::ShapeType& right); 235 236 /** 237 \brief 238 computes a symmetric matrix from your square matrix A: (A + transpose(A)) / 2 239 240 \param in - vector containing the matrix A 241 \param inShape - shape of the matrix A 242 \param inOffset - the beginning of A within the vector in 243 \param ev - vector to store the output matrix 244 \param evShape - expected shape of the output matrix 245 \param evOffset - starting location for storing ev in vector ev 246 */ 247 ESCRIPT_DLL_API 248 inline 249 void 250 symmetric(const DataTypes::RealVectorType& in, 251 const DataTypes::ShapeType& inShape, 252 DataTypes::RealVectorType::size_type inOffset, 253 DataTypes::RealVectorType& ev, 254 const DataTypes::ShapeType& evShape, 255 DataTypes::RealVectorType::size_type evOffset) 256 { 257 if (DataTypes::getRank(inShape) == 2) { 258 int i0, i1; 259 int s0=inShape[0]; 260 int s1=inShape[1]; 261 for (i0=0; i0 868 inline 869 bool 870 checkOffset(const VEC& data, 871 const DataTypes::ShapeType& shape, 872 typename VEC::size_type offset) 873 { 874 return (data.size() >= (offset+DataTypes::noValues(shape))); 875 } 876 877 template 878 inline 879 void 880 unaryOp(DataTypes::RealVectorType& data, const DataTypes::ShapeType& shape, 881 DataTypes::RealVectorType::size_type offset, 882 UnaryFunction operation) 883 { 884 EsysAssert((data.size()>0)&&checkOffset(data,shape,offset), 885 "Error - Couldn't perform unaryOp due to insufficient storage."); 886 DataTypes::RealVectorType::size_type nVals=DataTypes::noValues(shape); 887 for (DataTypes::RealVectorType::size_type i=0;i 894 inline 895 void 896 binaryOp(DataTypes::RealVectorType& left, 897 const DataTypes::ShapeType& leftShape, 898 DataTypes::RealVectorType::size_type leftOffset, 899 const DataTypes::RealVectorType& right, 900 const DataTypes::ShapeType& rightShape, 901 DataTypes::RealVectorType::size_type rightOffset, 902 BinaryFunction operation) 903 { 904 EsysAssert(leftShape==rightShape, 905 "Error - Couldn't perform binaryOp due to shape mismatch,"); 906 EsysAssert(((left.size()>0)&&checkOffset(left,leftShape, leftOffset)), 907 "Error - Couldn't perform binaryOp due to insufficient storage in left object."); 908 EsysAssert(((right.size()>0)&&checkOffset(right,rightShape,rightOffset)), 909 "Error - Couldn't perform binaryOp due to insufficient storage in right object."); 910 for (DataTypes::RealVectorType::size_type i=0;i 916 inline 917 void 918 binaryOp(DataTypes::RealVectorType& left, 919 const DataTypes::ShapeType& leftShape, 920 DataTypes::RealVectorType::size_type offset, 921 double right, 922 BinaryFunction operation) 923 { 924 EsysAssert(((left.size()>0)&&checkOffset(left,leftShape,offset)), 925 "Error - Couldn't perform binaryOp due to insufficient storage in left object."); 926 for (DataTypes::RealVectorType::size_type i=0;i 937 inline 938 void 939 binaryOpVectorHelper(LVEC& left, 940 const DataTypes::ShapeType& leftShape, 941 typename LVEC::size_type leftOffset, 942 const RVEC& right, 943 const DataTypes::ShapeType& rightShape, 944 typename RVEC::size_type rightOffset, 945 BinaryFunction operation) 946 { 947 for (DataTypes::RealVectorType::size_type i=0;i 955 inline 956 void 957 binaryOpVector(LVEC& left, 958 const DataTypes::ShapeType& leftShape, 959 typename LVEC::size_type leftOffset, 960 const RVEC& right, 961 const DataTypes::ShapeType& rightShape, 962 typename RVEC::size_type rightOffset, 963 escript::ESFunction operation) 964 { 965 typedef typename LVEC::ElementType ltype; 966 typedef typename RVEC::ElementType rtype; 967 EsysAssert(leftShape==rightShape, 968 "Error - Couldn't perform binaryOp due to shape mismatch,"); 969 EsysAssert(((left.size()>0)&&checkOffset(left,leftShape, leftOffset)), 970 "Error - Couldn't perform binaryOp due to insufficient storage in left object."); 971 EsysAssert(((right.size()>0)&&checkOffset(right,rightShape,rightOffset)), 972 "Error - Couldn't perform binaryOp due to insufficient storage in right object."); 973 switch (operation) 974 { 975 case POWF:binaryOpVectorHelper(left, leftShape, leftOffset, right, rightShape, rightOffset, pow_func()); break; 976 case PLUSF: binaryOpVectorHelper(left, leftShape, leftOffset, right, rightShape, rightOffset, plus_func()); break; 977 case MINUSF:binaryOpVectorHelper(left, leftShape, leftOffset, right, rightShape, rightOffset, minus_func()); break; 978 case MULTIPLIESF:binaryOpVectorHelper(left, leftShape, leftOffset, right, rightShape, rightOffset, multiplies_func()); break; 979 case DIVIDESF:binaryOpVectorHelper(left, leftShape, leftOffset, right, rightShape, rightOffset, divides_func()); break; 980 case LESSF: 981 case GREATERF: 982 case GREATER_EQUALF: 983 case LESS_EQUALF: 984 default: 985 throw DataException("Unsupported binary operation"); 986 } 987 } 988 989 990 template 991 inline 992 void 993 binaryOpVectorHelper(LVEC& left, 994 const DataTypes::ShapeType& leftShape, 995 typename LVEC::size_type offset, 996 SCALAR right, 997 BinaryFunction operation) 998 { 999 for (DataTypes::RealVectorType::size_type i=0;i 1006 inline 1007 void 1008 binaryOpVector(LVEC& left, 1009 const DataTypes::ShapeType& leftShape, 1010 typename LVEC::size_type offset, 1011 SCALAR right, 1012 escript::ESFunction operation) 1013 { 1014 typedef typename LVEC::ElementType ltype; 1015 typedef SCALAR rtype; 1016 EsysAssert(((left.size()>0)&&checkOffset(left,leftShape,offset)), 1017 "Error - Couldn't perform binaryOp due to insufficient storage in left object."); 1018 switch (operation) 1019 { 1020 case POWF: binaryOpVectorHelper(left, leftShape, offset, right, pow_func()); break; 1021 case PLUSF: binaryOpVectorHelper(left, leftShape, offset, right, plus_func()); break; 1022 case MINUSF:binaryOpVectorHelper(left, leftShape, offset, right, minus_func()); break; 1023 case MULTIPLIESF:binaryOpVectorHelper(left, leftShape, offset, right, multiplies_func()); break; 1024 case DIVIDESF:binaryOpVectorHelper(left, leftShape, offset, right, divides_func()); break; 1025 case LESSF: 1026 case GREATERF: 1027 case GREATER_EQUALF: 1028 case LESS_EQUALF: 1029 default: 1030 throw DataException("Unsupported binary operation"); 1031 } 1032 } 1033 1034 1035 // ------------------- 1036 1037 1038 template 1039 inline 1040 DataTypes::real_t 1041 reductionOp(const DataTypes::RealVectorType& left, 1042 const DataTypes::ShapeType& leftShape, 1043 DataTypes::RealVectorType::size_type offset, 1044 BinaryFunction operation, 1045 DataTypes::real_t initial_value) 1046 { 1047 EsysAssert(((left.size()>0)&&checkOffset(left,leftShape,offset)), 1048 "Error - Couldn't perform reductionOp due to insufficient storage."); 1049 DataTypes::real_t current_value=initial_value; 1050 for (DataTypes::RealVectorType::size_type i=0;i 1057 inline 1058 DataTypes::real_t 1059 reductionOp(const DataTypes::CplxVectorType& left, 1060 const DataTypes::ShapeType& leftShape, 1061 DataTypes::CplxVectorType::size_type offset, 1062 BinaryFunction operation, 1063 DataTypes::real_t initial_value) 1064 { 1065 EsysAssert(((left.size()>0)&&checkOffset(left,leftShape,offset)), 1066 "Error - Couldn't perform reductionOp due to insufficient storage."); 1067 DataTypes::real_t current_value=initial_value; 1068 for (DataTypes::RealVectorType::size_type i=0;i