"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "include/armadillo_bits/newarp_SymEigsSolver_meat.hpp" between
armadillo-10.8.2.tar.xz and armadillo-11.0.0.tar.xz

About: Armadillo is a C++ linear algebra library (matrix maths) aiming towards a good balance between speed and ease of use.

newarp_SymEigsSolver_meat.hpp  (armadillo-10.8.2.tar.xz):newarp_SymEigsSolver_meat.hpp  (armadillo-11.0.0.tar.xz)
skipping to change at line 24 skipping to change at line 24
// See the License for the specific language governing permissions and // See the License for the specific language governing permissions and
// limitations under the License. // limitations under the License.
// ------------------------------------------------------------------------ // ------------------------------------------------------------------------
namespace newarp namespace newarp
{ {
template<typename eT, int SelectionRule, typename OpType> template<typename eT, int SelectionRule, typename OpType>
inline inline
void void
SymEigsSolver<eT, SelectionRule, OpType>::fill_rand(eT* dest, const uword N, con
st uword seed_val)
{
arma_extra_debug_sigprint();
typedef typename std::mt19937_64::result_type seed_type;
local_rng.seed( seed_type(seed_val) );
std::uniform_real_distribution<double> dist(-1.0, +1.0);
for(uword i=0; i < N; ++i) { dest[i] = eT(dist(local_rng)); }
}
template<typename eT, int SelectionRule, typename OpType>
inline
void
SymEigsSolver<eT, SelectionRule, OpType>::factorise_from(uword from_k, uword to_ m, const Col<eT>& fk) SymEigsSolver<eT, SelectionRule, OpType>::factorise_from(uword from_k, uword to_ m, const Col<eT>& fk)
{ {
arma_extra_debug_sigprint(); arma_extra_debug_sigprint();
if(to_m <= from_k) { return; } if(to_m <= from_k) { return; }
fac_f = fk; fac_f = fk;
Col<eT> w(dim_n, arma_zeros_indicator()); Col<eT> w(dim_n, arma_zeros_indicator());
// Norm of f // Norm of f
skipping to change at line 48 skipping to change at line 64
fac_H.tail_cols(ncv - from_k).zeros(); fac_H.tail_cols(ncv - from_k).zeros();
fac_H.submat(span(from_k, ncv - 1), span(0, from_k - 1)).zeros(); fac_H.submat(span(from_k, ncv - 1), span(0, from_k - 1)).zeros();
for(uword i = from_k; i <= to_m - 1; i++) for(uword i = from_k; i <= to_m - 1; i++)
{ {
bool restart = false; bool restart = false;
// If beta = 0, then the next V is not full rank // If beta = 0, then the next V is not full rank
// We need to generate a new residual vector that is orthogonal // We need to generate a new residual vector that is orthogonal
// to the current V, which we call a restart // to the current V, which we call a restart
if(beta < near0) if(beta < near0)
{ {
// // Generate new random vector for fac_f
// blas_int idist = 2;
// blas_int iseed[4] = {1, 3, 5, 7};
// iseed[0] = (i + 100) % 4095;
// blas_int n = dim_n;
// lapack::larnv(&idist, &iseed[0], &n, fac_f.memptr());
// Generate new random vector for fac_f // Generate new random vector for fac_f
blas_int idist = 2; fill_rand(fac_f.memptr(), dim_n, i+1);
blas_int iseed[4] = {1, 3, 5, 7};
iseed[0] = (i + 100) % 4095;
blas_int n = dim_n;
lapack::larnv(&idist, &iseed[0], &n, fac_f.memptr());
// f <- f - V * V' * f, so that f is orthogonal to V // f <- f - V * V' * f, so that f is orthogonal to V
Mat<eT> Vs(fac_V.memptr(), dim_n, i, false); // First i columns Mat<eT> Vs(fac_V.memptr(), dim_n, i, false); // First i columns
Col<eT> Vf = Vs.t() * fac_f; Col<eT> Vf = Vs.t() * fac_f;
fac_f -= Vs * Vf; fac_f -= Vs * Vf;
// beta <- ||f|| // beta <- ||f||
beta = norm(fac_f); beta = norm(fac_f);
restart = true; restart = true;
} }
skipping to change at line 216 skipping to change at line 236
arma_extra_debug_sigprint(); arma_extra_debug_sigprint();
uword nev_new = nev; uword nev_new = nev;
for(uword i = nev; i < ncv; i++) for(uword i = nev; i < ncv; i++)
{ {
if(std::abs(ritz_est(i)) < near0) { nev_new++; } if(std::abs(ritz_est(i)) < near0) { nev_new++; }
} }
// Adjust nev_new, according to dsaup2.f line 677~684 in ARPACK // Adjust nev_new, according to dsaup2.f line 677~684 in ARPACK
nev_new += (std::min)(nconv, (ncv - nev_new) / 2); nev_new += (std::min)(nconv, (ncv - nev_new) / 2);
if(nev_new >= ncv) { nev_new = ncv - 1; } if(nev_new >= ncv) { nev_new = ncv - 1; }
if(nev_new == 1 && ncv >= 6)
{ if(nev_new == 1)
nev_new = ncv / 2;
}
else
if(nev_new == 1 && ncv > 2)
{ {
nev_new = 2; if(ncv >= 6) { nev_new = ncv / 2; }
else if(ncv > 2) { nev_new = 2; }
} }
return nev_new; return nev_new;
} }
template<typename eT, int SelectionRule, typename OpType> template<typename eT, int SelectionRule, typename OpType>
inline inline
void void
SymEigsSolver<eT, SelectionRule, OpType>::retrieve_ritzpair() SymEigsSolver<eT, SelectionRule, OpType>::retrieve_ritzpair()
{ {
skipping to change at line 259 skipping to change at line 277
// We keep this order since the first k values will always be // We keep this order since the first k values will always be
// the wanted collection, no matter k is nev_updated (used in restart()) // the wanted collection, no matter k is nev_updated (used in restart())
// or is nev (used in sort_ritzpair()) // or is nev (used in sort_ritzpair())
if(SelectionRule == EigsSelect::BOTH_ENDS) if(SelectionRule == EigsSelect::BOTH_ENDS)
{ {
std::vector<uword> ind_copy(ind); std::vector<uword> ind_copy(ind);
for(uword i = 0; i < ncv; i++) for(uword i = 0; i < ncv; i++)
{ {
// If i is even, pick values from the left (large values) // If i is even, pick values from the left (large values)
// If i is odd, pick values from the right (small values) // If i is odd, pick values from the right (small values)
if(i % 2 == 0) { ind[i] = ind_copy[i / 2]; } else { ind[i] = ind_copy[ncv
- 1 - i / 2]; } ind[i] = (i % 2 == 0) ? ind_copy[i / 2] : ind_copy[ncv - 1 - i / 2];
} }
} }
// Copy the ritz values and vectors to ritz_val and ritz_vec, respectively // Copy the ritz values and vectors to ritz_val and ritz_vec, respectively
for(uword i = 0; i < ncv; i++) for(uword i = 0; i < ncv; i++)
{ {
ritz_val(i) = evals(ind[i]); ritz_val(i) = evals(ind[i]);
ritz_est(i) = evecs(ncv - 1, ind[i]); ritz_est(i) = evecs(ncv - 1, ind[i]);
} }
for(uword i = 0; i < nev; i++) for(uword i = 0; i < nev; i++)
skipping to change at line 369 skipping to change at line 388
if(abs(fac_f).max() < eps) { fac_f.zeros(); } if(abs(fac_f).max() < eps) { fac_f.zeros(); }
} }
template<typename eT, int SelectionRule, typename OpType> template<typename eT, int SelectionRule, typename OpType>
inline inline
void void
SymEigsSolver<eT, SelectionRule, OpType>::init() SymEigsSolver<eT, SelectionRule, OpType>::init()
{ {
arma_extra_debug_sigprint(); arma_extra_debug_sigprint();
// podarray<eT> init_resid(dim_n);
// blas_int idist = 2; // Uniform(-1, 1)
// blas_int iseed[4] = {1, 3, 5, 7}; // Fixed random seed
// blas_int n = dim_n;
// lapack::larnv(&idist, &iseed[0], &n, init_resid.memptr());
// init(init_resid.memptr());
podarray<eT> init_resid(dim_n); podarray<eT> init_resid(dim_n);
blas_int idist = 2; // Uniform(-1, 1)
blas_int iseed[4] = {1, 3, 5, 7}; // Fixed random seed fill_rand(init_resid.memptr(), dim_n, 0);
blas_int n = dim_n;
lapack::larnv(&idist, &iseed[0], &n, init_resid.memptr());
init(init_resid.memptr()); init(init_resid.memptr());
} }
template<typename eT, int SelectionRule, typename OpType> template<typename eT, int SelectionRule, typename OpType>
inline inline
uword uword
SymEigsSolver<eT, SelectionRule, OpType>::compute(uword maxit, eT tol) SymEigsSolver<eT, SelectionRule, OpType>::compute(uword maxit, eT tol)
{ {
arma_extra_debug_sigprint(); arma_extra_debug_sigprint();
skipping to change at line 418 skipping to change at line 443
SymEigsSolver<eT, SelectionRule, OpType>::eigenvalues() SymEigsSolver<eT, SelectionRule, OpType>::eigenvalues()
{ {
arma_extra_debug_sigprint(); arma_extra_debug_sigprint();
uword nconv = std::count(ritz_conv.begin(), ritz_conv.end(), true); uword nconv = std::count(ritz_conv.begin(), ritz_conv.end(), true);
Col<eT> res(nconv, arma_zeros_indicator()); Col<eT> res(nconv, arma_zeros_indicator());
if(nconv > 0) if(nconv > 0)
{ {
uword j = 0; uword j = 0;
for(uword i = 0; i < nev; i++)
for(uword i=0; i < nev; i++)
{ {
if(ritz_conv[i]) if(ritz_conv[i]) { res(j) = ritz_val(i); j++; }
{
res(j) = ritz_val(i);
j++;
}
} }
} }
return res; return res;
} }
template<typename eT, int SelectionRule, typename OpType> template<typename eT, int SelectionRule, typename OpType>
inline inline
Mat<eT> Mat<eT>
SymEigsSolver<eT, SelectionRule, OpType>::eigenvectors(uword nvec) SymEigsSolver<eT, SelectionRule, OpType>::eigenvectors(uword nvec)
skipping to change at line 447 skipping to change at line 469
uword nconv = std::count(ritz_conv.begin(), ritz_conv.end(), true); uword nconv = std::count(ritz_conv.begin(), ritz_conv.end(), true);
nvec = (std::min)(nvec, nconv); nvec = (std::min)(nvec, nconv);
Mat<eT> res(dim_n, nvec); Mat<eT> res(dim_n, nvec);
if(nvec > 0) if(nvec > 0)
{ {
Mat<eT> ritz_vec_conv(ncv, nvec, arma_zeros_indicator()); Mat<eT> ritz_vec_conv(ncv, nvec, arma_zeros_indicator());
uword j = 0; uword j = 0;
for(uword i = 0; i < nev && j < nvec; i++)
for(uword i=0; i < nev && j < nvec; i++)
{ {
if(ritz_conv[i]) if(ritz_conv[i]) { ritz_vec_conv.col(j) = ritz_vec.col(i); j++; }
{
ritz_vec_conv.col(j) = ritz_vec.col(i);
j++;
}
} }
res = fac_V * ritz_vec_conv; res = fac_V * ritz_vec_conv;
} }
return res; return res;
} }
} // namespace newarp } // namespace newarp
 End of changes. 13 change blocks. 
30 lines changed or deleted 49 lines changed or added

Home  |  About  |  Features  |  All  |  Newest  |  Dox  |  Diffs  |  RSS Feeds  |  Screenshots  |  Comments  |  Imprint  |  Privacy  |  HTTP(S)