00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #ifndef OPM_AUTODIFF_HPP_HEADER
00022 #define OPM_AUTODIFF_HPP_HEADER
00023
00024 #include <cmath>
00025
00026 namespace Opm
00027 {
00028
00035 template <typename Scalar>
00036 class AutoDiff
00037 {
00038 public:
00041 static AutoDiff
00042 constant(const Scalar x)
00043 {
00044 return function(x, Scalar(0));
00045 }
00046
00049 static AutoDiff
00050 variable(const Scalar x)
00051 {
00052 return function(x, Scalar(1));
00053 }
00054
00057 static AutoDiff
00058 function(const Scalar x, const Scalar dx)
00059 {
00060 return AutoDiff(x, dx);
00061 }
00062
00063 void
00064 operator +=(const Scalar& rhs)
00065 {
00066 val_ += rhs;
00067 }
00068
00069 void
00070 operator +=(const AutoDiff& rhs)
00071 {
00072 val_ += rhs.val_;
00073 der_ += rhs.der_;
00074 }
00075
00076 void
00077 operator -=(const Scalar& rhs)
00078 {
00079 val_ -= rhs;
00080 }
00081
00082 void
00083 operator -=(const AutoDiff& rhs)
00084 {
00085 val_ -= rhs.val_;
00086 der_ -= rhs.der_;
00087 }
00088
00089 void
00090 operator *=(const Scalar& rhs)
00091 {
00092 val_ *= rhs;
00093 der_ *= rhs;
00094 }
00095
00096 void
00097 operator *=(const AutoDiff& rhs)
00098 {
00099 der_ = der_*rhs.val_ + val_*rhs.der_;
00100 val_ *= rhs.val_;
00101 }
00102
00103 void
00104 operator /=(const Scalar& rhs)
00105 {
00106 val_ /= rhs;
00107 der_ /= rhs;
00108 }
00109
00110 void
00111 operator /=(const AutoDiff& rhs)
00112 {
00113 der_ = (der_*rhs.val_ - val_*rhs.der_) / (rhs.val_ * rhs.val_);
00114 val_ /= rhs.val_;
00115 }
00116
00117 template <class Ostream>
00118 Ostream&
00119 print(Ostream& os) const
00120 {
00121 os << "(x,dx) = (" << val_ << ',' << der_ << ")";
00122
00123 return os;
00124 }
00125
00126 const Scalar val() const { return val_; }
00127 const Scalar der() const { return der_; }
00128
00129 private:
00130 AutoDiff(const Scalar x, const Scalar dx)
00131 : val_(x), der_(dx)
00132 {}
00133
00134 Scalar val_;
00135 Scalar der_;
00136 };
00137
00138
00139 template <class Ostream, typename Scalar>
00140 Ostream&
00141 operator<<(Ostream& os, const AutoDiff<Scalar>& fw)
00142 {
00143 return fw.print(os);
00144 }
00145
00146 template <typename Scalar>
00147 AutoDiff<Scalar>
00148 operator +(const AutoDiff<Scalar>& lhs,
00149 const AutoDiff<Scalar>& rhs)
00150 {
00151 AutoDiff<Scalar> ret = lhs;
00152 ret += rhs;
00153
00154 return ret;
00155 }
00156
00157 template <typename Scalar, typename T>
00158 AutoDiff<Scalar>
00159 operator +(const T lhs,
00160 const AutoDiff<Scalar>& rhs)
00161 {
00162 AutoDiff<Scalar> ret = rhs;
00163 ret += Scalar(lhs);
00164
00165 return ret;
00166 }
00167
00168 template <typename Scalar, typename T>
00169 AutoDiff<Scalar>
00170 operator +(const AutoDiff<Scalar>& lhs,
00171 const T rhs)
00172 {
00173 AutoDiff<Scalar> ret = lhs;
00174 ret += Scalar(rhs);
00175
00176 return ret;
00177 }
00178
00179 template <typename Scalar>
00180 AutoDiff<Scalar>
00181 operator -(const AutoDiff<Scalar>& lhs,
00182 const AutoDiff<Scalar>& rhs)
00183 {
00184 AutoDiff<Scalar> ret = lhs;
00185 ret -= rhs;
00186
00187 return ret;
00188 }
00189
00190 template <typename Scalar, typename T>
00191 AutoDiff<Scalar>
00192 operator -(const T lhs,
00193 const AutoDiff<Scalar>& rhs)
00194 {
00195 return AutoDiff<Scalar>::function(Scalar(lhs) - rhs.val(), -rhs.der());
00196 }
00197
00198 template <typename Scalar, typename T>
00199 AutoDiff<Scalar>
00200 operator -(const AutoDiff<Scalar>& lhs,
00201 const T rhs)
00202 {
00203 AutoDiff<Scalar> ret = lhs;
00204 ret -= Scalar(rhs);
00205
00206 return ret;
00207 }
00208
00209 template <typename Scalar>
00210 AutoDiff<Scalar>
00211 operator *(const AutoDiff<Scalar>& lhs,
00212 const AutoDiff<Scalar>& rhs)
00213 {
00214 AutoDiff<Scalar> ret = lhs;
00215 ret *= rhs;
00216
00217 return ret;
00218 }
00219
00220 template <typename Scalar, typename T>
00221 AutoDiff<Scalar>
00222 operator *(const T lhs,
00223 const AutoDiff<Scalar>& rhs)
00224 {
00225 AutoDiff<Scalar> ret = rhs;
00226 ret *= Scalar(lhs);
00227
00228 return ret;
00229 }
00230
00231 template <typename Scalar, typename T>
00232 AutoDiff<Scalar>
00233 operator *(const AutoDiff<Scalar>& lhs,
00234 const T rhs)
00235 {
00236 AutoDiff<Scalar> ret = lhs;
00237 ret *= Scalar(rhs);
00238
00239 return ret;
00240 }
00241
00242 template <typename Scalar>
00243 AutoDiff<Scalar>
00244 operator /(const AutoDiff<Scalar>& lhs,
00245 const AutoDiff<Scalar>& rhs)
00246 {
00247 AutoDiff<Scalar> ret = lhs;
00248 ret /= rhs;
00249
00250 return ret;
00251 }
00252
00253 template <typename Scalar, typename T>
00254 AutoDiff<Scalar>
00255 operator /(const T lhs,
00256 const AutoDiff<Scalar>& rhs)
00257 {
00258 Scalar a = Scalar(lhs) / rhs.val();
00259 Scalar b = (-Scalar(lhs) / (rhs.val() * rhs.val())) * rhs.der();
00260
00261 return AutoDiff<Scalar>::function(a, b);
00262 }
00263
00264 template <typename Scalar, typename T>
00265 AutoDiff<Scalar>
00266 operator /(const AutoDiff<Scalar>& lhs,
00267 const T rhs)
00268 {
00269 Scalar a = lhs.val() / Scalar(rhs);
00270 Scalar b = lhs.der() / Scalar(rhs);
00271
00272 return AutoDiff<Scalar>::function(a, b);
00273 }
00274
00275 template <typename Scalar>
00276 AutoDiff<Scalar>
00277 cos(const AutoDiff<Scalar>& x)
00278 {
00279 Scalar a = std::cos(x.val());
00280 Scalar b = -std::sin(x.val()) * x.der();
00281
00282 return AutoDiff<Scalar>::function(a, b);
00283 }
00284
00285 template <typename Scalar>
00286 AutoDiff<Scalar>
00287 sqrt(const AutoDiff<Scalar>& x)
00288 {
00289 Scalar a = std::sqrt(x.val());
00290 Scalar b = (Scalar(1.0) / (Scalar(2.0) * a)) * x.der();
00291
00292 return AutoDiff<Scalar>::function(a, b);
00293 }
00294
00295 }
00296
00297 namespace std {
00298 using Opm::cos;
00299 using Opm::sqrt;
00300 }
00301
00302 #endif