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, <nore...@launchpad.net> wrote: > ------------------------------------------------------------ > revno: 4868 > committer: Garth N. Wells <gn...@cam.ac.uk> > 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) > { > > > _______________________________________________ Mailing list: https://launchpad.net/~dolfin Post to : dolfin@lists.launchpad.net Unsubscribe : https://launchpad.net/~dolfin More help : https://help.launchpad.net/ListHelp
