On Wed, 15 Sep 2010 14:53:18 +0200 (CEST), Tim Kroeger
<[email protected]> wrote:
> > This chooses the sizes of each part, but you will still have to move
> > indices around.
>
> This is done by VecScatter, isn't it? I think that the overhead of
> doing this once before and after the solve will be small against the
> gain in performance by having balanced vectors. Well, of course this
> depends on a lot of things, but in the mean...
VecScatter moves entries around, and isn't the issue here.
Redistributing indices to satisfy some a-priori size bound is somewhat
different (though not hard).
> Anyway, I have now attached the current state of my extension.
> Questions are:
>
> * (mainly to Jed:) Are the things in petsc_linear_solver.{h,C} now
> correct, sensible, and reasonably efficient? If yes, I will do the
> same things to the remaining overloaded solve() methods.
> Index: include/solvers/petsc_linear_solver.h
> ===================================================================
> --- include/solvers/petsc_linear_solver.h (Revision 3957)
> +++ include/solvers/petsc_linear_solver.h (Arbeitskopie)
> @@ -122,7 +122,15 @@
> * Initialize data structures if not done so already plus much more
> */
> void init (PetscMatrix<T>* matrix);
> +
> /**
> + * After calling this method, all successive solves will be limited
> + * to the given set of dofs. This mode can be disabled by calling
> + * this method with \p dofs being a \p NULL pointer.
> + */
> + virtual void solve_only_on (const std::vector<unsigned int>* const dofs);
> +
> + /**
> * Call the Petsc solver. It calls the method below, using the
> * same matrix for the system and preconditioner matrices.
> */
> @@ -253,13 +261,35 @@
> * Krylov subspace context
> */
> KSP _ksp;
> +
> + /**
> + * PETSc index set containing the dofs on which to solve (\p NULL
> + * means solve on all dofs).
> + */
> + IS _solve_only_on_is;
> +
> + /**
> + * PETSc index set containing the <em>local</em> dofs on which to
> + * solve. This index set is determined when a system is solved on a
> + * subset the first time. Before that point (or if no subset is
> + * selected), it will be a \p NULL pointer.
> + */
> + IS _solve_only_on_is_local;
I don't think you want or need this, see below.
> +
> + /**
> + * Internal method that returns the size of \p _solve_only_on_is.
> + */
> + size_t _solve_only_on_is_size(void)const;
Note that you only call this once.
> +
> };
>
>
> /*----------------------- functions ----------------------------------*/
> template <typename T>
> inline
> -PetscLinearSolver<T>::PetscLinearSolver ()
> + PetscLinearSolver<T>::PetscLinearSolver ():
> + _solve_only_on_is(NULL),
> + _solve_only_on_is_local(NULL)
> {
> if (libMesh::n_processors() == 1)
> this->_preconditioner_type = ILU_PRECOND;
> @@ -277,6 +307,23 @@
> }
>
>
> +
> +template <typename T>
> +inline size_t
> +PetscLinearSolver<T>::
> +_solve_only_on_is_size(void)const
> +{
> + libmesh_assert(_solve_only_on_is!=NULL);
> +
> + PetscInt s;
> + int ierr = ISGetSize(_solve_only_on_is,&s);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> +
> + return static_cast<size_t>(s);
> +}
> +
> +
> +
> } // namespace libMesh
>
>
> Index: src/solvers/petsc_linear_solver.C
> ===================================================================
> --- src/solvers/petsc_linear_solver.C (Revision 3957)
> +++ src/solvers/petsc_linear_solver.C (Arbeitskopie)
> @@ -98,6 +98,8 @@
> {
> if (this->initialized())
> {
> + this->solve_only_on(NULL);
> +
> this->_is_initialized = false;
>
> int ierr=0;
> @@ -294,7 +296,7 @@
> CHKERRABORT(libMesh::COMM_WORLD,ierr);
>
> // Set operators. The input matrix works as the preconditioning matrix
> - ierr = KSPSetOperators(_ksp, matrix->mat(),
> matrix->mat(),SAME_NONZERO_PATTERN);
> + ierr = KSPSetOperators(_ksp, matrix->mat(),
> matrix->mat(),DIFFERENT_NONZERO_PATTERN);
> CHKERRABORT(libMesh::COMM_WORLD,ierr);
>
> // Set user-specified solver and preconditioner types
> @@ -360,9 +362,46 @@
>
>
>
> +template <typename T>
> +void
> +PetscLinearSolver<T>::solve_only_on (const std::vector<unsigned int>* const
> dofs)
> +{
> + libmesh_assert(this->initialized());
>
> + int ierr=0;
>
> + if(_solve_only_on_is!=NULL)
> + {
> + ierr = ISDestroy(_solve_only_on_is);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + _solve_only_on_is = NULL;
> + }
>
> + if(_solve_only_on_is_local!=NULL)
> + {
> + ierr = ISDestroy(_solve_only_on_is_local);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + _solve_only_on_is_local = NULL;
> + }
> +
> + if(dofs!=NULL)
> + {
> + PetscInt* petsc_dofs;
> + ierr = PetscMalloc(dofs->size()*sizeof(PetscInt),&petsc_dofs);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> +
> + for(size_t i=0; i<dofs->size(); i++)
> + {
> + petsc_dofs[i] = (*dofs)[i];
> + }
> +
> + ierr =
> ISCreateGeneralNC(libMesh::COMM_WORLD,dofs->size(),petsc_dofs,&_solve_only_on_is);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + }
> +}
> +
> +
> +
> template <typename T>
> std::pair<unsigned int, Real>
> PetscLinearSolver<T>::solve (SparseMatrix<T>& matrix_in,
> @@ -406,10 +445,15 @@
>
> // 2.1.x & earlier style
> #if PETSC_VERSION_LESS_THAN(2,2,0)
> +
> + if(_solve_only_on_is!=NULL)
> + {
> + libmesh_not_implemented();
> + }
>
> // Set operators. The input matrix works as the preconditioning matrix
> ierr = SLESSetOperators(_sles, matrix->mat(), precond->mat(),
> - SAME_NONZERO_PATTERN);
> + DIFFERENT_NONZERO_PATTERN);
> CHKERRABORT(libMesh::COMM_WORLD,ierr);
>
> // Set the tolerances for the iterative solver. Use the user-supplied
> @@ -431,6 +475,11 @@
> // 2.2.0
> #elif PETSC_VERSION_LESS_THAN(2,2,1)
>
> + if(_solve_only_on_is!=NULL)
> + {
> + libmesh_not_implemented();
> + }
> +
> // Set operators. The input matrix works as the preconditioning matrix
> //ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(),
> // SAME_NONZERO_PATTERN);
> @@ -475,19 +524,77 @@
> // 2.2.1 & newer style
> #else
>
> - // Set operators. The input matrix works as the preconditioning matrix
> - if(!this->same_preconditioner)
> - {
> - ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(),
> - SAME_NONZERO_PATTERN);
> - CHKERRABORT(libMesh::COMM_WORLD,ierr);
> - }
> + Mat submat = NULL;
> + Mat subprecond = NULL;
> + Vec subrhs;
> + Vec subsolution;
> + VecScatter scatter;
> +
> + // Set operators. Also restrict rhs and solution vector to
> + // subdomain if neccessary.
> + if(_solve_only_on_is!=NULL)
> + {
> + size_t is_size = this->_solve_only_on_is_size();
> +
> + ierr = VecCreate(libMesh::COMM_WORLD,&subrhs);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + ierr = VecSetSizes(subrhs,PETSC_DECIDE,is_size);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + ierr = VecSetFromOptions(subrhs);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> +
> + ierr = VecCreate(libMesh::COMM_WORLD,&subsolution);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + ierr = VecSetSizes(subsolution,PETSC_DECIDE,is_size);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + ierr = VecSetFromOptions(subsolution);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> +
> + ierr = VecScatterCreate(rhs->vec(),_solve_only_on_is, subrhs,NULL,
> &scatter);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> +
> + ierr =
> VecScatterBegin(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + ierr =
> VecScatterEnd(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> +
> + ierr =
> VecScatterBegin(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + ierr =
> VecScatterEnd(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> +
> + /* Make the local index set if it doesn't exist yet. */
> + if(_solve_only_on_is_local==NULL)
> + {
> + PetscInt firstIndex = 0;
> + PetscInt lastIndex = 0;
> + ierr = VecGetOwnershipRange(subrhs,&firstIndex,&lastIndex);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> +
> + ierr =
> ISCreateStride(libMesh::COMM_WORLD,lastIndex-firstIndex,firstIndex,1,&_solve_only_on_is_local);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + }
> +
> + ierr = MatGetSubMatrix(matrix->mat(),
> + _solve_only_on_is_local,_solve_only_on_is,
> + PETSC_DECIDE,MAT_INITIAL_MATRIX,&submat);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
Why do you want this non-square matrix? I'm pretty sure you want to
extract the symmetric block (_solve_only_on_is,_solve_only_on_is).
> +
> + ierr = MatGetSubMatrix(precond->mat(),
> + _solve_only_on_is_local,_solve_only_on_is,
> + PETSC_DECIDE,MAT_INITIAL_MATRIX,&subprecond);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> +
> + ierr = KSPSetOperators(_ksp, submat, subprecond,
> + this->same_preconditioner ? SAME_PRECONDITIONER :
> DIFFERENT_NONZERO_PATTERN);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + }
> else
> - {
> - ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(),
> - SAME_PRECONDITIONER);
> - CHKERRABORT(libMesh::COMM_WORLD,ierr);
> - }
> + {
> + ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(),
> + this->same_preconditioner ? SAME_PRECONDITIONER :
> DIFFERENT_NONZERO_PATTERN);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + }
>
> // Set the tolerances for the iterative solver. Use the user-supplied
> // tolerance for the relative residual & leave the others at default
> values.
> @@ -496,8 +603,16 @@
> CHKERRABORT(libMesh::COMM_WORLD,ierr);
>
> // Solve the linear system
> - ierr = KSPSolve (_ksp, rhs->vec(), solution->vec());
> - CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + if(_solve_only_on_is!=NULL)
> + {
> + ierr = KSPSolve (_ksp, subrhs, subsolution);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + }
> + else
> + {
> + ierr = KSPSolve (_ksp, rhs->vec(), solution->vec());
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + }
>
> // Get the number of iterations required for convergence
> ierr = KSPGetIterationNumber (_ksp, &its);
> @@ -507,6 +622,26 @@
> ierr = KSPGetResidualNorm (_ksp, &final_resid);
> CHKERRABORT(libMesh::COMM_WORLD,ierr);
>
> + if(_solve_only_on_is!=NULL)
> + {
> + ierr =
> VecScatterBegin(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + ierr =
> VecScatterEnd(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
You probably want to set the rest of solution->vec() to the right side,
or zero it out, otherwise a solve is not a linear operation. The reason
for the first is that it makes the preconditioner look like the
submatrix was extended to operate on the full domain using the identity
(instead of zero which would make it very singular).
> +
> + ierr = VecScatterDestroy(scatter);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> +
> + ierr = VecDestroy(subsolution);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + ierr = VecDestroy(subrhs);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + ierr = MatDestroy(submat);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + ierr = MatDestroy(subprecond);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + }
> +
> #endif
>
> STOP_LOG("solve()", "PetscLinearSolver");
> @@ -574,7 +709,7 @@
>
> // Set operators. The input matrix works as the preconditioning matrix
> ierr = KSPSetOperators(_ksp, mat, mat,
> - SAME_NONZERO_PATTERN);
> + DIFFERENT_NONZERO_PATTERN);
> CHKERRABORT(libMesh::COMM_WORLD,ierr);
>
> // Set the tolerances for the iterative solver. Use the user-supplied
> Index: src/solvers/linear_solver.C
> ===================================================================
> --- src/solvers/linear_solver.C (Revision 3957)
> +++ src/solvers/linear_solver.C (Arbeitskopie)
> @@ -115,7 +115,18 @@
> _preconditioner = preconditioner;
> }
>
> +template <typename T>
> +void
> +LinearSolver<T>::solve_only_on(const std::vector<unsigned int>* const dofs)
> +{
> + if(dofs!=NULL)
> + {
> + libmesh_not_implemented();
> + }
> +}
>
> +
> +
> //------------------------------------------------------------------
> // Explicit instantiations
> template class LinearSolver<Number>;
------------------------------------------------------------------------------
Start uncovering the many advantages of virtual appliances
and start using them to simplify application deployment and
accelerate your shift to cloud computing.
http://p.sf.net/sfu/novell-sfdev2dev
_______________________________________________
Libmesh-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-users
