double u = fe.shape_value(j,k);
elmat(i,j) += fe.JxW(k) *
- (1.*u*v + du*dv)
+ (0.1*u*v + du*dv)
;
}
}
double u = fe.shape_value(j,k);
elmat(i,j) += fe.JxW(k) *
- (1.*u*v + du*dv)
+ (0.1*u*v + du*dv)
;
}
}
public MGSmootherBase
{
private:
- SmartPointer<MGMatrix<SparseMatrix<float> > >matrices;
+ SmartPointer<MGMatrix<SparseMatrix<double> > >matrices;
public:
- MGSmootherLAC(MGMatrix<SparseMatrix<float> >&);
+ MGSmootherLAC(MGMatrix<SparseMatrix<double> >&);
virtual void smooth (const unsigned int level,
- Vector<float> &u,
- const Vector<float> &rhs) const;
+ Vector<double> &u,
+ const Vector<double> &rhs) const;
};
u = 0.;
vector<SparseMatrixStruct> mgstruct(tr.n_levels());
- MGMatrix<SparseMatrix<float> > mgA(0,tr.n_levels()-1);
+ MGMatrix<SparseMatrix<double> > mgA(0,tr.n_levels()-1);
for (unsigned int i=0;i<tr.n_levels();++i)
{
mgstruct[i].reinit(mgdof.n_dofs(i), mgdof.n_dofs(i),
equation.build_mgmatrix(mgA, mgdof, quadrature);
SolverControl cgcontrol(10,0., false, false);
- PrimitiveVectorMemory<Vector<float> > cgmem;
- SolverCG<SparseMatrix<float>, Vector<float> > cgsolver(cgcontrol, cgmem);
- PreconditionIdentity<Vector<float> > cgprec;
- MGCoarseGridLACIteration<SolverCG<SparseMatrix<float>, Vector<float> >,
- SparseMatrix<float>, PreconditionIdentity<Vector<float> > >
+ PrimitiveVectorMemory<Vector<double> > cgmem;
+ SolverCG<SparseMatrix<double>, Vector<double> > cgsolver(cgcontrol, cgmem);
+ PreconditionIdentity<Vector<double> > cgprec;
+ MGCoarseGridLACIteration<SolverCG<SparseMatrix<double>, Vector<double> >,
+ SparseMatrix<double>, PreconditionIdentity<Vector<double> > >
coarse(cgsolver, mgA[0], cgprec);
MGSmootherLAC smoother(mgA);
return 1.;
}
-MGSmootherLAC::MGSmootherLAC(MGMatrix<SparseMatrix<float> >& matrix)
+MGSmootherLAC::MGSmootherLAC(MGMatrix<SparseMatrix<double> >& matrix)
:
matrices(&matrix)
{}
void
MGSmootherLAC::smooth (const unsigned int level,
- Vector<float> &u,
- const Vector<float> &rhs) const
+ Vector<double> &u,
+ const Vector<double> &rhs) const
{
- SolverControl control(1,1.e-300,false,false);
- PrimitiveVectorMemory<Vector<float> > mem;
- SolverRichardson<SparseMatrix<float> , Vector<float> > rich(control, mem);
- PreconditionRelaxation<SparseMatrix<float> , Vector<float> >
- prec((*matrices)[level], &SparseMatrix<float> ::template precondition_SSOR<float>, 1.);
+ 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> >
+ prec((*matrices)[level], &SparseMatrix<double> ::template precondition_SSOR<double>, 1.);
rich.solve((*matrices)[level], u, rhs, prec);
}
# Automatic generation of dependencies
############################################################
-all-cc-files = $(solver-cc-files) $(mgbase-cc-files) $(mg-cc-files) $(benchmark-cc-files)
+all-cc-files = $(wildcard *.cc)
Make.depend: $(all-cc-files)
@echo =====Dependecies=== Make.depend
#include <lac/mgbase.h>
#define TYPE long double
-#define ACCURACY 1.e-18
+#define ACCURACY 1.e-20
template<class VECTOR>
void print_vector(ostream& s, const VECTOR& v)
/**
* Pointer to the level matrices.
*/
- SmartPointer<MGMatrix<SparseMatrix<float> > >matrices;
+ SmartPointer<MGMatrix<SparseMatrix<double> > >matrices;
public:
- FDMG(unsigned int maxlevel, MGMatrix<SparseMatrix<float> >& matrices,
+ FDMG(unsigned int maxlevel, MGMatrix<SparseMatrix<double> >& matrices,
FDMGTransfer& transfer);
- virtual void level_residual (unsigned int level,
- Vector<float>& dst,
- const Vector<float>& src,
- const Vector<float>& rhs);
+ virtual void level_vmult (unsigned int level,
+ Vector<double>& dst,
+ const Vector<double>& src,
+ const Vector<double>& rhs);
void copy_to_mg(const Vector<TYPE>& rhs);
void copy_from_mg(Vector<TYPE>& lsg);
public MGSmootherBase
{
private:
- SmartPointer<MGMatrix<SparseMatrix<float> > >matrices;
+ SmartPointer<MGMatrix<SparseMatrix<double> > >matrices;
public:
- MGSmootherLAC(MGMatrix<SparseMatrix<float> >&);
+ MGSmootherLAC(MGMatrix<SparseMatrix<double> >&);
virtual void smooth (const unsigned int level,
- Vector<float> &u,
- const Vector<float> &rhs) const;
+ Vector<double> &u,
+ const Vector<double> &rhs) const;
};
-typedef MGCoarseGridLACIteration<SolverCG<SparseMatrix<float> , Vector<float> >,
-SparseMatrix<float>, /*PreconditionRelaxation<SparseMatrix<float> ,*/
- PreconditionIdentity<Vector<float> > >
+typedef MGCoarseGridLACIteration<SolverCG<SparseMatrix<double> , Vector<double> >,
+SparseMatrix<double>, /*PreconditionRelaxation<SparseMatrix<double> ,*/
+ PreconditionIdentity<Vector<double> > >
Coarse;
deallog.attach(logfile);
PrimitiveVectorMemory<Vector<TYPE> > mem;
- SolverControl control(100, ACCURACY, true);
+ SolverControl control(10000, ACCURACY);
const unsigned int base = 3;
const unsigned int maxlevel = 8;
FDMGTransfer transfer(maxsize, maxsize, maxlevel);
// coarse grid solver
- PrimitiveVectorMemory<Vector<float> > cgmem;
+ PrimitiveVectorMemory<Vector<double> > cgmem;
ReductionControl cgcontrol(100, 1.e-30, 1.e-2, false, false);
- SolverCG<SparseMatrix<float> , Vector<float> > cgcg(cgcontrol,cgmem);
+ SolverCG<SparseMatrix<double> , Vector<double> > cgcg(cgcontrol,cgmem);
for (unsigned int level = 0; level <= maxlevel; ++level)
// Make matrix
vector<SparseMatrixStruct > structure(maxlevel+1);
- MGMatrix<SparseMatrix<float> > A(minlevel,maxlevel);
+ MGMatrix<SparseMatrix<double> > A(minlevel,maxlevel);
FDMatrix testproblem(size, size);
FDMG multigrid(level, A, transfer);
-// PreconditionRelaxation<SparseMatrix<float> , Vector<float> >
- PreconditionIdentity<Vector<float> >
- cgprec;//(A[minlevel], &SparseMatrix<float> ::template precondition_SSOR<float>, 1.2);
+// PreconditionRelaxation<SparseMatrix<double> , Vector<double> >
+ PreconditionIdentity<Vector<double> >
+ cgprec;//(A[minlevel], &SparseMatrix<double> ::template precondition_SSOR<double>, 1.2);
Coarse coarsegrid(cgcg, A[minlevel], cgprec);
PreconditionMG<FDMG, Vector<TYPE> >
precondition(multigrid, smoother, smoother, coarsegrid);
-// SolverRichardson<SparseMatrix<float> , Vector<TYPE> > solver(control, mem);
- SolverCG<SparseMatrix<float> , Vector<TYPE> > solver(control, mem);
+// SolverRichardson<SparseMatrix<double> , Vector<TYPE> > solver(control, mem);
+ SolverCG<SparseMatrix<double> , Vector<TYPE> > solver(control, mem);
Vector<TYPE> u(dim);
Vector<TYPE> f(dim);
}
}
-FDMG::FDMG(unsigned int maxlevel, MGMatrix<SparseMatrix<float> >& matrices,
+FDMG::FDMG(unsigned int maxlevel, MGMatrix<SparseMatrix<double> >& matrices,
FDMGTransfer& transfer)
:
MGBase(transfer, 0, maxlevel),
void
-FDMG::level_residual (unsigned int level,
- Vector<float>& dst,
- const Vector<float>& src,
- const Vector<float>& rhs)
+FDMG::level_vmult (unsigned int level,
+ Vector<double>& dst,
+ const Vector<double>& src,
+ const Vector<double>&)
{
- (*matrices)[level].residual(dst, src, rhs);
+ (*matrices)[level].vmult(dst, src);
}
void
}
-MGSmootherLAC::MGSmootherLAC(MGMatrix<SparseMatrix<float> >& matrix)
+MGSmootherLAC::MGSmootherLAC(MGMatrix<SparseMatrix<double> >& matrix)
:
matrices(&matrix)
{}
void
MGSmootherLAC::smooth (const unsigned int level,
- Vector<float> &u,
- const Vector<float> &rhs) const
+ Vector<double> &u,
+ const Vector<double> &rhs) const
{
SolverControl control(1,1.e-300,false,false);
- PrimitiveVectorMemory<Vector<float> > mem;
- SolverRichardson<SparseMatrix<float> , Vector<float> > rich(control, mem);
- PreconditionRelaxation<SparseMatrix<float> , Vector<float> >
- prec((*matrices)[level], &SparseMatrix<float> ::template precondition_SSOR<float>, 1.);
+ PrimitiveVectorMemory<Vector<double> > mem;
+ SolverRichardson<SparseMatrix<double> , Vector<double> > rich(control, mem);
+ PreconditionRelaxation<SparseMatrix<double> , Vector<double> >
+ prec((*matrices)[level], &SparseMatrix<double> ::template precondition_SSOR<double>, 1.);
rich.solve((*matrices)[level], u, rhs, prec);
}
public MGTransferBase
{
virtual void prolongate(unsigned l,
- Vector<float>& dst,
- const Vector<float>& src) const;
+ Vector<double>& dst,
+ const Vector<double>& src) const;
virtual void restrict(unsigned l,
- Vector<float>& dst,
- const Vector<float>& src) const;
+ Vector<double>& dst,
+ const Vector<double>& src) const;
friend class MGBase;
};
{
public:
virtual void smooth (const unsigned int level,
- Vector<float> &u,
- const Vector<float>& rhs) const;
+ Vector<double> &u,
+ const Vector<double>& rhs) const;
};
class MGCoarseGridTest
{
public:
virtual void operator () (unsigned int level,
- Vector<float>& dst,
- const Vector<float>& src) const ;
+ Vector<double>& dst,
+ const Vector<double>& src) const ;
};
public:
MGTest(TransferTest& t, SmoothTest& sm, MGCoarseGridTest& cs)
- : MGBase(t, 9, 3),
+ : MGBase(t, 3, 9),
smoother(sm),
solver(cs)
{
}
}
- virtual void level_residual(unsigned level,
- Vector<float>& dst,
- const Vector<float>& src,
- const Vector<float>& rhs);
+ virtual void level_vmult(unsigned level,
+ Vector<double>& dst,
+ const Vector<double>& src,
+ const Vector<double>& rhs);
void doit()
{
void
TransferTest::prolongate(unsigned l,
- Vector<float>&,
- const Vector<float>&) const
+ Vector<double>&,
+ const Vector<double>&) const
{
cout << "Prolongating " << l-1 << " to " << l << endl;
}
void
TransferTest::restrict(unsigned l,
- Vector<float>&,
- const Vector<float>&) const
+ Vector<double>&,
+ const Vector<double>&) const
{
cout << "Restricting " << l << " to " << l-1 << endl;
}
void
SmoothTest::smooth (const unsigned int level,
- Vector<float> &,
- const Vector<float>&) const
+ Vector<double> &,
+ const Vector<double>&) const
{
cout << "Smoothing on " << level << endl;
}
void
-MGTest::level_residual(unsigned l,
- Vector<float>&,
- const Vector<float>&,
- const Vector<float>&)
+MGTest::level_vmult(unsigned l,
+ Vector<double>&,
+ const Vector<double>&,
+ const Vector<double>&)
{
cout << "Residual on " << l << endl;
}
void
MGCoarseGridTest::operator() (unsigned int level,
- Vector<float>&,
- const Vector<float>&) const
+ Vector<double>&,
+ const Vector<double>&) const
{
cout << "Solving on " << level << endl;
}
#include <cmath>
#include <fstream>
+#include <iomanip>
#include "testmatrix.h"
#include <base/logstream.h>
#include <lac/sparsematrix.h>
deallog.attach(logfile);
PrimitiveVectorMemory<Vector<double> > mem;
- SolverControl control(100, 1.e-5, true);
+ SolverControl control(50, 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);
- for (unsigned int size=10; size <= 10; size *= 10)
+ for (unsigned int size=4; size <= 40; size *= 3)
{
deallog << "Size " << size << endl;
Vector<double> f(dim);
Vector<double> u(dim);
+ Vector<double> res(dim);
+
+ f = 1.;
+ u = 1.;
+
+ A.residual(res,u,f);
+ A.SOR(res);
+ res.add(1.,u);
+ A.SOR_step(u,f);
+ res.add(-1.,u);
+
+ deallog << "SOR-diff:" << res*res << endl;
deallog.push("no");
void
FDMGTransfer::build_matrix(unsigned int nx, unsigned int ny,
- SparseMatrixStruct& structure, SparseMatrix<float>& matrix)
+ SparseMatrixStruct& structure, SparseMatrix<double>& matrix)
{
structure.reinit((nx-1)*(ny-1),(2*nx-1)*(2*ny-1),9);
void
FDMGTransfer::prolongate (const unsigned int to_level,
- Vector<float> &dst,
- const Vector<float> &src) const
+ Vector<double> &dst,
+ const Vector<double> &src) const
{
Assert((to_level>0) && (to_level<=matrices.size()),
ExcIndexRange(to_level, 0, matrices.size()+1));
void
FDMGTransfer::restrict (const unsigned int from_level,
- Vector<float> &dst,
- const Vector<float> &src) const
+ Vector<double> &dst,
+ const Vector<double> &src) const
{
Assert((from_level>0) && (from_level<=matrices.size()),
ExcIndexRange(from_level, 0, matrices.size()+1));
* function in #MGTranferBase#.
*/
virtual void prolongate (const unsigned int to_level,
- Vector<float> &dst,
- const Vector<float> &src) const;
+ Vector<double> &dst,
+ const Vector<double> &src) const;
/**
* Implementation of abstract
* function in #MGTranferBase#.
*/
virtual void restrict (const unsigned int from_level,
- Vector<float> &dst,
- const Vector<float> &src) const;
+ Vector<double> &dst,
+ const Vector<double> &src) const;
/**
* Exception.
/**
* Prolongation matrices.
*/
- vector<SparseMatrix<float> > matrices;
+ vector<SparseMatrix<double> > matrices;
/**
* Matrix generator.
* fine to coarse (#vmult#).
*/
void build_matrix(unsigned int nx, unsigned int ny,
- SparseMatrixStruct&, SparseMatrix<float>&);
+ SparseMatrixStruct&, SparseMatrix<double>&);
};