dolfin team mailing list archive
-
dolfin team
-
Mailing list archive
-
Message #18796
Re: [Branch ~dolfin-core/dolfin/main] Rev 4868: Improve re-use in PETScLUSolver.
This changeset broke the Poisson demo.
Solving linear variational problem
Reset matrix (B)
Applying boundary conditions to linear system.
Solving linear system of size 1089 x 1089 (PETSc LU solver, umfpack).
[0]PETSC ERROR: --------------------- Error Message
------------------------------------
[0]PETSC ERROR: Object is in wrong state!
[0]PETSC ERROR: Matrix must be set first!
The advection diffusion demo works fine.
Kristian
On 22 July 2010 19:08, <noreply@xxxxxxxxxxxxx> wrote:
> ------------------------------------------------------------
> revno: 4868
> committer: Garth N. Wells <gnw20@xxxxxxxxx>
> branch nick: dolfin-all
> timestamp: Thu 2010-07-22 19:05:26 +0100
> message:
> Improve re-use in PETScLUSolver.
> modified:
> demo/pde/advection-diffusion/cpp/main.cpp
> dolfin/la/PETScBaseMatrix.h
> dolfin/la/PETScKrylovSolver.cpp
> dolfin/la/PETScKrylovSolver.h
> dolfin/la/PETScLUSolver.cpp
> dolfin/la/PETScLUSolver.h
> dolfin/la/PETScMatrix.cpp
> dolfin/nls/NewtonSolver.cpp
>
>
> --
> lp:dolfin
> https://code.launchpad.net/~dolfin-core/dolfin/main
>
> Your team DOLFIN Core Team is subscribed to branch lp:dolfin.
> To unsubscribe from this branch go to https://code.launchpad.net/~dolfin-core/dolfin/main/+edit-subscription
>
> === modified file 'demo/pde/advection-diffusion/cpp/main.cpp'
> --- demo/pde/advection-diffusion/cpp/main.cpp 2010-07-17 22:18:58 +0000
> +++ demo/pde/advection-diffusion/cpp/main.cpp 2010-07-22 18:05:26 +0000
> @@ -19,8 +19,8 @@
> int main(int argc, char *argv[])
> {
> //parameters["linear_algebra_backend"] = "uBLAS";
> - //parameters["linear_algebra_backend"] = "PETSc";
> - parameters["linear_algebra_backend"] = "Epetra";
> + parameters["linear_algebra_backend"] = "PETSc";
> + //parameters["linear_algebra_backend"] = "Epetra";
>
> // Read mesh
> Mesh mesh("../mesh.xml.gz");
> @@ -69,7 +69,7 @@
> lu.parameters["reuse_factorization"] = true;
>
> // Parameters for time-stepping
> - double T = 2.0;
> + double T = 20.0;
> const double k = 0.05;
> double t = k;
>
> @@ -94,6 +94,7 @@
> p = t / T;
> t += k;
> }
> + std::cout << "Test " << u.vector().norm("l2") << std::endl;
>
> // Plot solution
> plot(u);
>
> === modified file 'dolfin/la/PETScBaseMatrix.h'
> --- dolfin/la/PETScBaseMatrix.h 2010-07-21 17:31:22 +0000
> +++ dolfin/la/PETScBaseMatrix.h 2010-07-22 18:05:26 +0000
> @@ -50,11 +50,14 @@
> /// Return number of rows (dim = 0) or columns (dim = 1) along dimension dim
> uint size(uint dim) const
> {
> - assert(A);
> - int M = 0;
> - int N = 0;
> - MatGetSize(*A, &M, &N);
> - return (dim == 0 ? M : N);
> + if (A)
> + {
> + int M(0), N(0);
> + MatGetSize(*A, &M, &N);
> + return (dim == 0 ? M : N);
> + }
> + else
> + return 0;
> }
>
> /// Return PETSc Mat pointer
>
> === modified file 'dolfin/la/PETScKrylovSolver.cpp'
> --- dolfin/la/PETScKrylovSolver.cpp 2010-07-21 17:31:22 +0000
> +++ dolfin/la/PETScKrylovSolver.cpp 2010-07-22 18:05:26 +0000
> @@ -112,10 +112,6 @@
> //-----------------------------------------------------------------------------
> void PETScKrylovSolver::set_operator(const PETScBaseMatrix& A)
> {
> - // Check dimensions
> - if (A.size(0) != A.size(1))
> - error("Cannot solve non-square PETSc matrix.");
> -
> this->A = reference_to_no_delete_pointer(A);
> assert(this->A);
>
>
> === modified file 'dolfin/la/PETScKrylovSolver.h'
> --- dolfin/la/PETScKrylovSolver.h 2010-07-21 17:31:22 +0000
> +++ dolfin/la/PETScKrylovSolver.h 2010-07-22 18:05:26 +0000
> @@ -77,12 +77,12 @@
> /// Return informal string representation (pretty-print)
> std::string str(bool verbose) const;
>
> + /// Return PETSc KSP pointer
> + boost::shared_ptr<KSP> ksp() const;
> +
> /// Default parameter values
> static Parameters default_parameters();
>
> - /// Return PETSc KSP pointer
> - boost::shared_ptr<KSP> ksp() const;
> -
> private:
>
> /// Initialize KSP solver
>
> === modified file 'dolfin/la/PETScLUSolver.cpp'
> --- dolfin/la/PETScLUSolver.cpp 2010-07-16 11:45:36 +0000
> +++ dolfin/la/PETScLUSolver.cpp 2010-07-22 18:05:26 +0000
> @@ -75,7 +75,7 @@
> init_solver(lu_package);
> }
> //-----------------------------------------------------------------------------
> -PETScLUSolver::PETScLUSolver(boost::shared_ptr<const GenericMatrix> A,
> +PETScLUSolver::PETScLUSolver(boost::shared_ptr<const PETScMatrix> A,
> std::string lu_package) : A(A)
> {
> // Check dimensions
> @@ -96,31 +96,18 @@
> //-----------------------------------------------------------------------------
> void PETScLUSolver::set_operator(const GenericMatrix& A)
> {
> - // Check dimensions
> - if (A.size(0) != A.size(1))
> - error("Cannot LU factorize non-square PETSc matrix.");
> -
> + set_operator(A.down_cast<PETScMatrix>());
> +}
> +//-----------------------------------------------------------------------------
> +void PETScLUSolver::set_operator(const PETScMatrix& A)
> +{
> this->A = reference_to_no_delete_pointer(A);
> assert(this->A);
> -
> - // Downcast matrix
> - const PETScMatrix& _A = A.down_cast<PETScMatrix>();
> -
> - // Get some parameters
> - const bool reuse_fact = parameters["reuse_factorization"];
> - const bool same_pattern = parameters["same_nonzero_pattern"];
> -
> - // Set operators with appropriate option
> - if (reuse_fact)
> - KSPSetOperators(*ksp, *_A.mat(), *_A.mat(), SAME_PRECONDITIONER);
> - else if (same_pattern)
> - KSPSetOperators(*ksp, *_A.mat(), *_A.mat(), SAME_NONZERO_PATTERN);
> - else
> - KSPSetOperators(*ksp, *_A.mat(), *_A.mat(), DIFFERENT_NONZERO_PATTERN);
> }
> //-----------------------------------------------------------------------------
> dolfin::uint PETScLUSolver::solve(GenericVector& x, const GenericVector& b)
> {
> + assert(_ksp);
> assert(A);
>
> // Downcast matrix and vectors
> @@ -135,14 +122,13 @@
> x.resize(A->size(1));
>
> // Set PETSc operators (depends on factorization re-use options);
> - const PETScMatrix& _A = this->A->down_cast<PETScMatrix>();
> - set_petsc_operators(_A);
> + set_petsc_operators();
>
> // Write a pre-solve message
> pre_report(A->down_cast<PETScMatrix>());
>
> // Solve linear system
> - KSPSolve(*ksp, *_b.vec(), *_x.vec());
> + KSPSolve(*_ksp, *_b.vec(), *_x.vec());
>
> return 1;
> }
> @@ -167,7 +153,7 @@
> pre_report(A);
>
> // Solve linear system
> - KSPSolve(*ksp, *b.vec(), *x.vec());
> + KSPSolve(*_ksp, *b.vec(), *x.vec());
>
> return 1;
> }
> @@ -179,7 +165,7 @@
> if (verbose)
> {
> warning("Verbose output for PETScLUSolver not implemented, calling PETSc KSPView directly.");
> - KSPView(*ksp, PETSC_VIEWER_STDOUT_WORLD);
> + KSPView(*_ksp, PETSC_VIEWER_STDOUT_WORLD);
> }
> else
> s << "<PETScLUSolver>";
> @@ -187,6 +173,11 @@
> return s.str();
> }
> //-----------------------------------------------------------------------------
> +boost::shared_ptr<KSP> PETScLUSolver::ksp() const
> +{
> + return _ksp;
> +}
> +//-----------------------------------------------------------------------------
> const MatSolverPackage PETScLUSolver::select_solver(std::string& lu_package) const
> {
> // Check package string
> @@ -221,25 +212,25 @@
> const MatSolverPackage solver_package = select_solver(lu_package);
>
> // Destroy old solver environment if necessary
> - if (ksp)
> + if (_ksp)
> {
> - if (!ksp.unique())
> + if (!_ksp.unique())
> error("Cannot create new KSP Krylov solver. More than one object points to the underlying PETSc object.");
> }
> - ksp.reset(new KSP, PETScKSPDeleter());
> + _ksp.reset(new KSP, PETScKSPDeleter());
>
> // Create solver
> if (MPI::num_processes() > 1)
> - KSPCreate(PETSC_COMM_WORLD, ksp.get());
> + KSPCreate(PETSC_COMM_WORLD, _ksp.get());
> else
> - KSPCreate(PETSC_COMM_SELF, ksp.get());
> + KSPCreate(PETSC_COMM_SELF, _ksp.get());
>
> // Make solver preconditioner only
> - KSPSetType(*ksp, KSPPREONLY);
> + KSPSetType(*_ksp, KSPPREONLY);
>
> // Set preconditioner to LU factorization
> PC pc;
> - KSPGetPC(*ksp, &pc);
> + KSPGetPC(*_ksp, &pc);
> PCSetType(pc, PCLU);
>
> // Set solver package
> @@ -254,49 +245,28 @@
> #endif
> }
> //-----------------------------------------------------------------------------
> -void PETScLUSolver::set_petsc_operators(const PETScMatrix& A)
> +void PETScLUSolver::set_petsc_operators()
> {
> + assert(A->mat());
> +
> // Get some parameters
> const bool reuse_fact = parameters["reuse_factorization"];
> const bool same_pattern = parameters["same_nonzero_pattern"];
>
> - // Check if operators have been set
> - PetscTruth matset, pmatset;
> - KSPGetOperatorsSet(*ksp, &matset, &pmatset);
> - if (matset == PETSC_TRUE && pmatset == PETSC_TRUE)
> - {
> - // Get preconditioner re-use flag
> - MatStructure mat_struct;
> - Mat Amat, Pmat;
> - KSPGetOperators(*ksp, &Amat, &Pmat, &mat_struct);
> -
> - // Set operators if necessary
> - if (reuse_fact && mat_struct != SAME_PRECONDITIONER)
> - KSPSetOperators(*ksp, *A.mat(), *A.mat(), SAME_PRECONDITIONER);
> - else if (same_pattern && mat_struct != SAME_NONZERO_PATTERN)
> - KSPSetOperators(*ksp, *A.mat(), *A.mat(), SAME_NONZERO_PATTERN);
> - else if (!reuse_fact && !same_pattern && mat_struct != DIFFERENT_NONZERO_PATTERN)
> - KSPSetOperators(*ksp, *A.mat(), *A.mat(), DIFFERENT_NONZERO_PATTERN);
> - }
> - else if (matset != PETSC_TRUE && pmatset != PETSC_TRUE)
> - {
> - // Set operators with appropriate option
> - if (reuse_fact)
> - KSPSetOperators(*ksp, *A.mat(), *A.mat(), SAME_PRECONDITIONER);
> - else if (same_pattern)
> - KSPSetOperators(*ksp, *A.mat(), *A.mat(), SAME_NONZERO_PATTERN);
> - else
> - KSPSetOperators(*ksp, *A.mat(), *A.mat(), DIFFERENT_NONZERO_PATTERN);
> - }
> + // Set operators with appropriate preconditioner option
> + if (reuse_fact)
> + KSPSetOperators(*_ksp, *A->mat(), *A->mat(), SAME_PRECONDITIONER);
> + else if (same_pattern)
> + KSPSetOperators(*_ksp, *A->mat(), *A->mat(), SAME_NONZERO_PATTERN);
> else
> - error("Operators incorrectly set for PETSc.");
> + KSPSetOperators(*_ksp, *A->mat(), *A->mat(), DIFFERENT_NONZERO_PATTERN);
> }
> //-----------------------------------------------------------------------------
> void PETScLUSolver::pre_report(const PETScMatrix& A) const
> {
> const MatSolverPackage solver_type;
> PC pc;
> - KSPGetPC(*ksp, &pc);
> + KSPGetPC(*_ksp, &pc);
> PCFactorGetMatSolverPackage(pc, &solver_type);
>
> // Get parameter
>
> === modified file 'dolfin/la/PETScLUSolver.h'
> --- dolfin/la/PETScLUSolver.h 2010-07-16 11:45:36 +0000
> +++ dolfin/la/PETScLUSolver.h 2010-07-22 18:05:26 +0000
> @@ -42,7 +42,7 @@
> PETScLUSolver(const GenericMatrix& A, std::string lu_package="default");
>
> /// Constructor
> - PETScLUSolver(boost::shared_ptr<const GenericMatrix> A,
> + PETScLUSolver(boost::shared_ptr<const PETScMatrix> A,
> std::string lu_package="default");
>
> /// Destructor
> @@ -51,6 +51,9 @@
> /// Set operator (matrix)
> void set_operator(const GenericMatrix& A);
>
> + /// Set operator (matrix)
> + void set_operator(const PETScMatrix& A);
> +
> /// Solve linear system Ax = b
> uint solve(GenericVector& x, const GenericVector& b);
>
> @@ -63,6 +66,9 @@
> /// Return informal string representation (pretty-print)
> std::string str(bool verbose) const;
>
> + /// Return PETSc KSP pointer
> + boost::shared_ptr<KSP> ksp() const;
> +
> /// Default parameter values
> static Parameters default_parameters();
>
> @@ -78,16 +84,16 @@
> void init_solver(std::string& lu_package);
>
> // Set PETSc operators
> - void set_petsc_operators(const PETScMatrix& A);
> + void set_petsc_operators();
>
> - // Pritn pre-solve report
> + // Print pre-solve report
> void pre_report(const PETScMatrix& A) const;
>
> /// PETSc solver pointer
> - boost::shared_ptr<KSP> ksp;
> + boost::shared_ptr<KSP> _ksp;
>
> // Operator (the matrix)
> - boost::shared_ptr<const GenericMatrix> A;
> + boost::shared_ptr<const PETScMatrix> A;
>
> };
>
>
> === modified file 'dolfin/la/PETScMatrix.cpp'
> --- dolfin/la/PETScMatrix.cpp 2010-07-21 17:31:22 +0000
> +++ dolfin/la/PETScMatrix.cpp 2010-07-22 18:05:26 +0000
> @@ -70,6 +70,8 @@
> if (A && !A.unique())
> error("Cannot resize PETScMatrix. More than one object points to the underlying PETSc object.");
> A.reset(new Mat, PETScMatrixDeleter());
> + std::cout << "Reset matrix (A)" << std::endl;
> +
>
> // FIXME: maybe 50 should be a parameter?
> // FIXME: it should definitely be a parameter
> @@ -112,6 +114,7 @@
> if (A && !A.unique())
> error("Cannot initialise PETScMatrix. More than one object points to the underlying PETSc object.");
> A.reset(new Mat, PETScMatrixDeleter());
> + std::cout << "Reset matrix (B)" << std::endl;
>
> // Initialize matrix
> if (row_range.first == 0 && row_range.second == M)
> @@ -389,13 +392,13 @@
> // Check for self-assignment
> if (this != &A)
> {
> - if (this->A && !this->A.unique())
> - error("Cannot assign PETScMatrix with different non-zero pattern because more than one object points to the underlying PETSc object.");
> - this->A.reset(new Mat, PETScMatrixDeleter());
> + if (this->A && !this->A.unique())
> + error("Cannot assign PETScMatrix because more than one object points to the underlying PETSc object.");
> + this->A.reset(new Mat, PETScMatrixDeleter());
> + std::cout << "Reset matrix (C)" << std::endl;
>
> - // Duplicate with the same pattern as A.A
> - MatDuplicate(*A.mat(), MAT_COPY_VALUES, this->A.get());
> - //}
> + // Duplicate with the same pattern as A.A
> + MatDuplicate(*A.mat(), MAT_COPY_VALUES, this->A.get());
> }
> return *this;
> }
>
> === modified file 'dolfin/nls/NewtonSolver.cpp'
> --- dolfin/nls/NewtonSolver.cpp 2010-07-22 12:05:22 +0000
> +++ dolfin/nls/NewtonSolver.cpp 2010-07-22 18:05:26 +0000
> @@ -26,7 +26,7 @@
> {
> Parameters p("newton_solver");
>
> - p.add("maximum_iterations", 50);
> + p.add("maximum_iterations", 10);
> p.add("relative_tolerance", 1e-9);
> p.add("absolute_tolerance", 1e-10);
> p.add("convergence_criterion", "residual");
> @@ -35,6 +35,8 @@
> p.add("report", true);
> p.add("error_on_nonconvergence", true);
>
> + //p.add("reuse_preconditioner", false);
> +
> return p;
> }
> //-----------------------------------------------------------------------------
> @@ -81,14 +83,12 @@
> newton_iteration = 0;
> bool newton_converged = false;
>
> - // FIXME: Should the operator be set in the constructor?
> - // Set linear solver operator
> - solver->set_operator(*A);
> -
> // Compute F(u)
> nonlinear_problem.form(*A, *b, x);
> nonlinear_problem.F(*b, x);
>
> + solver->set_operator(*A);
> +
> // Start iterations
> while (!newton_converged && newton_iteration < maxiter)
> {
>
>
>
Follow ups
