[petsc-users] Making a FD Jacobian

[Mark F. Adams]

I'm trying to make a Jacobian and have this code working:

  ierr = DMCreateMatrix(dm, &A);CHKERRQ(ierr);

  { /* finite difference Jacobian A */
    SNES snes;
    MatStructure flag;

    ierr = SNESCreate(comm, &snes);CHKERRQ(ierr);
    ierr = SNESSetDM(snes,dm);CHKERRQ(ierr);
    ierr = SNESSetFunction(snes,NULL,ComputeFunctionLap,user);CHKERRQ(ierr);
    ierr = SNESSetJacobian(snes,A,A,SNESComputeJacobianDefaultColor,0);CHKERRQ(ierr);

    ierr = SNESComputeJacobian(snes,X,&A,&A,&flag);CHKERRQ(ierr);
    /* ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */
    ierr = SNESDestroy(&snes);CHKERRQ(ierr);
  }

except I need to comment out these line in snesj2.c line ~75:

  ierr = VecEqual(x1,snes->vec_sol,&solvec);CHKERRQ(ierr);
  if (!snes->vec_rhs && solvec) {
    ierr = MatFDColoringSetF(color,F);CHKERRQ(ierr);
  }
  
The problem is that snes->vec_sol is NULL and VecEqual asserts that these are valid arguments.  How should I go about doing this?

Mark