--- /dev/null
+// $Id$
+
+#include "../tests.h"
+#include <deal.II/base/exceptions.h>
+#include <iostream>
+#include <deal.II/lac/petsc_matrix_free.h>
+#include <deal.II/lac/vector.h>
+
+
+/**
+ * Finite difference PETSc matrix-free object on uniform grid.
+ * Generator for simple 5-point discretization of Laplace problem.
+ */
+
+class PetscFDMatrix : public dealii::PETScWrappers::MatrixFree
+{
+ public:
+ /**
+ * Constructor specifying grid resolution.
+ */
+ PetscFDMatrix(unsigned int size, unsigned int dim);
+
+ /**
+ * Matrix-vector multiplication:
+ * let <i>dst = M*src</i> with
+ * <i>M</i> being this matrix.
+ *
+ * Source and destination must
+ * not be the same vector.
+ */
+ void vmult (dealii::PETScWrappers::VectorBase &dst,
+ const dealii::PETScWrappers::VectorBase &src) const;
+
+ /**
+ * Matrix-vector multiplication: let
+ * <i>dst = M<sup>T</sup>*src</i> with
+ * <i>M</i> being this matrix. This
+ * function does the same as @p vmult()
+ * but takes the transposed matrix.
+ *
+ * Source and destination must
+ * not be the same vector.
+ */
+ void Tvmult (dealii::PETScWrappers::VectorBase &dst,
+ const dealii::PETScWrappers::VectorBase &src) const;
+
+ /**
+ * Adding Matrix-vector
+ * multiplication. Add
+ * <i>M*src</i> on <i>dst</i>
+ * with <i>M</i> being this
+ * matrix.
+ *
+ * Source and destination must
+ * not be the same vector.
+ */
+ void vmult_add (dealii::PETScWrappers::VectorBase &dst,
+ const dealii::PETScWrappers::VectorBase &src) const;
+
+ /**
+ * Adding Matrix-vector
+ * multiplication. Add
+ * <i>M<sup>T</sup>*src</i> to
+ * <i>dst</i> with <i>M</i> being
+ * this matrix. This function
+ * does the same as @p vmult_add()
+ * but takes the transposed
+ * matrix.
+ *
+ * Source and destination must
+ * not be the same vector.
+ */
+ void Tvmult_add (dealii::PETScWrappers::VectorBase &dst,
+ const dealii::PETScWrappers::VectorBase &src) const;
+
+private:
+ /**
+ * Number of gridpoints in x-direction.
+ */
+ unsigned int nx;
+
+ /**
+ * Number of gridpoints in y-direction.
+ */
+ unsigned int ny;
+};
+
+
+// --------------- inline and template functions -----------------
+
+inline
+PetscFDMatrix::PetscFDMatrix(unsigned int size, unsigned int dim)
+ : PETScWrappers::MatrixFree (dim, dim, dim, dim),
+ nx (size), ny (size)
+{}
+
+
+
+inline
+void
+PetscFDMatrix::vmult_add (dealii::PETScWrappers::VectorBase &dst,
+ const dealii::PETScWrappers::VectorBase &src) const
+{
+ for(unsigned int i=0;i<=ny-2;i++)
+ {
+ for(unsigned int j=0;j<=nx-2; j++)
+ {
+ // Number of the row to be entered
+ unsigned int row = j+(nx-1)*i;
+
+ dst(row) += 4. * src(row); // A.set(row, row, 4.);
+
+ if (j>0)
+ {
+ dst(row-1) += -1. * src(row); // A.set(row-1, row, -1.);
+ dst(row) += -1. * src(row-1); // A.set(row, row-1, -1.);
+ }
+ if (i>0)
+ {
+ dst(row-(nx-1)) += -1. * src(row); // A.set(row-(nx-1), row, -1.);
+ dst(row) += -1. * src(row-(nx-1)); // A.set(row, row-(nx-1), -1.);
+ }
+ }
+ }
+}
+
+
+
+inline
+void
+PetscFDMatrix::vmult (dealii::PETScWrappers::VectorBase &dst,
+ const dealii::PETScWrappers::VectorBase &src) const
+{
+ dst = 0;
+ vmult_add (dst, src);
+}
+
+
+
+inline
+void
+PetscFDMatrix::Tvmult (dealii::PETScWrappers::VectorBase &dst,
+ const dealii::PETScWrappers::VectorBase &src) const
+{
+ dst = 0;
+ vmult_add (dst, src);
+}
+
+
+
+inline
+void
+PetscFDMatrix::Tvmult_add (dealii::PETScWrappers::VectorBase &dst,
+ const dealii::PETScWrappers::VectorBase &src) const
+{
+ vmult_add (dst, src);
+}
+
+
--- /dev/null
+//---------------------------- petsc_solver_03_mf.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2004, 2005, 2010 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- petsc_solver_03_mf.cc ---------------------------
+
+// test the PETSc CG solver with PETSc MatrixFree class
+
+
+#include "../tests.h"
+#include "../lac/petsc_mf_testmatrix.h"
+#include <cmath>
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/petsc_sparse_matrix.h>
+#include <deal.II/lac/petsc_vector.h>
+#include <deal.II/lac/petsc_solver.h>
+#include <deal.II/lac/petsc_precondition.h>
+#include <deal.II/lac/vector_memory.h>
+#include <typeinfo>
+
+template<class SOLVER, class MATRIX, class VECTOR, class PRECONDITION>
+void
+check_solve( SOLVER& solver, const MATRIX& A,
+ VECTOR& u, VECTOR& f, const PRECONDITION& P)
+{
+ deallog << "Solver type: " << typeid(solver).name() << std::endl;
+
+ u = 0.;
+ f = 1.;
+ try
+ {
+ solver.solve(A,u,f,P);
+ }
+ catch (std::exception& e)
+ {
+ std::cout << e.what() << std::endl;
+ deallog << e.what() << std::endl;
+ abort ();
+ }
+
+ deallog << "Solver stopped after " << solver.control().last_step()
+ << " iterations" << std::endl;
+}
+
+
+int main(int argc, char **argv)
+{
+ std::ofstream logfile("solver_03_mf/output");
+ deallog.attach(logfile);
+ deallog << std::setprecision(4);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ PetscInitialize(&argc,&argv,0,0);
+ {
+ SolverControl control(100, 1.e-3);
+
+ const unsigned int size = 32;
+ unsigned int dim = (size-1)*(size-1);
+
+ deallog << "Size " << size << " Unknowns " << dim << std::endl;
+
+ PetscFDMatrix A(size, dim);
+
+ PETScWrappers::Vector f(dim);
+ PETScWrappers::Vector u(dim);
+ f = 1.;
+ A.compress ();
+ f.compress ();
+ u.compress ();
+
+ PETScWrappers::SolverCG solver(control);
+ PETScWrappers::PreconditionNone preconditioner(A);
+ check_solve (solver, A,u,f, preconditioner);
+ }
+ PetscFinalize ();
+}
+
projects / dealii.git / commitdiff