[petsc-users] strange and slow behavior when using KSPSolve
Alejandro Marcos Aragón
alejandro.aragon at gmail.com
Thu Mar 31 10:19:26 CDT 2011
Hi all,
I am currently using PETSc to solve a system of linear equations in parallel, so I created a wrapper class that interfaces my C++ code with the PETSc library.
Since the problem I'm trying to solve is nonlinear, I assemble the matrix once and I reuse the structure (the sparsity pattern is the same). The problem is that after a few times of solving the system very fast, it suddenly takes a considerable amount of time and I really don't understand why this is happening.
This is an extract of the code showing the timing issue that I have just before solving the system with 2 processes:
[0] Inside PETSc solver:7.15256e-06 s
[0] Got ownership range... 3.21865e-05 s
Inside PETSc solver:1.5974e-05 s
[1] Got ownership range... 4.00543e-05 s
[1] Zeroed-out sparse matrix entries... [0] Zeroed-out sparse matrix entries... 0.000444174 s
0.000442982 s
[1] Filled sparse matrix... 0.0772262 s
[0] Filled sparse matrix... 0.081254 s
[1] Sparse matrix assembled... 0.0897892 s
[0] Sparse matrix assembled... 0.0954142 s
[0] Vectors created... 0.0955811 s[1] Vectors created... 0.095592 s
[0] Right hand side set... 0.0965121 s
[1] Right hand side set... 0.0965312 s
[0] Right-hand-side vector assembled... 0.096647 s
[1] Right-hand-side vector assembled... 0.096657 s
[0] System solved... 32.6188 s
[1] System solved... 32.6188 s
[0] Solution vector scattered... [1] Solution vector scattered... 32.6192 s
32.6192 s
[1] Solution vector copied... 32.6215[0] Solution vector copied... 32.6215 s
s
[1] Temporary data structures destroyed... 32.6216 s[0] Temporary data structures destroyed... 32.6216 s
So you can see that it takes a considerable amount of time to solve the system at this point, when usually this is what I get:
[0] Right-hand-side vector assembled... 0.089427 s
[1] Right-hand-side vector assembled... 0.089438 s
[0] System solved... 0.383509 s
[1] System solved... 0.38352 s
Does anyone know why this is happening? Thanks in advance for any insight on this.
Alejandro M. Aragón
P.S.: The code I use is given below:
struct PETSc_solver {
typedef doVERBOSEuble value_type;
typedef sparse_vector<value_type> sparse_vector_type;
typedef sparse_matrix<value_type> sparse_matrix_type;
int argc_;
char **args_;
Mat A_; //!< linear system matrix
Vec x_; //!< Solution vector
KSP ksp_; //!< linear solver context
bool allocated_;
PETSc_solver() : argc_(), args_(NULL), allocated_(false) {
PetscInitialize(&argc_,&args_,NULL,NULL);
PetscErrorCode ierr = MatCreate(PETSC_COMM_WORLD,&A_);CHKERR(ierr);
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp_);CHKERR(ierr);
}
PETSc_solver(int argc, char** args) : argc_(argc), args_(args) {
PetscInitialize(&argc_,&args_,NULL,NULL);
PetscErrorCode ierr = MatCreate(PETSC_COMM_WORLD,&A_);CHKERR(ierr);
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp_);CHKERR(ierr);
}
//! Overload operator() to solve system of linear equations
sparse_vector_type operator()(const sparse_matrix_type& AA, const sparse_vector_type& bb);
~PETSc_solver() {
PetscErrorCode ierr = MatDestroy(A_);CHKERR(ierr);
ierr = KSPDestroy(ksp_);CHKERR(ierr);
ierr = VecDestroy(x_);CHKERR(ierr);
}
};
PETSc_solver::sparse_vector_type PETSc_solver::operator()(const PETSc_solver::sparse_matrix_type& AA, const PETSc_solver::sparse_vector_type& bb) {
assert(AA.rows() == bb.size());
Vec x,b; /* approx solution, RHS */
PetscInt Istart,Iend;
PetscErrorCode ierr;
#ifdef VERBOSE
cpputils::ctimer timer;
MPIPRINT("Inside PETSc solver:" );
cout<<timer<<endl;
#endif
ierr = MatSetFromOptions(A_);CHKERR(ierr);
if (!allocated_) {
#ifdef VERBOSE
MPIPRINT(" Allocating memory... ");
cout<<timer<<endl;
#endif
ierr = MatSetSizes(A_,PETSC_DECIDE,PETSC_DECIDE, AA.rows(), AA.columns());CHKERR(ierr);
// get number of non-zeros in both diagonal and non-diagonal portions of the matrix
std::pair<size_t,size_t> nz = AA.non_zero_off_diagonal();
ierr = MatMPIAIJSetPreallocation(A_, nz.first, PETSC_NULL, nz.second, PETSC_NULL); CHKERR(ierr);
ierr = MatSeqAIJSetPreallocation(A_, nz.first, PETSC_NULL); CHKERR(ierr);
#ifdef VERBOSE
cout<<"Ok! "<<timer<<endl;
#endif
}
/*
Currently, all PETSc parallel matrix formats are partitioned by
contiguous chunks of rows across the processors. Determine which
rows of the matrix are locally owned.
*/
ierr = MatGetOwnershipRange(A_,&Istart,&Iend);CHKERR(ierr);
#ifdef VERBOSE
MPIPRINT(" Got ownership range... ");
cout<<timer<<endl;
#endif
// zero-out matrix
MatZeroEntries(A_);
#ifdef VERBOSE
MPIPRINT(" Zeroed-out sparse matrix entries... ");
cout<<timer<<endl;
#endif
for (sparse_matrix_type::const_hash_iterator it = AA.map_.begin(); it != AA.map_.end(); ++it) {
// get subscripts
std::pair<size_t,size_t> subs = AA.unhash(it->first);
int row = subs.first;
int col = subs.second;
ierr = MatSetValues(A_, 1, &row, 1, &col, &it->second, ADD_VALUES);CHKERR(ierr);
}
#ifdef VERBOSE
MPIPRINT(" Filled sparse matrix... ");
cout<<timer<<endl;
#endif
/*
Assemble matrix, using the 2-step process:
MatAssemblyBegin(), MatAssemblyEnd()
Computations can be done while messages are in transition
by placing code between these two statements.
*/
ierr = MatAssemblyBegin(A_,MAT_FINAL_ASSEMBLY);CHKERR(ierr);
ierr = MatAssemblyEnd(A_,MAT_FINAL_ASSEMBLY);CHKERR(ierr);
#ifdef VERBOSE
MPIPRINT(" Sparse matrix assembled... ");
cout<<timer<<endl;
#endif
ierr = MatSetOption(A_, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE);CHKERR(ierr);
/* Create parallel vectors. */
ierr = VecCreate(PETSC_COMM_WORLD,&b);CHKERR(ierr);
ierr = VecSetSizes(b,PETSC_DECIDE, bb.size());CHKERR(ierr);
ierr = VecSetFromOptions(b);CHKERR(ierr);
if (!allocated_) {
ierr = VecDuplicate(b,&x_);CHKERR(ierr);
// set flag to true
allocated_ = true;
}
VecZeroEntries(x_);
#ifdef VERBOSE
MPIPRINT(" Vectors created... ");
cout<<timer<<endl;
#endif
/* Set hight-hand-side vector */
for (sparse_vector_type::const_hash_iterator it = bb.map_.begin(); it != bb.map_.end(); ++it) {
int row = it->first;
ierr = VecSetValues(b, 1, &row, &it->second, ADD_VALUES);
}
#ifdef VERBOSE
MPIPRINT(" Right hand side set... ");
cout<<timer<<endl;
#endif
/*
Assemble vector, using the 2-step process:
VecAssemblyBegin(), VecAssemblyEnd()
Computations can be done while messages are in transition
by placing code between these two statements.
*/
ierr = VecAssemblyBegin(b);CHKERR(ierr);
ierr = VecAssemblyEnd(b);CHKERR(ierr);
#ifdef VERBOSE
MPIPRINT(" Right-hand-side vector assembled... ");
cout<<timer<<endl;
#endif
/*
Set operators. Here the matrix that defines the linear system
also serves as the preconditioning matrix.
*/
ierr = KSPSetOperators(ksp_,A_,A_,SAME_NONZERO_PATTERN);CHKERR(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Solve the linear system
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = KSPSolve(ksp_,b,x_);CHKERR(ierr);
#ifdef VERBOSE
MPIPRINT(" System solved... ");
cout<<timer<<endl;
#endif
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Check solution and clean up
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
VecScatter ctx;
Vec x_all;
ierr = VecScatterCreateToAll(x_, &ctx, &x_all);CHKERR(ierr);
ierr = VecScatterBegin(ctx,x_,x_all,INSERT_VALUES,SCATTER_FORWARD);CHKERR(ierr);
ierr = VecScatterEnd(ctx,x_,x_all,INSERT_VALUES,SCATTER_FORWARD);CHKERR(ierr);
#ifdef VERBOSE
MPIPRINT(" Solution vector scattered... ");
cout<<timer<<endl;
#endif
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Get values from solution and store them in the object that will be
returned
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
sparse_vector_type xx(bb.size());
/* Set solution vector */
for (sparse_vector_type::const_hash_iterator it = bb.map_.begin(); it != bb.map_.end(); ++it) {
int row = it->first;
ierr = VecGetValues(x_all, 1, &row, &xx[row]);
}
#ifdef VERBOSE
MPIPRINT(" Solution vector copied... ");
cout<<timer<<endl;
#endif
ierr = VecDestroy(b);CHKERR(ierr);
ierr = VecDestroy(x_all);CHKERR(ierr);
#ifdef VERBOSE
MPIPRINT(" Temporary data structures destroyed... ");
cout<<timer<<endl;
#endif
return xx;
}
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20110331/aa4bf1d9/attachment-0001.htm>
More information about the petsc-users
mailing list