class MGMatrix
:
public Subscriptor,
- public Vector<MATRIX>
+ public vector<MATRIX>
{
public:
/**
* preconditioning method for later
* use.
*/
- MGCoarseGridLACIteration(const SOLVER&, const MATRIX&, const PRECOND&);
+ MGCoarseGridLACIteration( SOLVER&, const MATRIX&, const PRECOND&);
/**
* Implementation of the abstract
/**
* Reference to the solver.
*/
- const SOLVER& solver;
+ SOLVER& solver;
/**
* Reference to the matrix.
*/
template<class SOLVER, class MATRIX, class PRECOND>
MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND>
-::MGCoarseGridLACIteration(const SOLVER& s, const MATRIX& m, const PRECOND& p)
+::MGCoarseGridLACIteration(SOLVER& s, const MATRIX& m, const PRECOND& p)
:
solver(s), matrix(m), precondition(p)
{}
template<class SOLVER, class MATRIX, class PRECOND>
void
MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND>::operator()
-(unsigned int level, Vector<float>& dst, const Vector<float>& src) const
+(unsigned int, Vector<float>& dst, const Vector<float>& src) const
{
solver.solve(matrix, dst, src, precondition);
}
post_smooth.smooth(level, s[level], d[level]);
}
+MGTransferBase::~MGTransferBase()
+{}
+
+MGSmootherBase::~MGSmootherBase()
+{}
class_declaration: class_head
| class_head
{
+ cout << "@@@" << endl << setw(OUTW) << "@Class-Definition:" << $1 << endl
+ << setw(OUTW) << "@In-Class:" << class_stack.top() << endl;
class_name = class_stack.top() + string("::") + class_name;
class_stack.push(class_name);
access_stack.push(string("private"));
{
class_stack.pop();
access_stack.pop();
- cout << "@@@" << endl << setw(OUTW) << "@Class-Definition:" << $1 << endl
- << setw(OUTW) << "@In-Class:" << class_stack.top() << endl;
}
;
class MGMatrix
:
public Subscriptor,
- public Vector<MATRIX>
+ public vector<MATRIX>
{
public:
/**
* preconditioning method for later
* use.
*/
- MGCoarseGridLACIteration(const SOLVER&, const MATRIX&, const PRECOND&);
+ MGCoarseGridLACIteration( SOLVER&, const MATRIX&, const PRECOND&);
/**
* Implementation of the abstract
/**
* Reference to the solver.
*/
- const SOLVER& solver;
+ SOLVER& solver;
/**
* Reference to the matrix.
*/
template<class SOLVER, class MATRIX, class PRECOND>
MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND>
-::MGCoarseGridLACIteration(const SOLVER& s, const MATRIX& m, const PRECOND& p)
+::MGCoarseGridLACIteration(SOLVER& s, const MATRIX& m, const PRECOND& p)
:
solver(s), matrix(m), precondition(p)
{}
template<class SOLVER, class MATRIX, class PRECOND>
void
MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND>::operator()
-(unsigned int level, Vector<float>& dst, const Vector<float>& src) const
+(unsigned int, Vector<float>& dst, const Vector<float>& src) const
{
solver.solve(matrix, dst, src, precondition);
}
post_smooth.smooth(level, s[level], d[level]);
}
+MGTransferBase::~MGTransferBase()
+{}
+
+MGSmootherBase::~MGSmootherBase()
+{}
--- /dev/null
+// $Id$
+
+// deal_II_libraries.g=-ldeal_II_2d.g
+// deal_II_libraries=-ldeal_II_2d
+
+#include <numerics/multigrid.h>
+#include <numerics/mg_smoother.h>
+#include <grid/tria.h>
+#include <grid/mg_dof.h>
+#include <fe/fe_lib.lagrange.h>
+
+main()
+{
+ Triangulation<2> tr;
+ MGDoFHandler<2> dof(&tr);
+ FELinear<2> fe;
+
+ tr.create_hypercube(-1.,1.);
+ tr.refine_global(3);
+ dof.distribute_dofs(fe);
+}
+
flags = $(CXXFLAGS.o) $(CXXFLAGS)
endif
-all: solver
+all: solver mgbase mg
############################################################
# Typical block for building a running program
solver: $(solver-o-files) $(libraries)
$(CXX) $(flags) -o $@ $^
+
+############################################################
+
+
+mgbase-cc-files = mgbase.cc
+
+ifeq ($(debug-mode),on)
+mgbase-o-files = $(mgbase-cc-files:.cc=.go)
+else
+mgbase-o-files = $(mgbase-cc-files:.cc=.o)
+endif
+
+mgbase: $(mgbase-o-files) $(libraries)
+ $(CXX) $(flags) -o $@ $^
+
+
+
+############################################################
+
+
+mg-cc-files = mg.cc testmatrix.cc
+
+ifeq ($(debug-mode),on)
+mg-o-files = $(mg-cc-files:.cc=.go)
+else
+mg-o-files = $(mg-cc-files:.cc=.o)
+endif
+
+mg: $(mg-o-files) $(libraries)
+ $(CXX) $(flags) -o $@ $^
+
+
+
############################################################
# Continue with other targets if needed
############################################################
# Automatic generation of dependencies
############################################################
-all-cc-files = $(solver-cc-files) # $(target2-cc-files) ...
+all-cc-files = $(solver-cc-files) $(mgbase-cc-files) $(mg-cc-files) # $(target2-cc-files) ...
Make.depend: $(all-cc-files)
@echo =====Dependecies=== Make.depend
--- /dev/null
+// $Id$
+
+// Test program for multigrid with test matrices (fd)
+
+#include <cmath>
+#include <fstream>
+#include "testmatrix.h"
+#include <base/logstream.h>
+#include <lac/sparsematrix.h>
+#include <lac/vector.h>
+#include <lac/vector_memory.h>
+#include <lac/solver_control.h>
+#include <lac/solver_cg.h>
+#include <lac/solver_gmres.h>
+#include <lac/solver_bicgstab.h>
+#include <lac/solver_richardson.h>
+#include <lac/precondition.h>
+#include <lac/mgbase.h>
+
+class FDMG
+ :
+ public MGBase
+{
+ private:
+ SmartPointer<MGMatrix<SparseMatrix<float> > >matrices;
+ public:
+ FDMG(unsigned int maxlevel, MGMatrix<SparseMatrix<float> >& matrices,
+ FDMGTransfer& transfer)
+ :
+ MGBase(transfer, maxlevel, 0),
+ matrices(&matrices)
+ {}
+
+
+ virtual void level_residual (unsigned int level,
+ Vector<float>& dst,
+ const Vector<float>& src,
+ const Vector<float>& rhs);
+};
+
+class MGSmootherLAC
+ :
+ MGSmootherBase
+{
+ private:
+ SmartPointer<MGMatrix<SparseMatrix<float> > >matrices;
+ public:
+ MGSmootherLAC(MGMatrix<SparseMatrix<float> >&);
+
+ virtual void smooth (const unsigned int level,
+ Vector<float> &u,
+ const Vector<float> &rhs) const;
+
+};
+
+typedef MGCoarseGridLACIteration<SolverCG<SparseMatrix<float> , Vector<float> >,
+SparseMatrix<float>, PreconditionRelaxation<SparseMatrix<float> , Vector<float> > >
+Coarse;
+
+class MGPrecondition
+{
+ private:
+ FDMG& mg;
+ MGSmootherLAC& smooth;
+ Coarse& coarse;
+ public:
+ MGPrecondition(FDMG& m, MGSmootherLAC& s, Coarse& c)
+ :
+ mg(m), smooth(s), coarse(c)
+ {}
+
+
+ void operator () (Vector<double>& dst, const Vector<double>& src) const
+ {
+
+ }
+
+
+};
+
+
+main()
+{
+ ofstream logfile("mg.output");
+ deallog.attach(logfile);
+
+ PrimitiveVectorMemory<Vector<double> > mem;
+ SolverControl control(100, 1.e-5, true);
+
+ const unsigned int base = 3;
+ const unsigned int maxlevel = 6;
+
+ const unsigned int maxsize = base * (1 << maxlevel);
+
+ // grid transfer
+ FDMGTransfer transfer(maxsize, maxsize, maxlevel);
+
+ // coarse grid solver
+ PrimitiveVectorMemory<Vector<float> > cgmem;
+ SolverControl cgcontrol(100, 1.e-12);
+ SolverCG<SparseMatrix<float> , Vector<float> > cgcg(cgcontrol,cgmem);
+
+ vector<SparseMatrixStruct > structure(maxlevel+1);
+ MGMatrix<SparseMatrix<float> > A(maxlevel+1);
+
+ for (unsigned int level = 0; level <= maxlevel; ++level)
+ {
+ const unsigned int size = base * (1 << level);
+
+ const unsigned int dim = (size+1)*(size+1);
+
+ // Make matrix
+ FDMatrix testproblem(size, size);
+ structure[level].reinit(dim, dim, 5);
+ testproblem.build_structure(structure[level]);
+ structure[level].compress();
+ A[level].reinit(structure[level]);
+ testproblem.laplacian(A[level]);
+
+ FDMG multigrid(level, A, transfer);
+
+ PreconditionRelaxation<SparseMatrix<float> , Vector<float> >
+ cgprec(A[0], &SparseMatrix<float> ::template precondition_SSOR<float>, 1.2);
+
+ Coarse coarsegrid(cgcg, A[0], cgprec);
+
+ MGSmootherLAC smoother(A);
+
+// MGPrecondition prec(multigrid);
+ }
+}
+
+void
+FDMG::level_residual (unsigned int level,
+ Vector<float>& dst,
+ const Vector<float>& src,
+ const Vector<float>& rhs)
+{
+ (*matrices)[level].residual(dst, src, rhs);
+}
+
+MGSmootherLAC::MGSmootherLAC(MGMatrix<SparseMatrix<float> >& matrix)
+ :
+ matrices(&matrix)
+{}
+
+
+void
+MGSmootherLAC::smooth (const unsigned int level,
+ Vector<float> &u,
+ const Vector<float> &rhs) const
+{
+ SolverControl control(0.,1);
+ 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.2);
+
+ rich.solve((*matrices)[level], u, rhs, prec);
+}
// deal_II_libraries.g=-ldeal_II_2d.g
// deal_II_libraries=-ldeal_II_2d
-#include <numerics/multigrid.h>
-#include <numerics/mg_smoother.h>
+#include <lac/mgbase.h>
class TransferTest
:
virtual void restrict(unsigned l,
Vector<float>& dst,
const Vector<float>& src) const;
- friend class MultiGridBase;
+ friend class MGBase;
};
class SmoothTest
const Vector<float>& rhs) const;
};
+class MGCoarseGridTest
+ :
+ public MGCoarseGridSolver
+{
+ public:
+ virtual void operator () (unsigned int level,
+ Vector<float>& dst,
+ const Vector<float>& src) const ;
+};
+
class MGTest
:
- public MultiGridBase
+ public MGBase
{
MGSmootherBase& smoother;
-
+ MGCoarseGridSolver& solver;
+
public:
- MGTest(TransferTest& t, SmoothTest& s)
- : MultiGridBase(t, 9, 3),
- smoother(s)
+ MGTest(TransferTest& t, SmoothTest& s, MGCoarseGridTest& c)
+ : MGBase(t, 9, 3),
+ smoother(s),
+ solver(c)
{}
virtual void level_residual(unsigned level,
void doit()
{
- level_mgstep(9, smoother, smoother);
+ level_mgstep(9, smoother, smoother, solver);
}
};
{
TransferTest tr;
SmoothTest s;
- MGTest mg(tr, s);
+ MGCoarseGridTest c;
+
+ MGTest mg(tr, s, c);
mg.doit();
}
Vector<float> &,
const Vector<float>&) const
{
- cout << "Smoothing on" << level << endl;
+ cout << "Smoothing on " << level << endl;
}
void
{
cout << "Residual on " << l << endl;
}
+
+void
+MGCoarseGridTest::operator() (unsigned int level,
+ Vector<float>&,
+ const Vector<float>&) const
+{
+ cout << "Solving on " << level << endl;
+}
// Test program for linear solvers.
#include <cmath>
+#include <fstream>
#include "testmatrix.h"
#include <base/logstream.h>
#include <lac/sparsematrix.h>
#include <lac/solver_richardson.h>
#include <lac/precondition.h>
-#define MATRIX SparseMatrix<float>
-#define VECTOR Vector<double>
+template<class SOLVER, class MATRIX, class VECTOR, class PRECONDITION>
+void
+check_method( SOLVER& solver, const MATRIX& A,
+ VECTOR& u, VECTOR& f, const PRECONDITION& P)
+{
+ u = 0.;
+ f = 1.;
+ solver.solve(A,u,f,P);
+}
main()
{
- PrimitiveVectorMemory<VECTOR > mem;
+ ofstream logfile("solver.output");
+ deallog.attach(logfile);
+
+ PrimitiveVectorMemory<Vector<double> > mem;
SolverControl control(100, 1.e-5, true);
- SolverCG<MATRIX, VECTOR > cg(control, mem);
- SolverGMRES<MATRIX, VECTOR > gmres(control, mem,20);
- SolverBicgstab<MATRIX, VECTOR > bicgstab(control, mem);
- SolverRichardson<MATRIX, VECTOR > rich(control, mem);
+ 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)
{
SparseMatrixStruct structure(dim, dim, 5);
testproblem.build_structure(structure);
structure.compress();
- MATRIX A(structure);
+ SparseMatrix<float> A(structure);
testproblem.laplacian(A);
- PreconditionIdentity<VECTOR >
+ PreconditionIdentity<Vector<double> >
prec_no;
- PreconditionRelaxation<MATRIX, VECTOR>
- prec_ssor(A, &MATRIX::template precondition_SSOR<double>, 1.2);
+ 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);
- VECTOR f(dim);
- VECTOR u(dim);
+ Vector<double> f(dim);
+ Vector<double> u(dim);
deallog.push("no");
- f = 1.;
- u = 0.;
- cg.solve(A,u,f,prec_no);
-
- f = 1.;
- u = 0.;
- bicgstab.solve(A,u,f,prec_no);
-
- f = 1.;
- u = 0.;
- gmres.solve(A,u,f,prec_no);
+ check_method(cg,A,u,f,prec_no);
+ check_method(bicgstab,A,u,f,prec_no);
+ check_method(gmres,A,u,f,prec_no);
deallog.pop();
deallog.push("ssor");
- f = 1.;
- u = 0.;
- rich.solve(A,u,f,prec_ssor);
-
- f = 1.;
- u = 0.;
- cg.solve(A,u,f,prec_ssor);
+ check_method(rich,A,u,f,prec_ssor);
+ check_method(cg,A,u,f,prec_ssor);
+ check_method(bicgstab,A,u,f,prec_ssor);
+ check_method(gmres,A,u,f,prec_ssor);
- f = 1.;
- u = 0.;
- bicgstab.solve(A,u,f,prec_ssor);
+ deallog.pop();
+ deallog.push("sor");
- f = 1.;
- u = 0.;
- gmres.solve(A,u,f,prec_ssor);
+ check_method(rich,A,u,f,prec_sor);
+ check_method(cg,A,u,f,prec_sor);
+ check_method(bicgstab,A,u,f,prec_sor);
+ check_method(gmres,A,u,f,prec_sor);
deallog.pop();
}
}
}
+FDMGTransfer::FDMGTransfer(unsigned int nx, unsigned int ny,
+ unsigned int nlevels)
+ :
+ structures(nlevels-1), matrices(nlevels-1)
+{
+ unsigned int power = 1 << nlevels;
+
+ Assert ((nx%power)==0, ExcDivide(nx,power));
+ Assert ((ny%power)==0, ExcDivide(ny,power));
+
+ nx /= power;
+ ny /= power;
+
+ for (unsigned int level = 0; level < nlevels-1; ++ level)
+ {
+ build_matrix(nx,ny,structures[level],matrices[level]);
+ nx *= 2;
+ ny *= 2;
+ }
+}
+
+void
+FDMGTransfer::build_matrix(unsigned int nx, unsigned int ny,
+ SparseMatrixStruct& structure, SparseMatrix<float>& matrix)
+{
+ structure.reinit((nx+1)*(ny+1),(2*nx+1)*(2*ny+1),9);
+
+ // Traverse all points of coarse grid
+ for (unsigned int i=0 ; i<=ny; ++i)
+ for (unsigned int j=0 ; j<=nx; ++j)
+ {
+ // coarse grid point number
+ unsigned int ncoarse =j+(nx+1)*i;
+ // same point on fine grid
+ unsigned int nfine = 2*j+(4*nx+2)*i;
+
+ structure.add(ncoarse,nfine);
+ // left
+ if (j>0)
+ {
+ structure.add(ncoarse,nfine-1);
+ // lower left
+ if (i>0)
+ structure.add(ncoarse,nfine-(2*nx+1)-1);
+ // upper left
+ if (i<ny)
+ structure.add(ncoarse,nfine+(2*nx+1)-1);
+ }
+ // right
+ if (j<nx)
+ {
+ structure.add(ncoarse, nfine+1);
+ // lower right
+ if (i>0)
+ structure.add(ncoarse,nfine-(2*nx+1)+1);
+ // upper right
+ if (i<ny)
+ structure.add(ncoarse,nfine+(2*nx+1)+1);
+ }
+
+ // lower
+ if (i>0)
+ structure.add(ncoarse,nfine-(2*nx+1));
+ // upper
+ if (i<ny)
+ structure.add(ncoarse,nfine+(2*nx+1));
+ }
+
+ structure.compress();
+ matrix.reinit(structure);
+
+ for (unsigned int i=0 ; i<=ny; ++i)
+ for (unsigned int j=0 ; j<=nx; ++j)
+ {
+ // coarse grid point number
+ unsigned int ncoarse =j+(nx+1)*i;
+ // same point on fine grid
+ unsigned int nfine = 2*j+(4*nx+2)*i;
+
+ matrix.set(ncoarse,nfine,1.);
+ // left
+ if (j>0)
+ {
+ matrix.set(ncoarse,nfine-1,.5);
+ // lower left
+ if (i>0)
+ matrix.set(ncoarse,nfine-(2*nx+1)-1,.25);
+ // upper left
+ if (i<ny)
+ matrix.set(ncoarse,nfine+(2*nx+1)-1,.25);
+ }
+ // right
+ if (j<nx)
+ {
+ matrix.set(ncoarse,nfine+1,.5);
+ // lower right
+ if (i>0)
+ matrix.set(ncoarse,nfine-(2*nx+1)+1,.25);
+ // upper right
+ if (i<ny)
+ matrix.set(ncoarse,nfine+(2*nx+1)+1,.25);
+ }
+
+ // lower
+ if (i>0)
+ matrix.set(ncoarse,nfine-(2*nx+1),.5);
+ // upper
+ if (i<ny)
+ matrix.set(ncoarse,nfine+(2*nx+1),.5);
+ }
+
+}
+
+void
+FDMGTransfer::prolongate (const unsigned int to_level,
+ Vector<float> &dst,
+ const Vector<float> &src) const
+{
+ Assert((to_level>0) && (to_level<=matrices.size()),
+ ExcIndexRange(to_level, 0, matrices.size()+1));
+
+ matrices[to_level-1].Tvmult(dst,src);
+}
+
+
+void
+FDMGTransfer::restrict (const unsigned int from_level,
+ Vector<float> &dst,
+ const Vector<float> &src) const
+{
+ Assert((from_level>0) && (from_level<=matrices.size()),
+ ExcIndexRange(from_level, 0, matrices.size()+1));
+
+ matrices[from_level-1].vmult(dst,src);
+}
+
+
template void FDMatrix::laplacian(SparseMatrix<double>&) const;
template void FDMatrix::laplacian(SparseMatrix<float>&) const;
// $Id$
#include <lac/forward-declarations.h>
+#include <base/exceptions.h>
+#include <lac/mgbase.h>
/**
* Finite difference matrix on uniform grid.
*/
unsigned int ny;
};
+
+/**
+ * Grid transfer for finite differences on uniform grid.
+ */
+
+class FDMGTransfer
+ :
+ public MGTransferBase
+{
+ public:
+ /**
+ * Constructor. Prepares grid
+ * transfer matrices for #nlevels#
+ * levels on an #nx# times #ny#
+ * grid.
+ */
+ FDMGTransfer(unsigned int nx, unsigned int ny,
+ unsigned int nlevels);
+
+ /**
+ * Implementation of abstract
+ * function in #MGTranferBase#.
+ */
+ virtual void prolongate (const unsigned int to_level,
+ Vector<float> &dst,
+ const Vector<float> &src) const;
+
+ /**
+ * Implementation of abstract
+ * function in #MGTranferBase#.
+ */
+ virtual void restrict (const unsigned int from_level,
+ Vector<float> &dst,
+ const Vector<float> &src) const;
+
+ /**
+ * Exception.
+ */
+ DeclException2(ExcDivide, unsigned int, unsigned int,
+ << "Cannot divide " << arg1 << " by " << arg2);
+
+ private:
+ /**
+ * Prolongation matrix structures.
+ */
+ vector<SparseMatrixStruct > structures;
+ /**
+ * Prolongation matrices.
+ */
+ vector<SparseMatrix<float> > matrices;
+
+ /**
+ * Matrix generator.
+ * The arguments #nx# and #ny#
+ * are the numbers on the COARSE level.
+ * Builds a transfer matrix from
+ * fine to coarse (#vmult#).
+ */
+ void build_matrix(unsigned int nx, unsigned int ny,
+ SparseMatrixStruct&, SparseMatrix<float>&);
+};