dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

Here's some sample code that illustrates the problem.

We're using a bunch of wrappers in DOLFIN to make life easier for us
and the GMRES solver is used repeatedly in the ODE solver where the
size of the system may change between different time steps (time
slabs):

    Matrix A;    // Created once, resized when necessary
    Vector x;    // Created once, resized when necessary
    Vector b;    // Created once, resized when necessary
    GMRES gmres; // Created once

    gmres.solve(A, x, b); // Called repeatedly

The attached code summarizes the PETSc calls that this gives rise to.

/Anders

On Tue, Oct 18, 2005 at 11:13:34AM -0500, Anders Logg wrote:

> Do I have to call KSPDestroy(), KSPCreate() when I change the size of
> a linear system?
> 
> If I set up a Krylov solver using KSPCreate() and then solve a linear
> system Ax = b with
> 
>     KSPSetOperators(ksp, A, A, DIFFERENT_NONZERO_PATTERN);
>     KSPSolve(ksp, b, x);
> 
> and then change the size of A, x, b and try to do
> 
>     KSPSetOperators(ksp, A, A, DIFFERENT_NONZERO_PATTERN);
>     KSPSolve(ksp, b, x);
> 
> again, I get an error message: "Nonconforming object sizes!"
> 
> From what I can see, the vector VEC_TEMP still has the old size when
> KSPInitialResidual gets called in gmres.c.
> 
> It works if I destroy and again create KSP pointer when the size
> changes.
> 
> /Anders

-- 
Anders Logg
Research Assistant Professor
Toyota Technological Institute at Chicago
http://www.tti-c.org/logg/

int usermult(Mat A, Vec x, Vec y)
{
  VecCopy(x, y);
}

int main(int argc, char* argv[])
{
  PetscInitialize(&argc, &argv, PETSC_NULL, PETSC_NULL);
  
  // Create vector x
  Vec x;
  VecCreate(PETSC_COMM_SELF, &x);
  VecSetSizes(x, PETSC_DECIDE, 10);
  VecSetFromOptions(x);
  VecSet(x, 0.0);
  
  // Create vector b
  Vec b;
  VecCreate(PETSC_COMM_SELF, &b);
  VecSetSizes(b, PETSC_DECIDE, 10);
  VecSetFromOptions(b);
  VecSet(b, 1.0);
  
  VecView(x, PETSC_VIEWER_STDOUT_SELF);
  VecView(b, PETSC_VIEWER_STDOUT_SELF);
  
  // Create matrix A
  int m, n, M, N;
  VecGetLocalSize(b, &m);
  VecGetLocalSize(x, &n);
  VecGetSize(b, &M);
  VecGetSize(x, &N);  
  Mat A;
  MatCreateShell(PETSC_COMM_SELF, m, n, M, N, NULL, &A);
  MatShellSetOperation(A, MATOP_MULT, (void (*)()) usermult);
  
  // Set up GMRES solver
  KSP ksp;
  KSPCreate(PETSC_COMM_SELF, &ksp);
  KSPSetFromOptions(ksp);  
  KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
  
  // Chooose preconditioner (none)
  PC pc;
  KSPGetPC(ksp, &pc);
  PCSetType(pc, PCNONE);
  
  // Solve
  KSPSetOperators(ksp, A, A, DIFFERENT_NONZERO_PATTERN);
  KSPSolve(ksp, b, x);
  
  VecView(x, PETSC_VIEWER_STDOUT_SELF);
  VecView(b, PETSC_VIEWER_STDOUT_SELF);
  
  // Change size of vector x
  VecDestroy(x);
  VecCreate(PETSC_COMM_SELF, &x);
  VecSetSizes(x, PETSC_DECIDE, 11);
  VecSetFromOptions(x);
  VecSet(x, 0.0);
  
  // Change size of vector b
  VecDestroy(b);
  VecCreate(PETSC_COMM_SELF, &b);
  VecSetSizes(b, PETSC_DECIDE, 11);
  VecSetFromOptions(b);
  VecSet(b, 0.0);
  
  // Change size of matrix A
  VecGetLocalSize(b, &m);
  VecGetLocalSize(x, &n);
  VecGetSize(b, &M);
  VecGetSize(x, &N);  
  MatDestroy(A);
  MatCreateShell(PETSC_COMM_SELF, m, n, M, N, NULL, &A);
  MatShellSetOperation(A, MATOP_MULT, (void (*)()) usermult);
  
  /*
    KSPDestroy(ksp);
    KSPCreate(PETSC_COMM_SELF, &ksp);
    KSPSetFromOptions(ksp);  
    KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
    KSPGetPC(ksp, &pc);
    PCSetType(pc, PCNONE);
  */
  
  // Solve again
  KSPSetOperators(ksp, A, A, DIFFERENT_NONZERO_PATTERN);
  KSPSolve(ksp, b, x);
  
  VecView(x, PETSC_VIEWER_STDOUT_SELF);
  VecView(b, PETSC_VIEWER_STDOUT_SELF);
  
  // Destroy matrix and vectors
  MatDestroy(A);
  VecDestroy(x);
  VecDestroy(b);

  return 0;
}