eigen_util.h
Go to the documentation of this file.
#ifndef _RHEOLEF_EIGEN_UTIL_H
#define _RHEOLEF_EIGEN_UTIL_H
//
// convert a dense matrix to sparse one
//
// author: Pierre.Saramito@imag.fr
//
// date: 2 september 2017
//
//#include "rheolef/point.h"
#include "rheolef/compiler_eigen.h"
namespace rheolef {
// -----------------------------------------------------
// eigen utilities
// -----------------------------------------------------
// convert to sparse
template <class T>
void
eigen_dense2sparse (
const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& a,
Eigen::SparseMatrix<T,Eigen::RowMajor>& as)
{
T a_max = 0;
for (size_t i = 0; i < size_t(a.rows()); ++i) {
for (size_t j = 0; j < size_t(a.cols()); ++j) {
a_max = std::max (a_max, abs(a(i,j)));
}}
T eps = a_max*std::numeric_limits<T>::epsilon();
as.resize (a.rows(), a.cols());
as.setZero();
for (size_t i = 0; i < size_t(a.rows()); ++i) {
for (size_t j = 0; j < size_t(a.cols()); ++j) {
if (abs(a(i,j)) > eps) {
as.coeffRef (i,j) = a(i,j);
}
}}
}
template <class T>
T
cond (const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& a) {
Eigen::JacobiSVD<Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> > svd(a);
return svd.singularValues()(0)
/ svd.singularValues()(svd.singularValues().size()-1);
}
template <class T, int Nrows, int Ncols>
bool
invert (const Eigen::Matrix<T,Nrows,Ncols>& a, Eigen::Matrix<T,Nrows,Ncols>& inv_a) {
using namespace Eigen;
FullPivLU <Matrix<T,Nrows,Ncols> > lu_a (a);
if (! lu_a.isInvertible()) return false;
Matrix<T,Nrows,Ncols> id = Matrix<T,Nrows,Ncols>::Identity (a.rows(),a.cols());
inv_a = lu_a.solve (id);
return true;
}
template <class T>
void
put_matrix_market (
std::ostream& out,
const Eigen::SparseMatrix<T,Eigen::RowMajor>& a)
{
out << "%%MatrixMarket matrix coordinate real general" << std::endl
<< a.rows() << " " << a.cols() << " " << a.nonZeros() << std::endl;
for (size_t i = 0; i < size_t(a.rows()); ++i) {
for (typename Eigen::SparseMatrix<T,Eigen::RowMajor>::InnerIterator p(a,i); p; ++p) {
out << i << " " << p.index() << " " << p.value() << std::endl;
}
}
}
template <class T>
void
put_matrix_market (
std::ostream& out,
const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& a)
{
out << "%%MatrixMarket matrix coordinate real general" << std::endl
<< a.rows() << " " << a.cols() << " " << a.rows()*a.cols() << std::endl;
for (size_t i = 0; i < size_t(a.rows()); ++i) {
for (size_t j = 0; j < size_t(a.cols()); ++j) {
out << i << " " << j << " " << a(i,j) << std::endl;
}}
}
}// namespace rheolef
#endif // _RHEOLEF_EIGEN_UTIL_H
T
Expr1::float_type T
Definition: field_expr.h:261
Eigen
Definition: field_eigen.h:35
mkgeo_ball.a
int a
Definition: mkgeo_ball.sh:151
rheolef::io::out
@ out
Definition: rheostream.h:167
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::eigen_dense2sparse
void eigen_dense2sparse(const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &a, Eigen::SparseMatrix< T, Eigen::RowMajor > &as)
Definition: eigen_util.h:41
rheolef::put_matrix_market
void put_matrix_market(std::ostream &out, const Eigen::SparseMatrix< T, Eigen::RowMajor > &a)
Definition: eigen_util.h:79
rheolef::cond
T cond(const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &a)
Definition: eigen_util.h:62
rheolef::invert
void invert(tiny_matrix< T > &a, tiny_matrix< T > &inv_a)
Definition: tiny_lu.h:127
a
Definition: diffusion_isotropic.h:25
p
Definition: sphere.icc:25
epsilon
Float epsilon
Definition: transmission_error.cc:25