From 9c36d71621fe64c1c0ed3785840fbc49a4edcc42 Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Wed, 29 May 2002 10:58:56 +0000 Subject: [PATCH] vector initialization in 1.47 fixed git-svn-id: https://svn.dealii.org/trunk@5924 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/solver_gmres.h | 192 +++++++++++++++++++------ 1 file changed, 146 insertions(+), 46 deletions(-) diff --git a/deal.II/lac/include/lac/solver_gmres.h b/deal.II/lac/include/lac/solver_gmres.h index e20129f2bc..9509af21e2 100644 --- a/deal.II/lac/include/lac/solver_gmres.h +++ b/deal.II/lac/include/lac/solver_gmres.h @@ -158,25 +158,142 @@ class SolverGMRES : public Solver void givens_rotation (Vector &h, Vector &b, Vector &ci, Vector &si, int col) const; - /** - * Projected system matrix - */ - FullMatrix H; - /** - * Auxiliary matrix for inverting @p{H} - */ - FullMatrix H1; - + /** + * Projected system matrix + */ + FullMatrix H; + /** + * Auxiliary matrix for inverting @p{H} + */ + FullMatrix H1; + private: - /** - * No copy constructor. - */ - SolverGMRES (const SolverGMRES&); + /** + * Class to hold temporary + * vectors. This class + * automatically allocates a new + * vector, once it is needed. + * + * A future version should also + * be able to shift through + * vectors automatically, + * avoiding restart. + */ + + class TmpVectors + { + public: + /** + * Constructor. Prepares an + * array of @p{VECTOR} of + * length @p{max_size}. + */ + TmpVectors(unsigned int max_size, + VectorMemory& vmem); + /** + * Delete all allocated vectors. + */ + ~TmpVectors(); + + /** + * Get vector number + * @p{i}. If this vector was + * unused before, an error + * occurs. + */ + VECTOR& operator[] (unsigned int i) const; + + /** + * Get vector number + * @p{i}. Allocate it if + * necessary. + * + * If a vector must be + * allocated, @p{temp} is + * used to reinit it to the + * proper dimensions. + */ + VECTOR& operator() (unsigned int i, + const VECTOR& temp); + + private: + /** + * Pool were vectors are obtained from. + */ + VectorMemory& mem; + /** + * Field for storing the vectors. + */ + std::vector data; + /** + * Offset of the first + * vector. This is for later + * when vector rotation will + * be implemented. + */ + unsigned int offset; + }; + /** + * No copy constructor. + */ + SolverGMRES (const SolverGMRES&); }; /* --------------------- Inline and template functions ------------------- */ +template +inline +SolverGMRES::TmpVectors::TmpVectors (unsigned int max_size, + VectorMemory& vmem) + : + mem(vmem), + data (max_size, 0), + offset(0) +{} + + +template +inline +SolverGMRES::TmpVectors::~TmpVectors () +{ + for (typename std::vector::iterator v = data.begin(); + v != data.end(); ++v) + if (*v != 0) + mem.free(*v); +} + + +template +inline VECTOR& +SolverGMRES::TmpVectors::operator[] (unsigned int i) const +{ + Assert (i+offset +inline VECTOR& +SolverGMRES::TmpVectors::operator() (unsigned int i, + const VECTOR& temp) +{ + Assert (i+offsetreinit(temp); + } + return *data[i]; +} + template SolverGMRES::SolverGMRES (SolverControl &cn, @@ -238,22 +355,10 @@ SolverGMRES::solve (const MATRIX &A, deallog.push("GMRES"); const unsigned int n_tmp_vectors = additional_data.max_n_tmp_vectors; - // allocate an array of n_tmp_vectors - // temporary vectors from the VectorMemory - // object; resize them but do not set their - // values since they will be overwritten soon - // anyway. - - // This is really bad since vectors - // should only be allocated if - // really needed. (GK) - std::vector tmp_vectors (n_tmp_vectors, 0); - for (unsigned int tmp=0; tmpreinit (x, true); - }; - + // Generate an object where basis + // vectors are stored. + TmpVectors tmp_vectors (n_tmp_vectors, memory); + // number of the present iteration; this // number is not reset to zero upon a // restart @@ -282,8 +387,8 @@ unsigned int dim = 0; const bool left_precondition = ! additional_data.right_preconditioning; // define two aliases - VECTOR &v = *tmp_vectors[0]; - VECTOR &p = *tmp_vectors[n_tmp_vectors-1]; + VECTOR &v = tmp_vectors(0, x); + VECTOR &p = tmp_vectors(n_tmp_vectors-1, x); /////////////////////////////////// @@ -339,14 +444,14 @@ unsigned int dim = 0; { ++accumulated_iterations; // yet another alias - VECTOR& vv = *tmp_vectors[inner_iteration+1]; + VECTOR& vv = tmp_vectors(inner_iteration+1, x); if (left_precondition) { - A.vmult(p, *tmp_vectors[inner_iteration]); + A.vmult(p, tmp_vectors[inner_iteration]); precondition.vmult(vv,p); } else { - precondition.vmult(p,*tmp_vectors[inner_iteration]); + precondition.vmult(p, tmp_vectors[inner_iteration]); A.vmult(vv,p); }; @@ -355,8 +460,8 @@ unsigned int dim = 0; /* Orthogonalization */ for (unsigned int i=0 ; i