1 #ifndef _RHEOLEF_EIGEN_UTIL_H
2 #define _RHEOLEF_EIGEN_UTIL_H
32 #include "rheolef/compiler_eigen.h"
42 const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
a,
43 Eigen::SparseMatrix<T,Eigen::RowMajor>& as)
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)));
51 as.resize (
a.rows(),
a.cols());
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);
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);
67 template <
class T,
int Nrows,
int Ncols>
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);
81 const Eigen::SparseMatrix<T,Eigen::RowMajor>&
a)
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;
95 const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
a)
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;
106 #endif // _RHEOLEF_EIGEN_UTIL_H