//
// algorithme d'Uzawa pour fonctionnelle quadratique
// contrainte linéaire
//
n=7;
//
minv=3; // valeur propre min de A
maxv=5; // valeur propre max de A
D=diag(minv*ones(n,1)+(maxv-minv)*rand(n,1));
P=rand(n,n);
A=P*D*inv(P);
b=rand(n,1);
//
m=2;
C=rand(m,n);
//
// solution exacte
//
M=[A,C';C,zeros(m,m)];bb=[b;zeros(m,1)];
X=inv(M)*bb;
Xex=X(1:n);
//
disp('solution exacte');disp(Xex);
//
// solution Uzawa
//
Niter=150;rho=0.1;alpha=0.1;
// 
X=rand(n,1);
mu=rand(m,1);
//
for i=1:Niter
  X=X-rho*(A*X-b+C'*mu);
  mu=mu+alpha*C*X;
end
disp('solution approchée Uzawa');disp(X)

