function [xo, vout, hist] = cycle(xi,df)
% CYCLE is lens optimization routine
% xi is input vector
% df is defect function df(x)
% xo is solution vector
% vout is solution vector at each iteration
% hist is [merit damping rank chg] for each iteration

nv = length(xi);
xo = xi(:);
fmt = '%3d %10.4g %10.3g %10.3g %2d\n';
f = df(xo);
fprev = norm(f);
vout(1,:)=xo';
hist(1,:)=[fprev 0 0 0];
sing = 1e-8;
for k=1:1000
    [A fz] = calculate_derivatives(df,xo);
    [U W V] = svd(A,0);
    d = diag(W);
    T = zeros(size(V));
    dmax = sing*max(d);
    for j=1:nv
        if (d(j)>dmax)
            T(j,j)=1/d(j);
        end
    end
    dv = -V*T'*U'*fz;
    [pmin fval]=fminbnd(@func,0,1,[],df,xo,dv);
    xo = xo+pmin*dv;
    f = df(xo);
    fval = norm(f);
    nr = length(find(d>dmax));
    chg = 100.*pmin*norm(dv)/norm(xo);
    vout(k+1,:)=xo';
    hist(k+1,:)=[fval pmin nr chg];
    fprintf(fmt,k,fval,pmin,chg,nr);
    % This code releases nearly sigular values
    % It may be better to include them from the beginning.
    if (fval<1e-10 | abs(fval-fprev)<1e-10)
        if (nr<nv)
            sing = 0.01*sing;
        else
            break;
        end
    end
    fprev = fval;
    if (pmin<1e-10)
        if (nr<nv)
            sing=0.01*sing;
        else
            break;
        end
    end
end

Published with MATLAB® 7.6