From: kanschat Date: Wed, 6 Sep 2006 21:23:45 +0000 (+0000) Subject: add condition number calvulation X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=3b5e14e4f5293c76c80576877bc66c612745a7aa;p=dealii-svn.git add condition number calvulation git-svn-id: https://svn.dealii.org/trunk@13842 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/solver_cg.h b/deal.II/lac/include/lac/solver_cg.h index e8162e26ff..1cf2f8b78b 100644 --- a/deal.II/lac/include/lac/solver_cg.h +++ b/deal.II/lac/include/lac/solver_cg.h @@ -15,6 +15,7 @@ #include +#include #include #include #include @@ -81,13 +82,40 @@ class SolverCG : public Solver * eigenvalue estimates. */ bool log_coefficients; - + /** + * Compute the condition + * number of the projected + * matrix. + * + * @note Requires LAPACK support. + */ + bool compute_condition_number; + /** + * Compute the condition + * number of the projected + * matrix in each step. + * + * @note Requires LAPACK support. + */ + bool compute_all_condition_numbers; + + /** + * Compute all eigenvalues of + * the projected matrix. + * + * @note Requires LAPACK support. + */ + bool compute_eigenvalues; + /** * Constructor. Initialize data * fields. Confer the description of * those. */ - AdditionalData (const bool log_coefficients = false); + AdditionalData (const bool log_coefficients = false, + const bool compute_condition_number = false, + const bool compute_all_condition_numbers = false, + const bool compute_eigenvalues = false); }; /** @@ -185,9 +213,15 @@ class SolverCG : public Solver template inline SolverCG::AdditionalData:: -AdditionalData (const bool log_coefficients) +AdditionalData (const bool log_coefficients, + const bool compute_condition_number, + const bool compute_all_condition_numbers, + const bool compute_eigenvalues) : - log_coefficients (log_coefficients) + log_coefficients (log_coefficients), + compute_condition_number(compute_condition_number), + compute_all_condition_numbers(compute_all_condition_numbers), + compute_eigenvalues(compute_eigenvalues) {} @@ -268,7 +302,18 @@ SolverCG::solve (const MATRIX &A, Vp = this->memory.alloc(); Vz = this->memory.alloc(); VAp = this->memory.alloc(); - + // Should we build the matrix for + // eigenvalue computations? + bool do_eigenvalues = additional_data.compute_condition_number + | additional_data.compute_all_condition_numbers + | additional_data.compute_eigenvalues; + double eigen_beta_alpha = 0; + + // vectors used for eigenvalue + // computations + std::vector diagonal; + std::vector offdiagonal; + try { // define some aliases for simpler access VECTOR& g = *Vr; @@ -339,6 +384,30 @@ SolverCG::solve (const MATRIX &A, if (additional_data.log_coefficients) deallog << "alpha-beta:" << alpha << '\t' << beta << std::endl; + // set up the vectors + // containing the diagonal + // and the off diagonal of + // the projected matrix. + if (do_eigenvalues) + { + diagonal.push_back(1./alpha + eigen_beta_alpha); + eigen_beta_alpha = beta/alpha; + offdiagonal.push_back(std::sqrt(beta)/alpha); + } + + if (additional_data.compute_all_condition_numbers) + { + TridiagonalMatrix T(diagonal.size()); + for (unsigned int i=0;i::solve (const MATRIX &A, cleanup(); throw; } + + // Write eigenvalues or condition number + if (do_eigenvalues) + { + TridiagonalMatrix T(diagonal.size()); + for (unsigned int i=0;i