1 #ifndef _RHEOLEF_EIGEN_UTIL_H
2 #define _RHEOLEF_EIGEN_UTIL_H
30 #include "rheolef/compiler_eigen.h"
40 const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
a,
41 Eigen::SparseMatrix<T,Eigen::RowMajor>& as)
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)));
49 as.resize (
a.rows(),
a.cols());
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);
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);
66 template <
class T,
int Nrows,
int Ncols>
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);
80 const Eigen::SparseMatrix<T,Eigen::RowMajor>&
a)
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;
94 const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
a)
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;
This file is part of Rheolef.
void eigen_dense2sparse(const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &a, Eigen::SparseMatrix< T, Eigen::RowMajor > &as)
void put_matrix_market(std::ostream &out, const Eigen::SparseMatrix< T, Eigen::RowMajor > &a)
T cond(const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &a)
void invert(tiny_matrix< T > &a, tiny_matrix< T > &inv_a)