% This uses multiplicative algorithm to solve optimal design with p-criterion
% for a general coefficient matrix K.
% 
% Jan 12, 2012

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

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

w = w0;
iter = 1;

[V,D] = svd(smat(blk,M*w));
invX = V*diag(diag(D).^(-1))*V';
KtinvX = K'*invX;
W = KtinvX*K;
[V,D] = svd(W);

VtKtinvX = V'*KtinvX;
swd = sum(diag(D).^(-p));
if p < -1
  fval = swd;
end
pM = VtKtinvX'*diag(diag(D).^(-p-1)/swd)*VtKtinvX;

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;
  
  [V,D] = svd(smat(blk,M*w));
  invX = V*diag(diag(D).^(-1))*V';
  KtinvX = K'*invX;
  W = KtinvX*K;
  [V,D] = svd(W);
  
  VtKtinvX = V'*KtinvX;
  swd = sum(diag(D).^(-p));
  if p < -1
    fval = min(swd, fval);
  end
  pM = VtKtinvX'*diag(diag(D).^(-p-1)/swd)*VtKtinvX;
  
  iter = iter + 1;
end

if p >= -1
  fval = swd;
end

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