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)));
48 T eps = a_max*std::numeric_limits<T>::epsilon();
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);
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);
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);
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;