Rheolef  7.2
an efficient C++ finite element environment
eigen_util.h
Go to the documentation of this file.
1 #ifndef _RHEOLEF_EIGEN_UTIL_H
2 #define _RHEOLEF_EIGEN_UTIL_H
23 //
24 // convert a dense matrix to sparse one
25 //
26 // author: Pierre.Saramito@imag.fr
27 //
28 // date: 2 september 2017
29 //
30 #include "rheolef/compiler_eigen.h"
31 
32 namespace rheolef {
33 // -----------------------------------------------------
34 // eigen utilities
35 // -----------------------------------------------------
36 // convert to sparse
37 template <class T>
38 void
40  const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& a,
41  Eigen::SparseMatrix<T,Eigen::RowMajor>& as)
42 {
43  T a_max = 0;
44  for (size_t i = 0; i < size_t(a.rows()); ++i) {
45  for (size_t j = 0; j < size_t(a.cols()); ++j) {
46  a_max = std::max (a_max, abs(a(i,j)));
47  }}
48  T eps = a_max*std::numeric_limits<T>::epsilon();
49  as.resize (a.rows(), a.cols());
50  as.setZero();
51  for (size_t i = 0; i < size_t(a.rows()); ++i) {
52  for (size_t j = 0; j < size_t(a.cols()); ++j) {
53  if (abs(a(i,j)) > eps) {
54  as.coeffRef (i,j) = a(i,j);
55  }
56  }}
57 }
58 template <class T>
59 T
60 cond (const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& a) {
61  if (a.rows() == 0 || a.cols() == 0) return 0;
62  Eigen::JacobiSVD<Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> > svd(a);
63  return svd.singularValues()(0)
64  / svd.singularValues()(svd.singularValues().size()-1);
65 }
66 template <class T, int Nrows, int Ncols>
67 bool
68 invert (const Eigen::Matrix<T,Nrows,Ncols>& a, Eigen::Matrix<T,Nrows,Ncols>& inv_a) {
69  using namespace Eigen;
70  FullPivLU <Matrix<T,Nrows,Ncols> > lu_a (a);
71  if (! lu_a.isInvertible()) return false;
72  Matrix<T,Nrows,Ncols> id = Matrix<T,Nrows,Ncols>::Identity (a.rows(),a.cols());
73  inv_a = lu_a.solve (id);
74  return true;
75 }
76 template <class T>
77 void
79  std::ostream& out,
80  const Eigen::SparseMatrix<T,Eigen::RowMajor>& a)
81 {
82  out << "%%MatrixMarket matrix coordinate real general" << std::endl
83  << a.rows() << " " << a.cols() << " " << a.nonZeros() << std::endl;
84  for (size_t i = 0; i < size_t(a.rows()); ++i) {
85  for (typename Eigen::SparseMatrix<T,Eigen::RowMajor>::InnerIterator p(a,i); p; ++p) {
86  out << i+1 << " " << p.index()+1 << " " << p.value() << std::endl;
87  }
88  }
89 }
90 template <class T>
91 void
93  std::ostream& out,
94  const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& a)
95 {
96  out << "%%MatrixMarket matrix coordinate real general" << std::endl
97  << a.rows() << " " << a.cols() << " " << a.rows()*a.cols() << std::endl;
98  for (size_t i = 0; i < size_t(a.rows()); ++i) {
99  for (size_t j = 0; j < size_t(a.cols()); ++j) {
100  out << i+1 << " " << j+1 << " " << a(i,j) << std::endl;
101  }}
102 }
103 
104 }// namespace rheolef
105 #endif // _RHEOLEF_EIGEN_UTIL_H
Expr1::float_type T
Definition: field_expr.h:230
This file is part of Rheolef.
void eigen_dense2sparse(const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &a, Eigen::SparseMatrix< T, Eigen::RowMajor > &as)
Definition: eigen_util.h:39
void put_matrix_market(std::ostream &out, const Eigen::SparseMatrix< T, Eigen::RowMajor > &a)
Definition: eigen_util.h:78
T cond(const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &a)
Definition: eigen_util.h:60
void invert(tiny_matrix< T > &a, tiny_matrix< T > &inv_a)
Definition: tiny_lu.h:127
Definition: sphere.icc:25
Float epsilon