00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #ifndef OPM_AUTODIFFMATRIX_HEADER_INCLUDED
00021 #define OPM_AUTODIFFMATRIX_HEADER_INCLUDED
00022
00023 #include <opm/common/utility/platform_dependent/disable_warnings.h>
00024
00025 #include <Eigen/Eigen>
00026 #include <Eigen/Sparse>
00027
00028 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
00029
00030 #include <opm/common/ErrorMacros.hpp>
00031 #include <opm/autodiff/fastSparseOperations.hpp>
00032 #include <vector>
00033
00034
00035 namespace Opm
00036 {
00037
00043 class AutoDiffMatrix
00044 {
00045 public:
00046 typedef std::vector<double> DiagRep;
00047 typedef Eigen::SparseMatrix<double> SparseRep;
00048
00049
00053 AutoDiffMatrix()
00054 : type_(Zero),
00055 rows_(0),
00056 cols_(0),
00057 diag_(),
00058 sparse_()
00059 {
00060 }
00061
00062
00066 AutoDiffMatrix(const int num_rows, const int num_cols)
00067 : type_(Zero),
00068 rows_(num_rows),
00069 cols_(num_cols),
00070 diag_(),
00071 sparse_()
00072 {
00073 }
00074
00075
00076
00077
00081 static AutoDiffMatrix createIdentity(const int num_rows_cols)
00082 {
00083 return AutoDiffMatrix(Identity, num_rows_cols, num_rows_cols);
00084 }
00085
00086
00087
00091 explicit AutoDiffMatrix(const Eigen::DiagonalMatrix<double, Eigen::Dynamic>& d)
00092 : type_(Diagonal),
00093 rows_(d.rows()),
00094 cols_(d.cols()),
00095 diag_(d.diagonal().array().data(), d.diagonal().array().data() + d.rows()),
00096 sparse_()
00097 {
00098 assert(rows_ == cols_);
00099 }
00100
00101
00102
00106 explicit AutoDiffMatrix(const Eigen::SparseMatrix<double>& s)
00107 : type_(Sparse),
00108 rows_(s.rows()),
00109 cols_(s.cols()),
00110 diag_(),
00111 sparse_(s)
00112 {
00113 }
00114
00115
00116
00117 AutoDiffMatrix(const AutoDiffMatrix& other) = default;
00118 AutoDiffMatrix& operator=(const AutoDiffMatrix& other) = default;
00119
00120
00121
00122 AutoDiffMatrix(AutoDiffMatrix&& other)
00123 : type_(Zero),
00124 rows_(0),
00125 cols_(0),
00126 diag_(),
00127 sparse_()
00128 {
00129 swap(other);
00130 }
00131
00132
00133
00134 AutoDiffMatrix& operator=(AutoDiffMatrix&& other)
00135 {
00136 swap(other);
00137 return *this;
00138 }
00139
00140
00141
00142 void swap(AutoDiffMatrix& other)
00143 {
00144 std::swap(type_, other.type_);
00145 std::swap(rows_, other.rows_);
00146 std::swap(cols_, other.cols_);
00147 diag_.swap(other.diag_);
00148 sparse_.swap(other.sparse_);
00149 }
00150
00151
00152
00159 AutoDiffMatrix operator+(const AutoDiffMatrix& rhs) const
00160 {
00161 assert(rows_ == rhs.rows_);
00162 assert(cols_ == rhs.cols_);
00163 switch (type_) {
00164 case Zero:
00165 return rhs;
00166 case Identity:
00167 switch (rhs.type_) {
00168 case Zero:
00169 return *this;
00170 case Identity:
00171 return addII(*this, rhs);
00172 case Diagonal:
00173 return addDI(rhs, *this);
00174 case Sparse:
00175 return addSI(rhs, *this);
00176 default:
00177 OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << rhs.type_);
00178 }
00179 case Diagonal:
00180 switch (rhs.type_) {
00181 case Zero:
00182 return *this;
00183 case Identity:
00184 return addDI(*this, rhs);
00185 case Diagonal:
00186 return addDD(*this, rhs);
00187 case Sparse:
00188 return addSD(rhs, *this);
00189 default:
00190 OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << rhs.type_);
00191 }
00192 case Sparse:
00193 switch (rhs.type_) {
00194 case Zero:
00195 return *this;
00196 case Identity:
00197 return addSI(*this, rhs);
00198 case Diagonal:
00199 return addSD(*this, rhs);
00200 case Sparse:
00201 return addSS(*this, rhs);
00202 default:
00203 OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << rhs.type_);
00204 }
00205 default:
00206 OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << rhs.type_);
00207 }
00208 }
00209
00210
00211
00212
00213
00214
00220 AutoDiffMatrix operator*(const AutoDiffMatrix& rhs) const
00221 {
00222 assert(cols_ == rhs.rows_);
00223 switch (type_) {
00224 case Zero:
00225 return AutoDiffMatrix(rows_, rhs.cols_);
00226 case Identity:
00227 return rhs;
00228 case Diagonal:
00229 switch (rhs.type_) {
00230 case Zero:
00231 return AutoDiffMatrix(rows_, rhs.cols_);
00232 case Identity:
00233 return *this;
00234 case Diagonal:
00235 return mulDD(*this, rhs);
00236 case Sparse:
00237 return mulDS(*this, rhs);
00238 default:
00239 OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << rhs.type_);
00240 }
00241 case Sparse:
00242 switch (rhs.type_) {
00243 case Zero:
00244 return AutoDiffMatrix(rows_, rhs.cols_);
00245 case Identity:
00246 return *this;
00247 case Diagonal:
00248 return mulSD(*this, rhs);
00249 case Sparse:
00250 return mulSS(*this, rhs);
00251 default:
00252 OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << rhs.type_);
00253 }
00254 default:
00255 OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << rhs.type_);
00256 }
00257 }
00258
00259
00260
00261
00262
00263
00264
00265 AutoDiffMatrix& operator+=(const AutoDiffMatrix& rhs)
00266 {
00267 if( type_ == Sparse && rhs.type_ == Sparse )
00268 {
00269 fastSparseAdd( sparse_, rhs.sparse_ );
00270 }
00271 else {
00272 *this = *this + rhs;
00273 }
00274 return *this;
00275 }
00276
00277
00278
00279
00280
00281
00282 AutoDiffMatrix& operator-=(const AutoDiffMatrix& rhs)
00283 {
00284 if( type_ == Sparse && rhs.type_ == Sparse )
00285 {
00286 fastSparseSubstract( sparse_, rhs.sparse_ );
00287 }
00288 else {
00289 *this = *this + (rhs * -1.0);
00290 }
00291 return *this;
00292 }
00293
00294
00295
00296
00297
00298
00304 AutoDiffMatrix operator*(const double rhs) const
00305 {
00306 switch (type_) {
00307 case Zero:
00308 return *this;
00309 case Identity:
00310 {
00311 AutoDiffMatrix retval(*this);
00312 retval.type_ = Diagonal;
00313 retval.diag_.assign(rows_, rhs);
00314 return retval;
00315 }
00316 case Diagonal:
00317 {
00318 AutoDiffMatrix retval(*this);
00319 for (double& elem : retval.diag_) {
00320 elem *= rhs;
00321 }
00322 return retval;
00323 }
00324 case Sparse:
00325 {
00326 AutoDiffMatrix retval(*this);
00327 retval.sparse_ *= rhs;
00328 return retval;
00329 }
00330 default:
00331 OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << type_);
00332 }
00333 }
00334
00335
00336
00337
00338
00339
00345 AutoDiffMatrix operator/(const double rhs) const
00346 {
00347 switch (type_) {
00348 case Zero:
00349 return *this;
00350 case Identity:
00351 {
00352 AutoDiffMatrix retval(*this);
00353 retval.type_ = Diagonal;
00354 retval.diag_.assign(rows_, 1.0/rhs);
00355 return retval;
00356 }
00357 case Diagonal:
00358 {
00359 AutoDiffMatrix retval(*this);
00360 for (double& elem : retval.diag_) {
00361 elem /= rhs;
00362 }
00363 return retval;
00364 }
00365 case Sparse:
00366 {
00367 AutoDiffMatrix retval(*this);
00368 retval.sparse_ /= rhs;
00369 return retval;
00370 }
00371 default:
00372 OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << type_);
00373 }
00374 }
00375
00376
00377
00378
00379
00380
00386 Eigen::VectorXd operator*(const Eigen::VectorXd& rhs) const
00387 {
00388 assert(cols_ == rhs.size());
00389 switch (type_) {
00390 case Zero:
00391 return Eigen::VectorXd::Zero(rows_);
00392 case Identity:
00393 return rhs;
00394 case Diagonal:
00395 {
00396 const Eigen::VectorXd d = Eigen::Map<const Eigen::VectorXd>(diag_.data(), rows_);
00397 return d.cwiseProduct(rhs);
00398 }
00399 case Sparse:
00400 return sparse_ * rhs;
00401 default:
00402 OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << type_);
00403 }
00404 }
00405
00406
00407
00408
00409
00410
00411 static AutoDiffMatrix addII(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
00412 {
00413 assert(lhs.type_ == Identity);
00414 assert(rhs.type_ == Identity);
00415 AutoDiffMatrix retval;
00416 retval.type_ = Diagonal;
00417 retval.rows_ = lhs.rows_;
00418 retval.cols_ = rhs.cols_;
00419 retval.diag_.assign(lhs.rows_, 2.0);
00420 return retval;
00421 }
00422
00423
00424 static AutoDiffMatrix addDI(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
00425 {
00426 static_cast<void>(rhs);
00427 assert(lhs.type_ == Diagonal);
00428 assert(rhs.type_ == Identity);
00429 AutoDiffMatrix retval = lhs;
00430 for (int r = 0; r < lhs.rows_; ++r) {
00431 retval.diag_[r] += 1.0;
00432 }
00433 return retval;
00434 }
00435
00436
00437 static AutoDiffMatrix addDD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
00438 {
00439 assert(lhs.type_ == Diagonal);
00440 assert(rhs.type_ == Diagonal);
00441 AutoDiffMatrix retval = lhs;
00442 for (int r = 0; r < lhs.rows_; ++r) {
00443 retval.diag_[r] += rhs.diag_[r];
00444 }
00445 return retval;
00446 }
00447
00448
00449 static AutoDiffMatrix addSI(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
00450 {
00451 static_cast<void>(rhs);
00452 assert(lhs.type_ == Sparse);
00453 assert(rhs.type_ == Identity);
00454 AutoDiffMatrix retval = lhs;
00455 retval.sparse_ += spdiag(Eigen::VectorXd::Ones(lhs.rows_));;
00456 return retval;
00457 }
00458
00459
00460 static AutoDiffMatrix addSD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
00461 {
00462 assert(lhs.type_ == Sparse);
00463 assert(rhs.type_ == Diagonal);
00464 AutoDiffMatrix retval = lhs;
00465 retval.sparse_ += spdiag(rhs.diag_);
00466 return retval;
00467 }
00468
00469
00470 static AutoDiffMatrix addSS(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
00471 {
00472 assert(lhs.type_ == Sparse);
00473 assert(rhs.type_ == Sparse);
00474 AutoDiffMatrix retval = lhs;
00475 retval.sparse_ += rhs.sparse_;
00476 return retval;
00477 }
00478
00479
00480
00481
00482
00483 static AutoDiffMatrix mulDD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
00484 {
00485 assert(lhs.type_ == Diagonal);
00486 assert(rhs.type_ == Diagonal);
00487 AutoDiffMatrix retval = lhs;
00488 for (int r = 0; r < lhs.rows_; ++r) {
00489 retval.diag_[r] *= rhs.diag_[r];
00490 }
00491 return retval;
00492 }
00493
00494
00495 static AutoDiffMatrix mulDS(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
00496 {
00497 assert(lhs.type_ == Diagonal);
00498 assert(rhs.type_ == Sparse);
00499 AutoDiffMatrix retval;
00500 retval.type_ = Sparse;
00501 retval.rows_ = lhs.rows_;
00502 retval.cols_ = rhs.cols_;
00503 fastDiagSparseProduct(lhs.diag_, rhs.sparse_, retval.sparse_);
00504 return retval;
00505 }
00506
00507
00508 static AutoDiffMatrix mulSD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
00509 {
00510 assert(lhs.type_ == Sparse);
00511 assert(rhs.type_ == Diagonal);
00512 AutoDiffMatrix retval;
00513 retval.type_ = Sparse;
00514 retval.rows_ = lhs.rows_;
00515 retval.cols_ = rhs.cols_;
00516 fastSparseDiagProduct(lhs.sparse_, rhs.diag_, retval.sparse_);
00517 return retval;
00518 }
00519
00520
00521 static AutoDiffMatrix mulSS(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
00522 {
00523 assert(lhs.type_ == Sparse);
00524 assert(rhs.type_ == Sparse);
00525 AutoDiffMatrix retval;
00526 retval.type_ = Sparse;
00527 retval.rows_ = lhs.rows_;
00528 retval.cols_ = rhs.cols_;
00529 fastSparseProduct(lhs.sparse_, rhs.sparse_, retval.sparse_);
00530 return retval;
00531 }
00532
00533
00534
00535
00541 template<class Scalar, int Options, class Index>
00542 void toSparse(Eigen::SparseMatrix<Scalar, Options, Index>& s) const
00543 {
00544 switch (type_) {
00545 case Zero:
00546 s = Eigen::SparseMatrix<Scalar, Options, Index>(rows_, cols_);
00547 return;
00548 case Identity:
00549 s = spdiag(Eigen::VectorXd::Ones(rows_));
00550 return;
00551 case Diagonal:
00552 s = spdiag(diag_);
00553 return;
00554 case Sparse:
00555 s = sparse_;
00556 return;
00557 }
00558 }
00559
00560
00561
00565 int rows() const
00566 {
00567 return rows_;
00568 }
00569
00570
00571
00575 int cols() const
00576 {
00577 return cols_;
00578 }
00579
00580
00581
00588 int nonZeros() const
00589 {
00590 switch (type_) {
00591 case Zero:
00592 return 0;
00593 case Identity:
00594 return rows_;
00595 case Diagonal:
00596 return rows_;
00597 case Sparse:
00598 return sparse_.nonZeros();
00599 default:
00600 OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << type_);
00601 }
00602 }
00603
00604
00605
00606
00610 double coeff(const int row, const int col) const
00611 {
00612 switch (type_) {
00613 case Zero:
00614 return 0.0;
00615 case Identity:
00616 return (row == col) ? 1.0 : 0.0;
00617 case Diagonal:
00618 return (row == col) ? diag_[row] : 0.0;
00619 case Sparse:
00620 return sparse_.coeff(row, col);
00621 default:
00622 OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << type_);
00623 }
00624 }
00625
00626
00627
00628
00629
00635 const SparseRep& getSparse() const {
00636 if (type_ != Sparse) {
00644 SparseRep& mutable_sparse = const_cast<SparseRep&>(sparse_);
00645 toSparse(mutable_sparse);
00646 }
00647 return sparse_;
00648 }
00649
00650
00651 private:
00652 enum AudoDiffMatrixType { Zero, Identity, Diagonal, Sparse };
00653
00654 AudoDiffMatrixType type_;
00655 int rows_;
00656 int cols_;
00657 DiagRep diag_;
00658 SparseRep sparse_;
00659
00660
00661
00665 AutoDiffMatrix(AudoDiffMatrixType type, int rows_arg, int cols_arg,
00666 DiagRep diag=DiagRep(), SparseRep sparse=SparseRep())
00667 : type_(type),
00668 rows_(rows_arg),
00669 cols_(cols_arg),
00670 diag_(diag),
00671 sparse_(sparse)
00672 {
00673 }
00674
00675
00676
00677
00678
00684 template <class V>
00685 static inline
00686 SparseRep
00687 spdiag(const V& d)
00688 {
00689 const int n = d.size();
00690 SparseRep mat(n, n);
00691 mat.reserve(Eigen::ArrayXi::Ones(n, 1));
00692 for (SparseRep::Index i = 0; i < n; ++i) {
00693 if (d[i] != 0.0) {
00694 mat.insert(i, i) = d[i];
00695 }
00696 }
00697
00698 return mat;
00699 }
00700
00701 };
00702
00703
00704
00708 inline void fastSparseProduct(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs, AutoDiffMatrix& res)
00709 {
00710 res = lhs * rhs;
00711 }
00712
00713
00717 inline void fastSparseProduct(const Eigen::SparseMatrix<double>& lhs, const AutoDiffMatrix& rhs, AutoDiffMatrix& res)
00718 {
00719 res = AutoDiffMatrix(lhs) * rhs;
00720 }
00721
00722
00726 inline AutoDiffMatrix operator*(const Eigen::SparseMatrix<double>& lhs, const AutoDiffMatrix& rhs)
00727 {
00728 AutoDiffMatrix retval;
00729 fastSparseProduct(lhs, rhs, retval);
00730 return retval;
00731 }
00732
00733 }
00734
00735
00736 #endif // OPM_AUTODIFFMATRIX_HEADER_INCLUDED