1: static char help[] = "Tests that SNESComputeJacobian() calls the user Jacobian function when a left NPC is active.\n\n";

  3: #include <petscsnes.h>

  5: typedef struct {
  6:   PetscInt jac_calls;
  7: } AppCtx;

  9: PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, PetscCtx ctx)
 10: {
 11:   PetscFunctionBeginUser;
 12:   PetscCall(VecZeroEntries(f));
 13:   PetscFunctionReturn(PETSC_SUCCESS);
 14: }

 16: PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, PetscCtx ctx)
 17: {
 18:   AppCtx *appctx = (AppCtx *)ctx;

 20:   PetscFunctionBeginUser;
 21:   appctx->jac_calls++;
 22:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
 23:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
 24:   if (jac != B) {
 25:     PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
 26:     PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
 27:   }
 28:   PetscFunctionReturn(PETSC_SUCCESS);
 29: }

 31: int main(int argc, char **argv)
 32: {
 33:   SNES   snes, npc;
 34:   Vec    x, r;
 35:   Mat    J;
 36:   AppCtx appctx;

 38:   PetscFunctionBeginUser;
 39:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));

 41:   appctx.jac_calls = 0;

 43:   PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
 44:   PetscCall(VecSetSizes(x, PETSC_DECIDE, 1));
 45:   PetscCall(VecSetFromOptions(x));
 46:   PetscCall(VecDuplicate(x, &r));

 48:   PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
 49:   PetscCall(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, 1, 1));
 50:   PetscCall(MatSetFromOptions(J));
 51:   PetscCall(MatSetUp(J));

 53:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
 54:   PetscCall(SNESSetFunction(snes, r, FormFunction, NULL));
 55:   PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, &appctx));
 56:   PetscCall(SNESSetNPCSide(snes, PC_LEFT));
 57:   PetscCall(SNESGetNPC(snes, &npc));
 58:   PetscCall(SNESSetType(npc, SNESNEWTONLS));

 60:   PetscCall(SNESSetFromOptions(snes));
 61:   PetscCall(SNESSetUp(snes));

 63:   PetscCall(VecSet(x, 1.0));
 64:   PetscCall(SNESComputeJacobian(snes, x, J, J));

 66:   PetscCheck(appctx.jac_calls > 0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Jacobian function was not called with left NPC");
 67:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Jacobian called with left NPC\n"));

 69:   PetscCall(VecDestroy(&x));
 70:   PetscCall(VecDestroy(&r));
 71:   PetscCall(MatDestroy(&J));
 72:   PetscCall(SNESDestroy(&snes));
 73:   PetscCall(PetscFinalize());
 74:   return 0;
 75: }

 77: /*TEST

 79:    test:
 80:       suffix: left_npc

 82: TEST*/