Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
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
32namespace rheolef {
33// -----------------------------------------------------
34// eigen utilities
35// -----------------------------------------------------
36// convert to sparse
37template <class T>
38void
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}
58template <class T>
59T
60cond (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}
66template <class T, int Nrows, int Ncols>
67bool
68invert (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}
76template <class T>
77void
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}
90template <class T>
91void
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 invert(tiny_matrix< T > &a, tiny_matrix< T > &inv_a)
Definition tiny_lu.h:127
T cond(const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &a)
Definition eigen_util.h:60
void put_matrix_market(std::ostream &out, const Eigen::SparseMatrix< T, Eigen::RowMajor > &a)
Definition eigen_util.h:78
Definition sphere.icc:25