% This uses multiplicative algorithm to solve p-criterion.
% 
% Sept 26, 2011

function [w fval] = multiplicativep(M,p,w0,freq,tol,maxiter)

n = sqrt(size(M,1));
M = rearrange(M);
blk{1,1} = 's';
blk{1,2} = n;

w = w0;
iter = 1;
Mw = smat(blk,M*w);
[V,D] = svd(Mw);
swd = sum(diag(D).^p);
if p < -1
  fval = swd;
end
pM = V*diag(diag(D).^(p-1)/swd)*V';

while iter < maxiter
  
  dw = M'*svec(blk,pM);  
  
  if mod(iter,freq)==0
    fprintf('iter sumw minw term1 %g %g %g %g\n', iter, sum(w), min(w) , max(dw))
  end
  
  if max(dw) <= (1+tol)
    break
  end
  
  w = w.*dw;
  
  Mw = smat(blk,M*w);
  [V,D] = svd(Mw);
  swd = sum(diag(D).^p);
  if p < -1
    fval = min(fval,swd);
  end
  pM = V*diag(diag(D).^(p-1)/swd)*V';
  
  iter = iter + 1;
end

if p >= -1
  fval = swd;
end

fprintf(' termination: iter, fval %g, %g\n',  iter, fval)