All Classes Namespaces Files Functions Variables Typedefs Enumerator Pages
AutoDiffMatrix.hpp
1 /*
2  Copyright 2014, 2015 SINTEF ICT, Applied Mathematics.
3 
4  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 3 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 */
19 
20 #ifndef OPM_AUTODIFFMATRIX_HEADER_INCLUDED
21 #define OPM_AUTODIFFMATRIX_HEADER_INCLUDED
22 
23 #include <opm/common/utility/platform_dependent/disable_warnings.h>
24 
25 #include <Eigen/Eigen>
26 #include <Eigen/Sparse>
27 
28 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
29 
30 #include <opm/common/ErrorMacros.hpp>
31 #include <opm/autodiff/fastSparseOperations.hpp>
32 #include <vector>
33 
34 
35 namespace Opm
36 {
37 
44  {
45  public:
46  typedef std::vector<double> DiagRep;
47  typedef Eigen::SparseMatrix<double> SparseRep;
48 
49 
54  : type_(Zero),
55  rows_(0),
56  cols_(0),
57  diag_(),
58  sparse_()
59  {
60  }
61 
62 
66  AutoDiffMatrix(const int num_rows, const int num_cols)
67  : type_(Zero),
68  rows_(num_rows),
69  cols_(num_cols),
70  diag_(),
71  sparse_()
72  {
73  }
74 
75 
76 
77 
81  static AutoDiffMatrix createIdentity(const int num_rows_cols)
82  {
83  return AutoDiffMatrix(Identity, num_rows_cols, num_rows_cols);
84  }
85 
86 
87 
91  explicit AutoDiffMatrix(const Eigen::DiagonalMatrix<double, Eigen::Dynamic>& d)
92  : type_(Diagonal),
93  rows_(d.rows()),
94  cols_(d.cols()),
95  diag_(d.diagonal().array().data(), d.diagonal().array().data() + d.rows()),
96  sparse_()
97  {
98  assert(rows_ == cols_);
99  }
100 
101 
102 
106  explicit AutoDiffMatrix(const Eigen::SparseMatrix<double>& s)
107  : type_(Sparse),
108  rows_(s.rows()),
109  cols_(s.cols()),
110  diag_(),
111  sparse_(s)
112  {
113  }
114 
115 
116 
117  AutoDiffMatrix(const AutoDiffMatrix& other) = default;
118  AutoDiffMatrix& operator=(const AutoDiffMatrix& other) = default;
119 
120 
121 
123  : type_(Zero),
124  rows_(0),
125  cols_(0),
126  diag_(),
127  sparse_()
128  {
129  swap(other);
130  }
131 
132 
133 
134  AutoDiffMatrix& operator=(AutoDiffMatrix&& other)
135  {
136  swap(other);
137  return *this;
138  }
139 
140 
141 
142  void swap(AutoDiffMatrix& other)
143  {
144  std::swap(type_, other.type_);
145  std::swap(rows_, other.rows_);
146  std::swap(cols_, other.cols_);
147  diag_.swap(other.diag_);
148  sparse_.swap(other.sparse_);
149  }
150 
151 
152 
160  {
161  assert(rows_ == rhs.rows_);
162  assert(cols_ == rhs.cols_);
163  switch (type_) {
164  case Zero:
165  return rhs;
166  case Identity:
167  switch (rhs.type_) {
168  case Zero:
169  return *this;
170  case Identity:
171  return addII(*this, rhs);
172  case Diagonal:
173  return addDI(rhs, *this);
174  case Sparse:
175  return addSI(rhs, *this);
176  default:
177  OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << rhs.type_);
178  }
179  case Diagonal:
180  switch (rhs.type_) {
181  case Zero:
182  return *this;
183  case Identity:
184  return addDI(*this, rhs);
185  case Diagonal:
186  return addDD(*this, rhs);
187  case Sparse:
188  return addSD(rhs, *this);
189  default:
190  OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << rhs.type_);
191  }
192  case Sparse:
193  switch (rhs.type_) {
194  case Zero:
195  return *this;
196  case Identity:
197  return addSI(*this, rhs);
198  case Diagonal:
199  return addSD(*this, rhs);
200  case Sparse:
201  return addSS(*this, rhs);
202  default:
203  OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << rhs.type_);
204  }
205  default:
206  OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << rhs.type_);
207  }
208  }
209 
210 
211 
212 
213 
214 
221  {
222  assert(cols_ == rhs.rows_);
223  switch (type_) {
224  case Zero:
225  return AutoDiffMatrix(rows_, rhs.cols_);
226  case Identity:
227  return rhs;
228  case Diagonal:
229  switch (rhs.type_) {
230  case Zero:
231  return AutoDiffMatrix(rows_, rhs.cols_);
232  case Identity:
233  return *this;
234  case Diagonal:
235  return mulDD(*this, rhs);
236  case Sparse:
237  return mulDS(*this, rhs);
238  default:
239  OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << rhs.type_);
240  }
241  case Sparse:
242  switch (rhs.type_) {
243  case Zero:
244  return AutoDiffMatrix(rows_, rhs.cols_);
245  case Identity:
246  return *this;
247  case Diagonal:
248  return mulSD(*this, rhs);
249  case Sparse:
250  return mulSS(*this, rhs);
251  default:
252  OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << rhs.type_);
253  }
254  default:
255  OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << rhs.type_);
256  }
257  }
258 
259 
260 
261 
262 
263 
264 
265  AutoDiffMatrix& operator+=(const AutoDiffMatrix& rhs)
266  {
267  if( type_ == Sparse && rhs.type_ == Sparse )
268  {
269  fastSparseAdd( sparse_, rhs.sparse_ );
270  }
271  else {
272  *this = *this + rhs;
273  }
274  return *this;
275  }
276 
277 
278 
279 
280 
281 
282  AutoDiffMatrix& operator-=(const AutoDiffMatrix& rhs)
283  {
284  if( type_ == Sparse && rhs.type_ == Sparse )
285  {
286  fastSparseSubstract( sparse_, rhs.sparse_ );
287  }
288  else {
289  *this = *this + (rhs * -1.0);
290  }
291  return *this;
292  }
293 
294 
295 
296 
297 
298 
304  AutoDiffMatrix operator*(const double rhs) const
305  {
306  switch (type_) {
307  case Zero:
308  return *this;
309  case Identity:
310  {
311  AutoDiffMatrix retval(*this);
312  retval.type_ = Diagonal;
313  retval.diag_.assign(rows_, rhs);
314  return retval;
315  }
316  case Diagonal:
317  {
318  AutoDiffMatrix retval(*this);
319  for (double& elem : retval.diag_) {
320  elem *= rhs;
321  }
322  return retval;
323  }
324  case Sparse:
325  {
326  AutoDiffMatrix retval(*this);
327  retval.sparse_ *= rhs;
328  return retval;
329  }
330  default:
331  OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << type_);
332  }
333  }
334 
335 
336 
337 
338 
339 
345  AutoDiffMatrix operator/(const double rhs) const
346  {
347  switch (type_) {
348  case Zero:
349  return *this;
350  case Identity:
351  {
352  AutoDiffMatrix retval(*this);
353  retval.type_ = Diagonal;
354  retval.diag_.assign(rows_, 1.0/rhs);
355  return retval;
356  }
357  case Diagonal:
358  {
359  AutoDiffMatrix retval(*this);
360  for (double& elem : retval.diag_) {
361  elem /= rhs;
362  }
363  return retval;
364  }
365  case Sparse:
366  {
367  AutoDiffMatrix retval(*this);
368  retval.sparse_ /= rhs;
369  return retval;
370  }
371  default:
372  OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << type_);
373  }
374  }
375 
376 
377 
378 
379 
380 
386  Eigen::VectorXd operator*(const Eigen::VectorXd& rhs) const
387  {
388  assert(cols_ == rhs.size());
389  switch (type_) {
390  case Zero:
391  return Eigen::VectorXd::Zero(rows_);
392  case Identity:
393  return rhs;
394  case Diagonal:
395  {
396  const Eigen::VectorXd d = Eigen::Map<const Eigen::VectorXd>(diag_.data(), rows_);
397  return d.cwiseProduct(rhs);
398  }
399  case Sparse:
400  return sparse_ * rhs;
401  default:
402  OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << type_);
403  }
404  }
405 
406 
407 
408 
409 
410  // Add identity to identity
411  static AutoDiffMatrix addII(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
412  {
413  assert(lhs.type_ == Identity);
414  assert(rhs.type_ == Identity);
415  AutoDiffMatrix retval;
416  retval.type_ = Diagonal;
417  retval.rows_ = lhs.rows_;
418  retval.cols_ = rhs.cols_;
419  retval.diag_.assign(lhs.rows_, 2.0);
420  return retval;
421  }
422 
423  // Add diagonal to identity
424  static AutoDiffMatrix addDI(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
425  {
426  static_cast<void>(rhs); // Silence release-mode warning.
427  assert(lhs.type_ == Diagonal);
428  assert(rhs.type_ == Identity);
429  AutoDiffMatrix retval = lhs;
430  for (int r = 0; r < lhs.rows_; ++r) {
431  retval.diag_[r] += 1.0;
432  }
433  return retval;
434  }
435 
436  // Add diagonal to diagonal
437  static AutoDiffMatrix addDD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
438  {
439  assert(lhs.type_ == Diagonal);
440  assert(rhs.type_ == Diagonal);
441  AutoDiffMatrix retval = lhs;
442  for (int r = 0; r < lhs.rows_; ++r) {
443  retval.diag_[r] += rhs.diag_[r];
444  }
445  return retval;
446  }
447 
448  // Add sparse to identity
449  static AutoDiffMatrix addSI(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
450  {
451  static_cast<void>(rhs); // Silence release-mode warning.
452  assert(lhs.type_ == Sparse);
453  assert(rhs.type_ == Identity);
454  AutoDiffMatrix retval = lhs;
455  retval.sparse_ += spdiag(Eigen::VectorXd::Ones(lhs.rows_));;
456  return retval;
457  }
458 
459  // Add sparse to diagonal
460  static AutoDiffMatrix addSD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
461  {
462  assert(lhs.type_ == Sparse);
463  assert(rhs.type_ == Diagonal);
464  AutoDiffMatrix retval = lhs;
465  retval.sparse_ += spdiag(rhs.diag_);
466  return retval;
467  }
468 
469  // Add sparse to sparse
470  static AutoDiffMatrix addSS(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
471  {
472  assert(lhs.type_ == Sparse);
473  assert(rhs.type_ == Sparse);
474  AutoDiffMatrix retval = lhs;
475  retval.sparse_ += rhs.sparse_;
476  return retval;
477  }
478 
479 
480 
481 
482  // Multiply diagonal with diagonal
483  static AutoDiffMatrix mulDD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
484  {
485  assert(lhs.type_ == Diagonal);
486  assert(rhs.type_ == Diagonal);
487  AutoDiffMatrix retval = lhs;
488  for (int r = 0; r < lhs.rows_; ++r) {
489  retval.diag_[r] *= rhs.diag_[r];
490  }
491  return retval;
492  }
493 
494  // Multiply diagonal with sparse
495  static AutoDiffMatrix mulDS(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
496  {
497  assert(lhs.type_ == Diagonal);
498  assert(rhs.type_ == Sparse);
499  AutoDiffMatrix retval;
500  retval.type_ = Sparse;
501  retval.rows_ = lhs.rows_;
502  retval.cols_ = rhs.cols_;
503  fastDiagSparseProduct(lhs.diag_, rhs.sparse_, retval.sparse_);
504  return retval;
505  }
506 
507  // Multiply sparse with diagonal
508  static AutoDiffMatrix mulSD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
509  {
510  assert(lhs.type_ == Sparse);
511  assert(rhs.type_ == Diagonal);
512  AutoDiffMatrix retval;
513  retval.type_ = Sparse;
514  retval.rows_ = lhs.rows_;
515  retval.cols_ = rhs.cols_;
516  fastSparseDiagProduct(lhs.sparse_, rhs.diag_, retval.sparse_);
517  return retval;
518  }
519 
520  // Multiply sparse with sparse
521  static AutoDiffMatrix mulSS(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
522  {
523  assert(lhs.type_ == Sparse);
524  assert(rhs.type_ == Sparse);
525  AutoDiffMatrix retval;
526  retval.type_ = Sparse;
527  retval.rows_ = lhs.rows_;
528  retval.cols_ = rhs.cols_;
529  fastSparseProduct(lhs.sparse_, rhs.sparse_, retval.sparse_);
530  return retval;
531  }
532 
533 
534 
535 
541  template<class Scalar, int Options, class Index>
542  void toSparse(Eigen::SparseMatrix<Scalar, Options, Index>& s) const
543  {
544  switch (type_) {
545  case Zero:
546  s = Eigen::SparseMatrix<Scalar, Options, Index>(rows_, cols_);
547  return;
548  case Identity:
549  s = spdiag(Eigen::VectorXd::Ones(rows_));
550  return;
551  case Diagonal:
552  s = spdiag(diag_);
553  return;
554  case Sparse:
555  s = sparse_;
556  return;
557  }
558  }
559 
560 
561 
565  int rows() const
566  {
567  return rows_;
568  }
569 
570 
571 
575  int cols() const
576  {
577  return cols_;
578  }
579 
580 
581 
588  int nonZeros() const
589  {
590  switch (type_) {
591  case Zero:
592  return 0;
593  case Identity:
594  return rows_;
595  case Diagonal:
596  return rows_;
597  case Sparse:
598  return sparse_.nonZeros();
599  default:
600  OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << type_);
601  }
602  }
603 
604 
605 
606 
610  double coeff(const int row, const int col) const
611  {
612  switch (type_) {
613  case Zero:
614  return 0.0;
615  case Identity:
616  return (row == col) ? 1.0 : 0.0;
617  case Diagonal:
618  return (row == col) ? diag_[row] : 0.0;
619  case Sparse:
620  return sparse_.coeff(row, col);
621  default:
622  OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << type_);
623  }
624  }
625 
626 
627 
628 
629 
635  const SparseRep& getSparse() const {
636  if (type_ != Sparse) {
644  SparseRep& mutable_sparse = const_cast<SparseRep&>(sparse_);
645  toSparse(mutable_sparse);
646  }
647  return sparse_;
648  }
649 
650 
651  private:
652  enum AudoDiffMatrixType { Zero, Identity, Diagonal, Sparse };
653 
654  AudoDiffMatrixType type_; //< Type of matrix
655  int rows_; //< Number of rows
656  int cols_; //< Number of columns
657  DiagRep diag_; //< Diagonal representation (only if type==Diagonal)
658  SparseRep sparse_; //< Sparse representation (only if type==Sparse)
659 
660 
661 
665  AutoDiffMatrix(AudoDiffMatrixType type, int rows_arg, int cols_arg,
666  DiagRep diag=DiagRep(), SparseRep sparse=SparseRep())
667  : type_(type),
668  rows_(rows_arg),
669  cols_(cols_arg),
670  diag_(diag),
671  sparse_(sparse)
672  {
673  }
674 
675 
676 
677 
678 
684  template <class V>
685  static inline
686  SparseRep
687  spdiag(const V& d)
688  {
689  const int n = d.size();
690  SparseRep mat(n, n);
691  mat.reserve(Eigen::ArrayXi::Ones(n, 1));
692  for (SparseRep::Index i = 0; i < n; ++i) {
693  if (d[i] != 0.0) {
694  mat.insert(i, i) = d[i];
695  }
696  }
697 
698  return mat;
699  }
700 
701  };
702 
703 
704 
708  inline void fastSparseProduct(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs, AutoDiffMatrix& res)
709  {
710  res = lhs * rhs;
711  }
712 
713 
717  inline void fastSparseProduct(const Eigen::SparseMatrix<double>& lhs, const AutoDiffMatrix& rhs, AutoDiffMatrix& res)
718  {
719  res = AutoDiffMatrix(lhs) * rhs;
720  }
721 
722 
726  inline AutoDiffMatrix operator*(const Eigen::SparseMatrix<double>& lhs, const AutoDiffMatrix& rhs)
727  {
728  AutoDiffMatrix retval;
729  fastSparseProduct(lhs, rhs, retval);
730  return retval;
731  }
732 
733 } // namespace Opm
734 
735 
736 #endif // OPM_AUTODIFFMATRIX_HEADER_INCLUDED
AutoDiffMatrix operator/(const double rhs) const
Divides an AutoDiffMatrix by a scalar.
Definition: AutoDiffMatrix.hpp:345
AutoDiffMatrix(const Eigen::SparseMatrix< double > &s)
Creates a sparse matrix from an Eigen sparse matrix.
Definition: AutoDiffMatrix.hpp:106
AutoDiffMatrix operator*(const AutoDiffMatrix &rhs) const
Multiplies two AutoDiffMatrices.
Definition: AutoDiffMatrix.hpp:220
static AutoDiffMatrix createIdentity(const int num_rows_cols)
Creates an identity matrix with num_rows_cols x num_rows_cols entries.
Definition: AutoDiffMatrix.hpp:81
AutoDiffMatrix is a wrapper class that optimizes matrix operations.
Definition: AutoDiffMatrix.hpp:43
const SparseRep & getSparse() const
Returns the sparse representation of this matrix.
Definition: AutoDiffMatrix.hpp:635
void toSparse(Eigen::SparseMatrix< Scalar, Options, Index > &s) const
Converts the AutoDiffMatrix to an Eigen SparseMatrix.This might be an expensive operation to perform ...
Definition: AutoDiffMatrix.hpp:542
AutoDiffMatrix(const Eigen::DiagonalMatrix< double, Eigen::Dynamic > &d)
Creates a diagonal matrix from an Eigen diagonal matrix.
Definition: AutoDiffMatrix.hpp:91
int cols() const
Returns number of columns in the matrix.
Definition: AutoDiffMatrix.hpp:575
Eigen::VectorXd operator*(const Eigen::VectorXd &rhs) const
Multiplies an AutoDiffMatrix with a vector.
Definition: AutoDiffMatrix.hpp:386
int rows() const
Returns number of rows in the matrix.
Definition: AutoDiffMatrix.hpp:565
AutoDiffMatrix()
Creates an empty zero matrix.
Definition: AutoDiffMatrix.hpp:53
AutoDiffMatrix operator+(const AutoDiffMatrix &rhs) const
Adds two AutoDiffMatrices.
Definition: AutoDiffMatrix.hpp:159
double coeff(const int row, const int col) const
Returns element (row, col) in the matrix.
Definition: AutoDiffMatrix.hpp:610
int nonZeros() const
Returns number of non-zero elements in the matrix.
Definition: AutoDiffMatrix.hpp:588
void fastSparseProduct(const AutoDiffMatrix &lhs, const AutoDiffMatrix &rhs, AutoDiffMatrix &res)
Utility function to lessen code changes required elsewhere.
Definition: AutoDiffMatrix.hpp:708
AutoDiffMatrix(const int num_rows, const int num_cols)
Creates a zero matrix with num_rows x num_cols entries.
Definition: AutoDiffMatrix.hpp:66
AutoDiffMatrix operator*(const double rhs) const
Multiplies an AutoDiffMatrix with a scalar.
Definition: AutoDiffMatrix.hpp:304