Rheolef  7.1
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/point.h"
31 
32 #include "rheolef/compiler_eigen.h"
33 
34 namespace rheolef {
35 // -----------------------------------------------------
36 // eigen utilities
37 // -----------------------------------------------------
38 // convert to sparse
39 template <class T>
40 void
42  const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& a,
43  Eigen::SparseMatrix<T,Eigen::RowMajor>& as)
44 {
45  T a_max = 0;
46  for (size_t i = 0; i < size_t(a.rows()); ++i) {
47  for (size_t j = 0; j < size_t(a.cols()); ++j) {
48  a_max = std::max (a_max, abs(a(i,j)));
49  }}
50  T eps = a_max*std::numeric_limits<T>::epsilon();
51  as.resize (a.rows(), a.cols());
52  as.setZero();
53  for (size_t i = 0; i < size_t(a.rows()); ++i) {
54  for (size_t j = 0; j < size_t(a.cols()); ++j) {
55  if (abs(a(i,j)) > eps) {
56  as.coeffRef (i,j) = a(i,j);
57  }
58  }}
59 }
60 template <class T>
61 T
62 cond (const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& a) {
63  Eigen::JacobiSVD<Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> > svd(a);
64  return svd.singularValues()(0)
65  / svd.singularValues()(svd.singularValues().size()-1);
66 }
67 template <class T, int Nrows, int Ncols>
68 bool
69 invert (const Eigen::Matrix<T,Nrows,Ncols>& a, Eigen::Matrix<T,Nrows,Ncols>& inv_a) {
70  using namespace Eigen;
71  FullPivLU <Matrix<T,Nrows,Ncols> > lu_a (a);
72  if (! lu_a.isInvertible()) return false;
73  Matrix<T,Nrows,Ncols> id = Matrix<T,Nrows,Ncols>::Identity (a.rows(),a.cols());
74  inv_a = lu_a.solve (id);
75  return true;
76 }
77 template <class T>
78 void
80  std::ostream& out,
81  const Eigen::SparseMatrix<T,Eigen::RowMajor>& a)
82 {
83  out << "%%MatrixMarket matrix coordinate real general" << std::endl
84  << a.rows() << " " << a.cols() << " " << a.nonZeros() << std::endl;
85  for (size_t i = 0; i < size_t(a.rows()); ++i) {
86  for (typename Eigen::SparseMatrix<T,Eigen::RowMajor>::InnerIterator p(a,i); p; ++p) {
87  out << i << " " << p.index() << " " << p.value() << std::endl;
88  }
89  }
90 }
91 template <class T>
92 void
94  std::ostream& out,
95  const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& a)
96 {
97  out << "%%MatrixMarket matrix coordinate real general" << std::endl
98  << a.rows() << " " << a.cols() << " " << a.rows()*a.cols() << std::endl;
99  for (size_t i = 0; i < size_t(a.rows()); ++i) {
100  for (size_t j = 0; j < size_t(a.cols()); ++j) {
101  out << i << " " << j << " " << a(i,j) << std::endl;
102  }}
103 }
104 
105 }// namespace rheolef
106 #endif // _RHEOLEF_EIGEN_UTIL_H
Expr1::float_type T
Definition: field_expr.h:261
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:41
void put_matrix_market(std::ostream &out, const Eigen::SparseMatrix< T, Eigen::RowMajor > &a)
Definition: eigen_util.h:79
T cond(const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &a)
Definition: eigen_util.h:62
void invert(tiny_matrix< T > &a, tiny_matrix< T > &inv_a)
Definition: tiny_lu.h:127
Definition: sphere.icc:25
Float epsilon