open (FILE, $filename);
while (<FILE>) {
-
# if the lines contains a "template" at the
# beginning and no semicolon at the end: join it
# with the next line.
if ( /^\s*template/ && !/;\s*$/ ) {
+ # if this is a multiline template, then join following lines as well
+ # find out by looking at the last char of the line
+ while ( ! />\s*$/ ) {
+ # more concise check: find out whether the number of '<' is
+ # not equal to the number of '>'
+ my ( $tmp1, $tmp2 ) = ($_, $_);
+ $tmp1 =~ s/[^<]//g;
+ $tmp2 =~ s/[^>]//g;
+ if ( length ($tmp1) == length ($tmp2) ) {
+ last;
+ }
+ else
+ {
+ s/\n//;
+ $_ = $_ . " " . <FILE>;
+ s/ \s+/ /g;
+ }
+ }
s/\n//;
+ s/ \s+/ /g;
$_ = $_ . " " . <FILE>;
}
- if ( /^\s*((template\s*<(([-\w,_\s]|<([-\w_,+\s])+>)+)>\s*)?(class|struct))(.*)/ ) {
+ if ( /^\s*((template\s*<
+ (
+ ([-\w,_=\s] |
+ <([-\w_,=\s])+> )+
+ )>\s*)?
+ (class|struct))(.*)/x ) {
# this is the declaration of a class, possibly a template
$basepart = $1;
$rest = $7;
# $rest contains the name of the class and what comes after that
#
# first extract the name of the class
- $rest =~ /([\w_]+(\s*<(([-\w,_\s]|<([-\w,+\s])+>)+)>)?)(.*)/;
+ $rest =~ /([\w_]+(\s*<(([-\w,_\s]|<([-\w,\s])+>)+)>)?)(.*)/;
$name = $1;
$rest = $6;
$declaration =~ s/^\s+//g;
$declaration =~ s/\s+$//g;
+ # strip default template parameters
+ while ( $declaration =~ s/<(.*)=\s[-\w,_\s]+(<[^.]*>)?(.*)>/<$1 $3 > / ) {
+ }
+
# impose a negativ-list of names we do not want to
# have in this file. this is due to a compiler bug in
# egcs 1.1.2, which does not gracefully handle local
* with a subscriptor for smart pointers.
* @author Guido Kanschat, 1999
*/
-template<class MATRIX>
+template<class MATRIX = SparseMatrix<double> >
class MGMatrix : public Subscriptor,
public vector<MATRIX>
{
{
Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected());
- SolverControl control(4000, 1e-16);
- PrimitiveVectorMemory<Vector<double> > memory;
- SolverCG<SparseMatrix<double>,Vector<double> > cg(control,memory);
+ SolverControl control(4000, 1e-16);
+ PrimitiveVectorMemory<> memory;
+ SolverCG<> cg(control,memory);
- PreconditionRelaxation<SparseMatrix<double>, Vector<double> >
+ PreconditionRelaxation<>
prec(system_matrix,
&SparseMatrix<double>::template precondition_SSOR<double>,
1.2);
MatrixTools<dim>::apply_boundary_values (boundary_values,
mass_matrix, vec, tmp);
- SolverControl control(1000,1e-16);
- PrimitiveVectorMemory<Vector<double> > memory;
- SolverCG<SparseMatrix<double>,Vector<double> > cg(control,memory);
+ SolverControl control(1000,1e-16);
+ PrimitiveVectorMemory<> memory;
+ SolverCG<> cg(control,memory);
- PreconditionRelaxation<SparseMatrix<double>, Vector<double> >
+ PreconditionRelaxation<>
prec(mass_matrix,
&SparseMatrix<double>::template precondition_SSOR<double>,
1.2);
Vector<double> boundary_projection (rhs.size());
- SolverControl control(1000, 1e-16);
- PrimitiveVectorMemory<Vector<double> > memory;
- SolverCG<SparseMatrix<double>,Vector<double> > cg(control,memory);
+ SolverControl control(1000, 1e-16);
+ PrimitiveVectorMemory<> memory;
+ SolverCG<> cg(control,memory);
- PreconditionRelaxation<SparseMatrix<double>, Vector<double> >
+ PreconditionRelaxation<>
prec(mass_matrix,
&SparseMatrix<double>::template precondition_SSOR<double>,
1.2);
* with a subscriptor for smart pointers.
* @author Guido Kanschat, 1999
*/
-template<class MATRIX>
+template<class MATRIX = SparseMatrix<double> >
class MGMatrix : public Subscriptor,
public vector<MATRIX>
{
*
* @author Guido Kanschat, Wolfgang Bangerth, 1999
*/
-template<class MATRIX, class VECTOR>
+template<class MATRIX = SparseMatrix<double>,
+ class VECTOR = Vector<double> >
class PreconditionUseMatrix
{
public:
*
* @author Guido Kanschat, Wolfgang Bangerth, 1999
*/
-template<class MATRIX, class VECTOR>
+template<class MATRIX = SparseMatrix<double>,
+ class VECTOR = Vector<double> >
class PreconditionRelaxation
{
public:
*
* @author Guido Kanschat, 1999
*/
-template<class SOLVER, class MATRIX, class PRECONDITION>
+template<class SOLVER,
+ class MATRIX = SparseMatrix<double>,
+ class PRECONDITION = PreconditionIdentity>
class PreconditionLACSolver
{
public:
* store the inverted blocks with less accuracy than the original matrix;
* for example, #number==double, inverse_type=float# might be a viable choice.
*/
-template<typename number, typename inverse_type>
+template<typename number,
+ typename inverse_type = number>
class PreconditionBlock: public Subscriptor
{
public:
* (of arbitray structure) below the diagonal blocks are used
* in the #operator ()# function of this class.
*/
-template<typename number, typename inverse_type>
+template<typename number,
+ typename inverse_type = number>
class PreconditionBlockSOR : public PreconditionBlock<number, inverse_type>
{
public:
*
* @author Ralf Hartmann, 1999
*/
-template <class Matrix, class Vector>
+template <class Matrix = SparseMatrix<double>,
+ class Vector = Vector<double> >
class PreconditionSelector
{
public:
-// forward declaration
-class SolverControl;
-template <class Vector> class VectorMemory;
+#include <lac/forward-declarations.h>
/**
* Base class for iterative solvers.
*
- * HAS TO BE UPDATED!
+//TODO: * HAS TO BE UPDATED!
*
* This class defines possible
* return states of linear solvers and provides interfaces to a memory
*
* @author Wolfgang Bangerth, Guido Kanschat, Ralf Hartmann, 1997, 1998, 1999
*/
-template <class Matrix, class Vector>
+template <class Matrix = SparseMatrix<double>,
+ class Vector = Vector<double> >
class Solver
{
public:
* has a default argument, so you may call it without the additional
* parameter.
*/
-template<class Matrix, class Vector>
+template <class Matrix = SparseMatrix<double>,
+ class Vector = Vector<double> >
class SolverBicgstab : public Solver<Matrix,Vector>
{
public:
*
* @author Original implementation by G. Kanschat, R. Becker and F.-T. Suttmeier, reworking and documentation by Wolfgang Bangerth
*/
-template<class Matrix, class Vector>
+template <class Matrix = SparseMatrix<double>,
+ class Vector = Vector<double> >
class SolverCG : public Solver<Matrix,Vector>
{
public:
*
* @author Wolfgang Bangerth
*/
-template<class Matrix, class VECTOR>
-class SolverGMRES : public Solver<Matrix, VECTOR>
+template <class MATRIX = SparseMatrix<double>,
+ class VECTOR = Vector<double> >
+class SolverGMRES : public Solver<MATRIX, VECTOR>
{
public:
/**
* Solver method.
*/
template<class Preconditioner>
- typename Solver<Matrix,VECTOR>::ReturnState solve (const Matrix &A,
+ typename Solver<MATRIX,VECTOR>::ReturnState solve (const MATRIX &A,
VECTOR &x,
const VECTOR &b,
const Preconditioner& precondition);
/* --------------------- Inline and template functions ------------------- */
-template <class Matrix, class VECTOR>
-SolverGMRES<Matrix,VECTOR>::SolverGMRES (SolverControl &cn,
+template <class MATRIX, class VECTOR>
+SolverGMRES<MATRIX,VECTOR>::SolverGMRES (SolverControl &cn,
VectorMemory<VECTOR> &mem,
const AdditionalData &data) :
- Solver<Matrix,VECTOR> (cn,mem),
+ Solver<MATRIX,VECTOR> (cn,mem),
additional_data(data)
{
Assert (additional_data.max_n_tmp_vectors >= 10,
-template <class Matrix, class VECTOR>
+template <class MATRIX, class VECTOR>
inline
void
-SolverGMRES<Matrix,VECTOR>::givens_rotation (Vector<double> &h,
+SolverGMRES<MATRIX,VECTOR>::givens_rotation (Vector<double> &h,
Vector<double> &b,
Vector<double> &ci,
Vector<double> &si,
-template<class Matrix, class VECTOR>
+template<class MATRIX, class VECTOR>
template<class Preconditioner>
-typename Solver<Matrix,VECTOR>::ReturnState
-SolverGMRES<Matrix,VECTOR>::solve (const Matrix& A,
+typename Solver<MATRIX,VECTOR>::ReturnState
+SolverGMRES<MATRIX,VECTOR>::solve (const MATRIX& A,
VECTOR & x,
const VECTOR& b,
const Preconditioner& precondition)
-template<class Matrix, class VECTOR>
+template<class MATRIX, class VECTOR>
double
-SolverGMRES<Matrix,VECTOR>::criterion ()
+SolverGMRES<MATRIX,VECTOR>::criterion ()
{
// dummy implementation. this function is
// not needed for the present implementation
*
* @author Guido Kanschat, 1999
*/
-template<class Matrix, class Vector>
+template <class Matrix = SparseMatrix<double>,
+ class Vector = Vector<double> >
class SolverQMRS : public Solver<Matrix,Vector>
{
public:
*
* @author Ralf Hartmann
*/
-template<class Matrix, class Vector>
+template <class Matrix = SparseMatrix<double>,
+ class Vector = Vector<double> >
class SolverRichardson : public Solver<Matrix, Vector>
{
public:
*
* @author Ralf Hartmann, 1999
*/
-template <class Matrix, class Vector>
+template <class Matrix = SparseMatrix<double>,
+ class Vector = Vector<double> >
class SolverSelector
{
public:
* sophisticated management of vectors. One of these has to be
* applied by the user according to his needs.
*/
-template<class Vector>
+template<class Vector = Vector<double> >
class VectorMemory : public Subscriptor
{
public:
* vectors as needed from the global heap, i.e. performs no
* specially adapted actions to the purpose of this class.
*/
-template<class Vector>
+template<class Vector = Vector<double> >
class PrimitiveVectorMemory : public VectorMemory<Vector>
{
public:
*
* @author Guido Kanschat, 1999
*/
-
-template<class Vector>
+template<class Vector = Vector<double> >
class GrowingVectorMemory : public VectorMemory<Vector>
{
public:
* of vectors.
*/
GrowingVectorMemory(unsigned int initial_size = 0);
+
/**
* Destructor.
* Release all vectors.
* are allocated vectors left.
*/
~GrowingVectorMemory();
+
/**
* Return new vector from the pool.
*/
* vector itself.
*/
typedef pair<bool, Vector* > entry_type;
+
/**
* Array of allocated vectors.
*/
vector<entry_type> pool;
+
/**
* Overall number of allocations.
*/
unsigned int n_alloc;
};
+
+
+/* --------------------- inline functions ---------------------- */
+
+
template<typename Vector>
-GrowingVectorMemory<Vector>::GrowingVectorMemory(unsigned int
- initial_size)
+GrowingVectorMemory<Vector>::GrowingVectorMemory(const unsigned int initial_size)
: pool(initial_size)
{
for (vector<entry_type>::iterator i=pool.begin();i != pool.end()
n_alloc = 0;
}
+
+
template<typename Vector>
GrowingVectorMemory<Vector>::~GrowingVectorMemory()
{
Vector<double> u;
u.reinit(f);
- PrimitiveVectorMemory<Vector<double> > mem;
+ PrimitiveVectorMemory<> mem;
SolverControl control(100, 1.e-12);
- SolverCG<SparseMatrix<double>, Vector<double> >solver(control, mem);
+ SolverCG<> solver(control, mem);
PreconditionIdentity precondition;
solver.solve(A,u,f,precondition);
equation.build_mgmatrix(mgA, mgdof, quadrature);
SolverControl cgcontrol(10,0., false, false);
- PrimitiveVectorMemory<Vector<double> > cgmem;
- SolverCG<SparseMatrix<double>, Vector<double> > cgsolver(cgcontrol, cgmem);
+ PrimitiveVectorMemory<> cgmem;
+ SolverCG<> cgsolver(cgcontrol, cgmem);
PreconditionIdentity cgprec;
- MGCoarseGridLACIteration<SolverCG<SparseMatrix<double>, Vector<double> >,
- SparseMatrix<double>, PreconditionIdentity>
+ MGCoarseGridLACIteration<SolverCG<>, SparseMatrix<double>, PreconditionIdentity>
coarse(cgsolver, mgA[0], cgprec);
MGSmootherLAC smoother(mgA);
MG<2> multigrid(mgdof, hanging_nodes, mgA, transfer);
- PreconditionMG<MG<2>, Vector<double> >
+ PreconditionMG<MG<2> >
mgprecondition(multigrid, smoother, smoother, coarse);
solver.solve(A, u, f, mgprecondition);
const Vector<double> &rhs) const
{
SolverControl control(2,1.e-300,false,false);
- PrimitiveVectorMemory<Vector<double> > mem;
- SolverRichardson<SparseMatrix<double> , Vector<double> > rich(control, mem);
- PreconditionRelaxation<SparseMatrix<double> , Vector<double> >
+ PrimitiveVectorMemory<> mem;
+ SolverRichardson<> rich(control, mem);
+ PreconditionRelaxation<>
prec((*matrices)[level], &SparseMatrix<double> ::template precondition_SSOR<double>, 1.);
rich.solve((*matrices)[level], u, rhs, prec);
Vector<double> u;
u.reinit(f);
- PrimitiveVectorMemory<Vector<double> > mem;
+ PrimitiveVectorMemory<> mem;
SolverControl control(20, 1.e-12, true);
- SolverRichardson<SparseMatrix<double>, Vector<double> >solver(control, mem);
+ SolverRichardson<> solver(control, mem);
vector<SparseMatrixStruct> mgstruct(tr.n_levels());
MGMatrix<SparseMatrix<double> > mgA(0,tr.n_levels()-1);
equation.build_mgmatrix(mgA, mgdof, quadrature);
SolverControl cgcontrol(20,0., false, false);
- PrimitiveVectorMemory<Vector<double> > cgmem;
- SolverCG<SparseMatrix<double>, Vector<double> > cgsolver(cgcontrol, cgmem);
+ PrimitiveVectorMemory<> cgmem;
+ SolverCG<> cgsolver(cgcontrol, cgmem);
PreconditionIdentity cgprec;
- MGCoarseGridLACIteration<SolverCG<SparseMatrix<double>, Vector<double> >,
- SparseMatrix<double>, PreconditionIdentity>
+ MGCoarseGridLACIteration<SolverCG<>, SparseMatrix<double>, PreconditionIdentity>
coarse(cgsolver, mgA[tr.n_levels()-2], cgprec);
MGSmootherRelaxation<double>
MG<2> multigrid(mgdof, hanging_nodes, mgA, transfer, tr.n_levels()-2);
- PreconditionMG<MG<2>, Vector<double> >
+ PreconditionMG<MG<2> >
mgprecondition(multigrid, smoother, smoother, coarse);
u = 0.;
unsigned int TimeStep_Wave<dim>::solve (const UserMatrix &matrix,
Vector<double> &solution,
const Vector<double> &rhs) const {
- SolverControl control(2000, 1.e-12);
- PrimitiveVectorMemory<Vector<double> > memory;
- SolverCG<UserMatrix,Vector<double> > pcg(control,memory);
+ SolverControl control(2000, 1.e-12);
+ PrimitiveVectorMemory<> memory;
+ SolverCG<UserMatrix> pcg(control,memory);
// solve
pcg.solve (matrix, solution, rhs,
- PreconditionUseMatrix<UserMatrix,Vector<double> >
+ PreconditionUseMatrix<UserMatrix>
(matrix,
&UserMatrix::precondition));
// distribute solution
};
-typedef MGCoarseGridLACIteration<SolverCG<SparseMatrix<double> , Vector<double> >,
-SparseMatrix<double>, /*PreconditionRelaxation<SparseMatrix<double> ,*/
- PreconditionIdentity >
-Coarse;
+typedef MGCoarseGridLACIteration<SolverCG<>,
+ SparseMatrix<double>,
+ PreconditionIdentity > Coarse;
int main()
FDMGTransfer transfer(maxsize, maxsize, maxlevel);
// coarse grid solver
- PrimitiveVectorMemory<Vector<double> > cgmem;
+ PrimitiveVectorMemory<> cgmem;
ReductionControl cgcontrol(100, 1.e-30, 1.e-2, false, false);
- SolverCG<SparseMatrix<double> , Vector<double> > cgcg(cgcontrol,cgmem);
+ SolverCG<> cgcg(cgcontrol,cgmem);
for (unsigned int level = 0; level <= maxlevel; ++level)
FDMG multigrid(level, A, transfer);
-// PreconditionRelaxation<SparseMatrix<double> , Vector<double> >
+// PreconditionRelaxation<>
PreconditionIdentity
cgprec;//(A[minlevel], &SparseMatrix<double> ::template precondition_SSOR<double>, 1.2);
{
for (unsigned int level = minlevel; level<=maxlevel ; ++level)
{
- s[level].reinit(matrices[level].m());
- d[level].reinit(matrices[level].m());
+ solution[level].reinit(matrices[level].m());
+ defect[level].reinit(matrices[level].m());
}
}
void
FDMG::copy_to_mg(const Vector<TYPE>& v)
{
- d[maxlevel] = v;
+ defect[maxlevel] = v;
}
void
FDMG::copy_from_mg(Vector<TYPE>& v)
{
- v = s[maxlevel];
+ v = solution[maxlevel];
}
const Vector<double> &rhs) const
{
SolverControl control(1,1.e-300,false,false);
- PrimitiveVectorMemory<Vector<double> > mem;
- SolverRichardson<SparseMatrix<double> , Vector<double> > rich(control, mem);
- PreconditionRelaxation<SparseMatrix<double> , Vector<double> >
+ PrimitiveVectorMemory<> mem;
+ SolverRichardson<> rich(control, mem);
+ PreconditionRelaxation<>
prec((*matrices)[level], &SparseMatrix<double> ::template precondition_SSOR<double>, 1.);
rich.solve((*matrices)[level], u, rhs, prec);
deallog.attach(logfile);
deallog.depth_console(0);
- GrowingVectorMemory<Vector<double> > mem;
+ GrowingVectorMemory<> mem;
SolverControl control(100, 1.e-5);
SolverControl verbose_control(100, 1.e-5, true);
- SolverCG<SparseMatrix<float> , Vector<double> > cg(control, mem);
- SolverGMRES<SparseMatrix<float> , Vector<double> > gmres(control, mem,20);
- SolverBicgstab<SparseMatrix<float> , Vector<double> > bicgstab(control, mem);
- SolverRichardson<SparseMatrix<float> , Vector<double> > rich(control, mem);
- SolverQMRS<SparseMatrix<float> , Vector<double> > qmrs(control, mem);
+ SolverCG<> cg(control, mem);
+ SolverGMRES<> gmres(control, mem,20);
+ SolverBicgstab<> bicgstab(control, mem);
+ SolverRichardson<> rich(control, mem);
+ SolverQMRS<> qmrs(control, mem);
for (unsigned int size=4; size <= 40; size *= 3)
{
SparseMatrixStruct structure(dim, dim, 5);
testproblem.build_structure(structure);
structure.compress();
- SparseMatrix<float> A(structure);
+ SparseMatrix<double> A(structure);
testproblem.laplacian(A);
PreconditionIdentity
prec_no;
- PreconditionRelaxation<SparseMatrix<float> , Vector<double> >
- prec_sor(A, &SparseMatrix<float> ::template precondition_SOR<double>, 1.2);
- PreconditionRelaxation<SparseMatrix<float> , Vector<double> >
- prec_ssor(A, &SparseMatrix<float> ::template precondition_SSOR<double>, 1.2);
+ PreconditionRelaxation<>
+ prec_sor(A, &SparseMatrix<double>::template precondition_SOR<double>, 1.2);
+ PreconditionRelaxation<>
+ prec_ssor(A, &SparseMatrix<double>::template precondition_SSOR<double>, 1.2);
Vector<double> f(dim);
Vector<double> u(dim);
DEAL:no:cg::Starting
DEAL:no:cg::Convergence step 58
DEAL:no:Bicgstab::Starting
-DEAL:no:Bicgstab::Convergence step 45
+DEAL:no:Bicgstab::Convergence step 44
DEAL:no:GMRES::Starting
DEAL:no:GMRES::Starting
DEAL:no:GMRES::Failure step 100
void
-FDMGTransfer::restrict (const unsigned int from_level,
- Vector<double> &dst,
- const Vector<double> &src) const
+FDMGTransfer::restrict_and_add (const unsigned int from_level,
+ Vector<double> &dst,
+ const Vector<double> &src) const
{
Assert((from_level>0) && (from_level<=matrices.size()),
ExcIndexRange(from_level, 0, matrices.size()+1));
* Implementation of abstract
* function in #MGTranferBase#.
*/
- virtual void restrict (const unsigned int from_level,
- Vector<double> &dst,
- const Vector<double> &src) const;
+ virtual void restrict_and_add (const unsigned int from_level,
+ Vector<double> &dst,
+ const Vector<double> &src) const;
/**
* Exception.