load "ilut"
mesh Th=square(10,10);
fespace Qh(Th,P1);
fespace Vh(Th,P2);
real nu=1;
varf mp(p,q)=int2d(Th)(p*q/nu);

matrix Mp=mp(Qh,Qh);
real[int] diagMp(Qh.ndof); diagMp=Mp.diag; diagMp=1/nu*diagMp;

// We are putting this here to workaround a bug in freefem++
real[int] ilutOUT(Vh.ndof+Qh.ndof);

// Uses ILUT for the velocity part and the mass matrix for the pressure part
func real[int] ilutMp(real[int] & xx) {
	real[int] uuin(Vh.ndof);
	real[int] ppin(Qh.ndof);

	real[int] uuout(Vh.ndof);
	real[int] ppout(Qh.ndof);

	[uuin,ppin]=xx;

	for(int k=0;k<Qh.ndof;++k)  {
		ppout[k]=ppin[k]/diagMp[k];
	}
	//ppout=Mp^-1*ppin;
	uuout=applyIlutPrecond(uuin);

	ilutOUT=[uuout,ppin];
	return ilutOUT;
}
//dumptable(cout);
// Uses ILUT for the velocity and I for the pressure
func real[int] ilut(real[int] & xx) {
	ilutOUT=applyIlutPrecond(xx);
	return ilutOUT;
}

// Dummy do-nothing preconditioner
func real[int] dummy(real[int] & xx) {
	ilutOUT=xx;
	return ilutOUT;
}
