From 6bf1a661166ebb844fe944dbdb6a760a59094778 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Mon, 23 Nov 2009 05:37:07 +0000 Subject: [PATCH] Improve description of set function and say that setting can only be done once if there is no sparsity pattern. git-svn-id: https://svn.dealii.org/trunk@20153 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/solver_cg.h | 44 +++++++++---------- .../lac/include/lac/trilinos_sparse_matrix.h | 33 +++++++++----- 2 files changed, 43 insertions(+), 34 deletions(-) diff --git a/deal.II/lac/include/lac/solver_cg.h b/deal.II/lac/include/lac/solver_cg.h index 81be4a62a7..a13dbd2afc 100644 --- a/deal.II/lac/include/lac/solver_cg.h +++ b/deal.II/lac/include/lac/solver_cg.h @@ -102,7 +102,7 @@ class SolverCG : public Solver * @note Requires LAPACK support. */ bool compute_all_condition_numbers; - + /** * Compute all eigenvalues of * the projected matrix. @@ -110,7 +110,7 @@ class SolverCG : public Solver * @note Requires LAPACK support. */ bool compute_eigenvalues; - + /** * Constructor. Initialize data * fields. Confer the description of @@ -186,7 +186,7 @@ class SolverCG : public Solver VECTOR *Vp; VECTOR *Vz; VECTOR *VAp; - + /** * Within the iteration loop, the * square of the residual vector is @@ -227,7 +227,7 @@ AdditionalData (const bool log_coefficients, compute_all_condition_numbers(compute_all_condition_numbers), compute_eigenvalues(compute_eigenvalues) {} - + template @@ -300,7 +300,7 @@ SolverCG::solve (const MATRIX &A, SolverControl::State conv=SolverControl::iterate; deallog.push("cg"); - + // Memory allocation Vr = this->memory.alloc(); Vp = this->memory.alloc(); @@ -312,17 +312,17 @@ SolverCG::solve (const MATRIX &A, | 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; - VECTOR& h = *Vp; - VECTOR& d = *Vz; + VECTOR& g = *Vr; + VECTOR& h = *Vp; + VECTOR& d = *Vz; VECTOR& Ad = *VAp; // resize the vectors, but do not set // the values since they'd be overwritten @@ -335,7 +335,7 @@ SolverCG::solve (const MATRIX &A, // library int it=0; double res,gh,alpha,beta; - + // compute residual. if vector is // zero, then short-circuit the // full computation @@ -349,42 +349,42 @@ SolverCG::solve (const MATRIX &A, res = g.l2_norm(); conv = this->control().check(0,res); - if (conv) + if (conv) { cleanup(); return; } precondition.vmult(h,g); - + d.equ(-1.,h); - + gh = g*h; - + while (conv == SolverControl::iterate) { it++; A.vmult(Ad,d); - + alpha = d*Ad; alpha = gh/alpha; - + g.add(alpha,Ad); x.add(alpha,d ); res = g.l2_norm(); - + print_vectors(it, x, g, d); - + conv = this->control().check(it,res); if (conv) break; precondition.vmult(h,g); - + beta = gh; gh = g*h; beta = gh/beta; - + if (additional_data.log_coefficients) deallog << "alpha-beta:" << alpha << '\t' << beta << std::endl; // set up the vectors @@ -411,7 +411,7 @@ SolverCG::solve (const MATRIX &A, deallog << "Condition number estimate: " << T.eigenvalue(T.n()-1)/T.eigenvalue(0) << std::endl; } - + d.sadd(beta,-1.,h); } } diff --git a/deal.II/lac/include/lac/trilinos_sparse_matrix.h b/deal.II/lac/include/lac/trilinos_sparse_matrix.h index 2b64416e3a..6fb4356d7f 100644 --- a/deal.II/lac/include/lac/trilinos_sparse_matrix.h +++ b/deal.II/lac/include/lac/trilinos_sparse_matrix.h @@ -1160,18 +1160,27 @@ namespace TrilinosWrappers * Set the element (i,j) * to @p value. * - * This function is able to insert - * new elements into the matrix as - * long as compress() has not been - * called, so the sparsity pattern - * will be extended. When compress() - * is called for the first time, then - * this is no longer possible and an - * insertion of elements at positions - * which have not been initialized - * will throw an exception. Moreover, - * if value is not a finite - * number an exception is thrown. + * This function is able to insert new + * elements into the matrix as long as + * compress() has not been called, so + * the sparsity pattern will be + * extended. When compress() is called + * for the first time, then this is no + * longer possible and an insertion of + * elements at positions which have not + * been initialized will throw an + * exception. Note that in case + * elements need to be inserted, it is + * mandatory that elements are inserted + * only once. Otherwise, the elements + * will actually be added in the end + * (since it is not possible to + * efficiently find values to the same + * entry before compress() has been + * called). In the case that an element + * is set more than once, initialize + * the matrix with a sparsity pattern + * first. */ void set (const unsigned int i, const unsigned int j, -- 2.39.5