function A=defpos(N,lmin,lmax) 
  // tirage aléatoire d'une matrice définie positive 
  // dont le spectre aléatoire est réel uniforme sur [lmin,lmax]
  A=diag( (lmax-lmin)*rand(N,1)+lmin);
  P=rand(N,N);
  A=P*A*inv(P);
endfunction 

function A=defpos_sym(N,lmin,lmax) 
  // tirage aléatoire d'une matrice définie positive et symetrique 
  // dont le spectre aléatoire est réel uniforme sur [lmin,lmax]
  A=diag((lmax-lmin)*rand(N,1)+lmin);
  [U,H]=hess(rand(N,N));
  A=U*A*U';
endfunction 

function test_rang(ne) 
// une matrice à coeficients aléatoires avec la loi uniforme 
// dans (0,1) est génériquement inversible 
  ne=100;
  for i=1:ne 
    A1=rand(N,N);
    if rank(A1)<>N then pause;
    end; 
  end
endfunction

function y=f(x,A,b)
  y = (1/2)* x'*A*x  + b'*x
endfunction

function y=df(x,A,b)
  y = (A+A')*x/2  + b;
endfunction

// calcul du alpha de alpha-convexité 

function rho_max=rmax(A) 
  // on suppose ici A symétrique 
  alpha = min(abs(spec((A+A')/2)))
  // valeur de Lipschitz du gradient 
  C= norm((A+A')/2,2)
  // contrainte sur le pas de gradient 
  rho_max= 2*alpha/C^2
endfunction 

function [c,xn,fxn]=gradient(x0,f,df,rho,stop,nmax,A,b) 
// une méthode de gradient 
  c=[];
  xn = x0;
  for i=1:nmax 
    xnp1 = xn - rho*df(xn,A,b);
    fxn=f(xn,A,b);
    c=[c,fxn];
    if norm(xnp1-xn) < stop then return;end 
    xn = xnp1;
  end
  mprintf('Arret sur nombre maximum d''itérations %d',nmax) 
endfunction 

// Rendre la fonction gradient plus générique 

function [c,xn,fxn]=gradient_1(x0,f,df,rho,stop,nmax,varargin) 
// une méthode de gradient 
  c=[];
  xn = x0;
  for i=1:nmax 
    xnp1 = xn - rho*df(xn,varargin(:));
    fxn=f(xn,varargin(:));
    c=[c,fxn];
    if norm(xnp1-xn) < stop then return;end 
    xn = xnp1;
  end
  mprintf('Arret sur nombre maximum d''itérations %d',nmax) 
endfunction 

function T=contrainte(N,p) 
// Une matrice aléatoire NxN de rang p
// pour rajouter des contraintes de la forme T*x = 0 
// T est une matrice NxN aléatoire dont le noyau est de dimension N/2 
  T=rand(N,N/2);
  T= T*T'; 
endfunction 

function y=fc(x,A,b,T,eps)
// la contrainte est de la forme T*x = 0 
// on utilise une penalisation 
  y = (1/2)* x'*A*x  + b'*x + (1/eps)*(1/2)*(T*x)'*T*x
endfunction

function y=dfc(x,A,b,T,eps)
    y = (A+A')/2*x  + b + (1/eps)*T*x
endfunction

function T1=contrainte_plein_rang(T) 
// On change les contraintes pour rendres les lignes linéairement inépendantes 
  [W,rk]=rowcomp(T);
  T1= W*T; 
  T1= T1(1:rk,:) ;
endfunction