Commit 1af8d689 authored by Tobias WEBER's avatar Tobias WEBER

pseudoinv

parent d5389098
......@@ -4252,7 +4252,7 @@ eigenvec(const t_mat& mat, bool only_evals=false, bool is_symmetric=false, bool
/**
* singular values of a real or complex matrix mat = U * diag{vals} * V^h
* returns [ ok, U, Vt, vals ]
* returns [ ok, U, Vh, vals ]
*/
template<class t_mat, class t_scalar = typename t_mat::value_type, class t_real = m::underlying_value_type<t_scalar>>
std::tuple<bool, t_mat, t_mat, std::vector<t_real>>
......@@ -4307,6 +4307,29 @@ requires m::is_mat<t_mat>
return std::make_tuple(err==0, U, Vh, vals);
}
template<class t_mat>
std::tuple<t_mat, bool> pseudoinv(const t_mat& mat)
requires m::is_mat<t_mat>
{
using t_scalar = typename t_mat::value_type;
using t_real = m::underlying_value_type<t_scalar>;
auto [ ok, U, Vh, vals ] = singval<t_mat>(mat);
auto V = m::herm(Vh);
auto Uh = m::herm(U);
for(t_real& d : vals)
{
if(!m::equals<t_real>(d, t_real(0)))
d = t_real(1)/d;
}
auto diag = m::diag<t_mat>(vals);
return std::make_tuple(m::prod<t_mat>(V, m::prod(diag, Uh)), ok);
}
}
#endif // USE_LAPACK
......
......@@ -45,6 +45,7 @@ int main()
for(std::size_t i=0; i<evals.size(); ++i)
std::cout << "eval: " << evals[i] << ", evec: " << evecs[i] << std::endl;
auto [ok2, U, Vh, vals] = m_la::singval<t_mat_cplx>(Z);
std::cout << "\nok = " << std::boolalpha << ok2 << std::endl;
std::cout << "singvals: ";
......@@ -54,6 +55,13 @@ int main()
std::cout << "U = " << U << "\nVh = " << Vh << std::endl;
std::cout << "diag{vals} * UVh = " << U*m::diag<t_mat_cplx>(vals)*Vh << std::endl;
auto [inva, ok3a] = m_la::pseudoinv<t_mat_cplx>(Z);
auto [invb, ok3b] = m::inv<t_mat_cplx>(Z);
std::cout << "\nok = " << std::boolalpha << ok3a << ", " << ok3b << std::endl;
std::cout << "pseudoinv = " << inva << std::endl;
std::cout << " inv = " << invb << std::endl;
}
{
......@@ -65,6 +73,7 @@ int main()
std::cout << "eval: " << evals_re[i] << " + i*" << evals_im[i]
<< ", evec: " << evecs_re[i] << " +i*" << evecs_im[i] << std::endl;
auto [ok2, U, Vt, vals] = m_la::singval<t_mat>(Z);
std::cout << "\nok = " << std::boolalpha << ok2 << std::endl;
std::cout << "singvals: ";
......@@ -74,7 +83,14 @@ int main()
std::cout << "U = " << U << "\nVt = " << Vt << std::endl;
std::cout << "diag{vals} * UVt = " << U*m::diag<t_mat>(vals)*Vt << std::endl;
}
auto [inva, ok3a] = m_la::pseudoinv<t_mat>(Z);
auto [invb, ok3b] = m::inv<t_mat>(Z);
std::cout << "\nok = " << std::boolalpha << ok3a << ", " << ok3b << std::endl;
std::cout << "pseudoinv = " << inva << std::endl;
std::cout << " inv = " << invb << std::endl;
}
return 0;
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment