%
% Name:                   nlsolve.m
%                         ------------
% Stand:                  22.6.03/sw
%
% reads the binary
% parameter file:         parameter_nlsolve.mat     
%
%                         this file contains either the original default
%                         parameters 
%                         or the specific user definitions 
% original    
% parameters are:         max_of_count=10
%                         max_of_bigcount=20
%                         fak_min=0.3    (control of speed of iteration)
%                         lambda=1
% **********************************************************************

%   Rene' Lamour 2002-01-16 
%   $Revisions$ $Date$
%
%   calculation of f(x) = 0 with an accuracy of eps with
%
%   f_k(x) = 0 with an accuracy of eps_k in every step and
%   f_k -> f and eps_k -> eps
%
%   the idea is to solve the intermediate steps with a lower accuracy as 
%   required
%
%   Parameter:  nlfile - file for the computation of f(x) and its Jacobian  
%               x0 - input: inital guess
%                    output: solution
%               ep - relative accuracy
%               lin_gauss - method to prepare Jac so that
%                           lin_elim can solve Jac*x=right
%               lin_elim - solve Jac*x=right (see lin_gauss)
%               print_level - controls the output: 
%                             0 no output
%                             1 standardoutput
%                             2 output to a file named: nlsolve_'date'.log
%               varargin - parameter of nlfile
% **********************************************************************


function [x0,ok] = nlsolve(nlfile,x0,ep,lin_gauss,lin_elim,print_level,varargin);

% read the actual user parameter version from the binary file
% and log the actual status
load parameter_nlsolve;

% ****************************************************
% default parameter definitions are:
% max_of_count=10;   
% max_of_bigcount=20; 
% fak_min=0.3;   (control of speed of iteration)
% lambda=1;
%
% ... will be loaded from parameter_nlsove.mat
% ****************************************************

%   output file
%
fid=file_open_date(print_level,'nlsolve');
fprintf(fid,'%s \n',headline(nlfile));
%
%   initial values
%
fprintf(fid,'%s \n',['relative accuracy = ep = ',num2str(ep)]);
print_mat(fid,'initial values',x0');
%
f_calls=0;
Jac_calls=0;
ok=1;
ep1=1/ep^(1/3);
normdyx=ep+1;
normright=ep+1;
bigcount=0;
x0_prev=[];
dyx_prev=[];
while  normright > ep %normdyx > ep
    % right side
    newJac=1;
    [ok_f,right,Jac]=feval(nlfile,x0,newJac,print_level,varargin{:});
    f_calls=f_calls+1;
    Jac_calls=Jac_calls+1;
    %
    if ok_f
        %
        %   current relative accuracy 
        %
        ep1=sqrt(ep*ep1);
        %
        normright=norm(full(right));
        % 
        fprintf(fid,'\n%s\n\n',['bigcout = ',num2str(bigcount),...
                '   relative accuracy of this step = ',num2str(ep1), '  norm(right) = ',num2str(normright)]);
        %
        if normright < ep   %  eps
            fprintf(fid,'%s \n','nlsove: normright < ep');
            break;
        end
        normright1=normright;
        normright1_old=normright1;
        normfirst=normright1;
        count=0;
        [L,U,sing]=feval(lin_gauss,Jac);
        if sing
            [ok_s,x0,lambda]=nlsolve_sing('fall singulaerer matrix',bigcount,lambda,x0_prev,dyx_prev);
            if ~ok_s
                ok=0;
                fclose(fid);
                return
            end;
        else
            x0_prev=x0; 
            lambda=min(2*lambda,1);
        end;
        lcopy=1;
        while   normright1>ep1 %normdyx > ep1 
            dyx=feval(lin_elim,L,U,right);
            if lcopy
                dyx_prev=dyx;
                lcopy=0;
            end;
            normdyx=rel_error_of_dx(dyx,x0);       
            if  normdyx < eps
                fprintf(fid,'%s \n',['very small correction = ',num2str(normdyx),' < eps']);
                ok=1;
                break;
            end;
            %%%%%%%  damping of Newton method?
            x0=x0-dyx;
            % right side
            newJac=0;
            [ok_f,right]=feval(nlfile,x0,newJac,print_level,varargin{:});
            f_calls=f_calls+1;
            normright1=norm(full(right));
            if normright1_old < normright1
                fprintf(fid,'%s \n','nlsolve: increasing defect');
                if normright1 <  normfirst*100   % factor ?
                    break;
                else
                    ok=0;
                    fprintf(fid,'%s \n','nlsolve: no convergence - choose better initial guess');
                    fclose(fid);
                    return
                end
            else
                fak=normright1/normright1_old;  % rate of convergence
                normright1_old = normright1; 
                if fak > fak_min
                    fprintf(fid,'%s \n',['rate of convergence = ',num2str(fak)]);
                    break;
                end;
            end
            %%%%%%%%%%
            count=count+1;
            fprintf(fid,'%s \n',['count = ',num2str(count),'  norm(right)= ',num2str(normright1),...
                    ' norm(dx) = ',num2str(normdyx)]);
            if count> max_of_count
                fprintf(fid,'%s \n','nlsolve: to many counts');
                break;
            end
        end
        bigcount=bigcount+1;
        if bigcount> max_of_bigcount
            fprintf(fid,'%s \n','nlsolve: to many bigcounts');
            ok=0;
            fclose(fid);
            return;
        end
    else
        [ok_s,x0,lambda]=nlsolve_sing('f(x) or Jac not computable',bigcount,lambda,x0_prev,dyx_prev,x0);
        if ~ok_s
            ok=0;
            fclose(fid);
            return
        end;
    end; 
end
fprintf(fid,'\n%s\n','solution');
print_mat(fid,'x =',x0');
fprintf(fid,'\n%s \n',['norm(dyx)    = ',num2str(normdyx)]);
fprintf(fid,'%s \n',['norm(right)  = ',num2str(normright)]);
fprintf(fid,'%s \n',['bigcount     = ',num2str(bigcount)]);
fprintf(fid,'%s \n',['nlfile_calls = ',num2str(f_calls)]);
fprintf(fid,'%s \n',['Jac_calls    = ',num2str(Jac_calls)]);
if print_level==2
    fclose(fid);
end;

% logfile for nlsolve

nlfile=['nlsolve_',date,'.log'];
log_parameter_nlsolve(nlfile);

