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

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

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

w = w0;

[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;
invKM = VtKtinvX'*diag(diag(D).^(-1)/k)*VtKtinvX;
iter = 1;

while iter < maxiter

  dw = M'*svec(blk,invKM);
  
  if mod(iter,freq) == 0
    fprintf('iter fval sumw minw term1 %g %g %g %g %g\n', iter, sum(log(diag(D))), 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;
  invKM = VtKtinvX'*diag(diag(D).^(-1)/k)*VtKtinvX;
  iter = iter + 1;
end

fval = sum(log(diag(D)));
fprintf(' termination: iter, fval %g, %g\n',  iter, fval)