Bilinear form with Operator

Hi,
I have the following issue. Let V be a discrete space and C: V → V a matrix mapping the coefficients of a function u in V to the coefficients Cu of a function in V. Now I want to solve the problem: Seek u in V with int Cu * Cv dx = int f Cv dx for all v in V. An example is given below. How can I do this?

from dolfin import *
mesh = UnitSquareMesh(4,4)
V = FunctionSpace(mesh,‘P’,1)
u = TrialFunction(V)
v = TestFunction(V)
#Create some PETScMatrix C: V → V (does not make so much sense here,
#but exemplifies what I want to do)
a = inner(u,v)*dx
dummy = Constant(0.0)vdx
C = PETScMatrix()
assemble_system(a,dummy,A_tensor=C)
#b = inner(C u,C v)*dx # this is how my mapping should look like
b = inner(u,v)*dx # just a dummy to have an executable code
L = Constant(1.0)vdx
uh = Function(V)
solve(b == L,uh)