[petsc-users] Question about TSComputeRHSJacobianConstant

[Smith, Barry F.]

  Sajid,

    This is a huge embarrassing performance bug in PETSc https://bitbucket.org/petsc/petsc/issues/293/refactoring-of-ts-handling-of-reuse-of

    It is using 74 percent of the time to perform MatAXPY() on two large sparse matrices, not knowing they have identical nonzero patterns and one of which has all zeros off of the diagonal. This despite the fact that a few lines higher in the code is special purpose code for exactly the case you have that  only stores one matrix and only ever shifts the diagonal of the matrix. 

   Please edit TSSetUp() and remove the lines 
  if (ts->rhsjacobian.reuse && rhsjac == TSComputeRHSJacobianConstant) {
    Mat Amat,Pmat;
    SNES snes;
    ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr);
    ierr = SNESGetJacobian(snes,&Amat,&Pmat,NULL,NULL);CHKERRQ(ierr);
    /* Matching matrices implies that an IJacobian is NOT set, because if it had been set, the IJacobian's matrix would
     * have displaced the RHS matrix */
    if (Amat && Amat == ts->Arhs) {
      /* we need to copy the values of the matrix because for the constant Jacobian case the user will never set the numerical values in this new location */
      ierr = MatDuplicate(ts->Arhs,MAT_COPY_VALUES,&Amat);CHKERRQ(ierr);
      ierr = SNESSetJacobian(snes,Amat,NULL,NULL,NULL);CHKERRQ(ierr);
      ierr = MatDestroy(&Amat);CHKERRQ(ierr);
    }
    if (Pmat && Pmat == ts->Brhs) {
      ierr = MatDuplicate(ts->Brhs,MAT_COPY_VALUES,&Pmat);CHKERRQ(ierr);
      ierr = SNESSetJacobian(snes,NULL,Pmat,NULL,NULL);CHKERRQ(ierr);
      ierr = MatDestroy(&Pmat);CHKERRQ(ierr);
    }
  }

You will be stunned by the improvement in time. 


> On May 16, 2019, at 3:06 PM, Sajid Ali via petsc-users <petsc-users at mcs.anl.gov> wrote:
> 
> Hi PETSc developers, 
> 
> I have a question about TSComputeRHSJacobianConstant. If I create a TS (of type linear) for a problem where the jacobian does not change with time (set with the aforementioned option) and run it for different number of time steps, why does the time it takes to evaluate the jacobian change (as indicated by TSJacobianEval) ? 
> 
> To clarify, I run with the example with different TSSetTimeStep, but the same jacobian matrix. I see that the time spent in KSPSolve increases with increasing number of steps (which is as expected as this is a KSPOnly SNES solver). But surprisingly, the time spent in TSJacobianEval also increases with decreasing time-step (or increasing number of steps). 
> 
> For reference, I attach the log files for two cases which were run with different time steps and the source code. 
> 
> Thank You,
> Sajid Ali
> Applied Physics
> Northwestern University
> <ex_dmda.c><out_50><out_100>