So I am attempting to use TAO's mixed complementarity solvers and it basically requires a formGradient/Constraints routine and a formJacobian routine.

I am solving a variational inequality via TAO and want the linear operator A to be my jacobian. However, I don't think it was correctly done because my semi smooth linear solver does not converge (but it should!). 

Attached is the code. Line 108 is what TaoSetJacobianRoutine() calls. Basically I want to assemble the linear operator into the matrix petsc_J. Is what I have the correct implementation?

run the code as:

> python MCP_poisson.py 100 100 2

Thanks,

Justin