portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-04-09 21:50:45 +01:00
Code Releases Activity

Implement analogon to test vector analysis in WMG codebase

Browse Source
This commit is contained in:
Daniel Richtmann 2017-11-29 15:04:32 +01:00
parent f260af546e
commit 4e965c168e
No known key found for this signature in database
GPG Key ID: B33C490AF0772057
1 changed files with 9 additions and 18 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

27
tests/solver/Test_wilson_ddalphaamg.cc
View File

@ -38,6 +38,9 @@ class TestVectorAnalyzer {
public: public:
void operator()(LinearOperatorBase<Field> &Linop, std::vector<Field> const & vectors) void operator()(LinearOperatorBase<Field> &Linop, std::vector<Field> const & vectors)
{ {
// this function corresponds to testvector_analysis_PRECISION from the
// DD-αAMG codebase
std::vector<Field> tmp(4, vectors[0]._grid); // bit hacky? std::vector<Field> tmp(4, vectors[0]._grid); // bit hacky?
Gamma g5(Gamma::Algebra::Gamma5); Gamma g5(Gamma::Algebra::Gamma5);
@ -45,32 +48,20 @@ public:
for (auto i = 0; i < vectors.size(); ++i) { for (auto i = 0; i < vectors.size(); ++i) {
Linop.Op(vectors[i], tmp[3]); // apply_operator_PRECISION( l->vbuf_PRECISION[3], test_vectors[i], &(l->p_PRECISION), l, no_threading ); // output, input Linop.Op(vectors[i], tmp[3]);
tmp[0] = g5 * tmp[3]; // is this the same as coarse_gamma5_PRECISION(tmp[0], tmp[3]) in WMG codebase??? tmp[0] = g5 * tmp[3]; // is this the same as coarse_gamma5_PRECISION?
// // use either these two
// auto lambda = innerProduct(vectors[i], l->vbuf_PRECISION[0]);
// lambda = lambda / innerProduct( vectors[i], vectors[i]);
// or this
auto lambda = innerProduct(vectors[i], tmp[0]) / innerProduct(vectors[i], vectors[i]); auto lambda = innerProduct(vectors[i], tmp[0]) / innerProduct(vectors[i], vectors[i]);
tmp[1] = tmp[0] - lambda * vectors[i]; // vector_PRECISION_saxpy(tmp[1], tmp[0], vectors[i], -lambda); tmp[1] = tmp[0] - lambda * vectors[i];
// auto mu = sqrt(norm2(tmp[1]) / norm2(vectors[i])); // mu = global_norm_PRECISION( l->vbuf_PRECISION[1], 0, l->inner_vector_size, l, no_threading )/global_norm_PRECISION( test_vectors[i], 0, l->inner_vector_size, l, no_threading ); auto mu = ::sqrt(norm2(tmp[1]) / norm2(vectors[i]));
// RealD mu = sqrt(norm2(tmp[1]));
// mu = mu / sqrt(norm2(vectors[i])); // mu = global_norm_PRECISION( l->vbuf_PRECISION[1], 0, l->inner_vector_size, l, no_threading )/global_norm_PRECISION( test_vectors[i], 0, l->inner_vector_size, l, no_threading );
RealD mu = norm2(tmp[1]);
mu = mu / norm2(vectors[i]); // mu = global_norm_PRECISION( l->vbuf_PRECISION[1], 0, l->inner_vector_size, l, no_threading )/global_norm_PRECISION( test_vectors[i], 0, l->inner_vector_size, l, no_threading );
mu = std::sqrt(mu);
std::cout << GridLogMessage << std::setprecision(2) << "vector " << i << ": " std::cout << GridLogMessage << std::setprecision(2) << "vector " << i << ": "
<< "singular value: " << lambda << " singular vector precision: " << mu << std::endl; // printf0("singular value: %+lf%+lfi, singular vector precision: %le\n", (double)creal(lambda), (double)cimag(lambda), (double)mu ); << "singular value: " << lambda
<< " singular vector precision: " << mu << std::endl;
} }
} }
}; };