Home > source > @penlab > solve.m

# solve

## PURPOSE

main solver loop ("pennon/pbm" iterations)

## SYNOPSIS

function [ifail] = solve(obj)

## DESCRIPTION

     main solver loop ("pennon/pbm" iterations)

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

## SUBFUNCTIONS

- function [b]=stop_crit(f1,f2,eps)
- function [u] = lmlt_update(u,u_updt,phiind,mu,umin)

## SOURCE CODE

    0001 % main solver loop ("pennon/pbm" iterations)
    0002 function [ifail] = solve(obj)
    0003
    0004   ifail = 0;
    0005
    0006   % right phase?
    0007   if (obj.phase<1)
    0008     error('The problem needs to be initialized first.');
    0009   else
    0010     obj.phase=2;  % solving
    0011   end
    0012
    0013   PENALTY_UPDT = obj.allopts.penalty_update;
    0014   PENALTY_UPDT_BAR = obj.allopts.penalty_update_bar;
    0015   MAX_PBMITER = obj.allopts.max_outer_iter;
    0016   PBMALPHA = obj.allopts.outer_stop_limit;
    0017   KKTALPHA = obj.allopts.kkt_stop_limit;
    0018   MU = obj.allopts.mlt_update;
    0019   UINIT = obj.allopts.uinit;
    0020   UINIT_BOX = obj.allopts.uinit_box;
    0021   UINIT_EQ = obj.allopts.uinit_eq;
    0022   UMIN = obj.allopts.umin;
    0023   PINIT = obj.allopts.pinit;
    0024   PINIT_BAR = obj.allopts.pinit_bar;
    0025   USEBARRIER = obj.allopts.usebarrier;
    0026   %visual = obj.allopts.visualisation;
    0027   GAP_NOPROGRESS = 1e20;
    0028   OBJ_UNBOUNDED = -1e20;
    0029   MIN_MAX_FAILED = 3;
    0030
    0031   % clear stats
    0032   obj.clearstats();
    0033
    0034   nFlag=length(penlab.solvermsg);  % last message: not finished yet
    0035   starttime=cputime;
    0036
    0037   %x=x0;  % xall? it gets projected automatically
    0038   nGapNoProgress=0;
    0039   nGrowingGap=0;
    0040   nMinFailed=0;
    0041
    0042   % do we have any constraints? (because of printing)
    0043   bConstrIneq=obj.Nxbox+obj.Nineq+obj.NYbox+obj.NA>0;
    0044   bConstrEq=obj.Neq>0;
    0045
    0046   % print out basic information for better orientation in log files
    0047   obj.print(3,Inf,'Problem name: %s',obj.probname);
    0048   if (~isempty(obj.comment))
    0049     obj.print(3,Inf,'Description:  %s',obj.comment);
    0050   end
    0051   obj.print(3,Inf,'Start time:   %s',datestr(now,0));
    0052   obj.print_opts(3,Inf);
    0053   obj.print(3,Inf,' ');
    0054
    0055   % Initialize penalty parameters and Lagrangian multipliers based on
    0056   % option settings if there are not valid values yet.
    0057   obj.init(false);
    0058
    0059   % automatic modification of the starting point
    0060   % reconstruct from Nxbox constraints lbx and ubx. This is a bit silly
    0061   % because in penlab() we have them, however, the point modification
    0062   % should be done before each solve() (if necessary/required) thus here.
    0063   lbx = -Inf(obj.Nx+obj.NYnnz,1);
    0064   ubx = Inf(obj.Nx+obj.NYnnz,1);
    0065   ind = find(obj.xboxmlt<0);
    0066   lbx(ind) = obj.xboxshift(ind);
    0067   ind = find(obj.xboxmlt>0);
    0068   ubx(ind) = -obj.xboxshift(ind);
    0069   if (obj.allopts.xinit_mod)
    0070     % push x within the box bounds even if not treated by barrier
    0071     % so that they are at least relative 'frac'tion from the boundary
    0072     frac = 0.01;
    0073     dist = frac*min(ubx - lbx,1);
    0074     obj.xall=max(lbx+dist,obj.xall);
    0075     obj.xall=min(ubx-dist,obj.xall);
    0076   elseif (~isempty(obj.xboxindbar))
    0077     % modify only x in barrier not satisfying the constraints
    0078     % or very close to the boundary to set up a valid point
    0079     xboxx = obj.xboxshift + obj.xboxmlt .* obj.xall(obj.xboxmap);
    0080     ind=find(xboxx(obj.xboxindbar)>=-1e-7);
    0081     if (~isempty(ind))
    0082       % from box constraint numbers find variable numbers x involved
    0083       xind=unique(obj.xboxmap(obj.xboxindbar(ind)));
    0084       frac = 1e-5;
    0085       dist = frac*min(ubx(xind) - lbx(xind),1);
    0086       obj.xall(xind)=max(lbx(xind)+dist,obj.xall(xind));
    0087       obj.xall(xind)=min(ubx(xind)-dist,obj.xall(xind));
    0088     end
    0089   end
    0090
    0091   % problem information output
    0092   nNLNineq=sum(obj.ineqmap<=obj.NgNLN);
    0093   nNLNeq=sum(obj.eqmap<=obj.NgNLN);
    0094   nNLNAineq=sum(obj.Amap<=obj.NANLN);
    0095   obj.print(2,Inf,'*******************************************************************************');
    0096   obj.print(2,Inf,penlab.solvername);
    0097   obj.print(2,Inf,'*******************************************************************************');
    0098   obj.print(2,Inf,'Number of variables                      %7d',obj.Nx);
    0099   obj.print(2,Inf,'Number of matrix variables               %7d',obj.NY);
    0100   obj.print(2,Inf,'   - degrees of freedom (var. elements)  %7d',obj.NYnnz);
    0101   obj.print(2,Inf,'(Function) constraints');
    0102   obj.print(2,Inf,'   - box inequalities                    %7d',obj.Nxbox);
    0103   obj.print(2,Inf,'   - linear inequalities                 %7d',obj.Nineq-nNLNineq);
    0104   obj.print(2,Inf,'   - nonlinear inequalities              %7d',nNLNineq);
    0105   obj.print(2,Inf,'   - linear equalities                   %7d',obj.Neq-nNLNeq);
    0106   obj.print(2,Inf,'   - nonlinear equalities                %7d',nNLNeq);
    0107   obj.print(2,Inf,'Matrix constraints');
    0108   obj.print(2,Inf,'   - box inequalities                    %7d',obj.NYbox);
    0109   obj.print(2,Inf,'   - linear inequalities                 %7d',obj.NA-nNLNAineq);
    0110   obj.print(2,Inf,'   - nonlinear inequalities              %7d',nNLNAineq);
    0111   obj.print(2,Inf,' ');
    0112   if (obj.Nxbox>0)
    0113     obj.print(2,Inf,'Min./Max. box-mult.:   %9f / %9f',min(obj.uxbox),max(obj.uxbox));
    0114   end
    0115   if (obj.Nineq>0)
    0116     obj.print(2,Inf,'Min./Max. ineq-mult.:  %9f / %9f',min(obj.uineq),max(obj.uineq));
    0117   end
    0118   if (obj.Neq>0)
    0119     obj.print(2,Inf,'Min./Max. equal-mult.: %9f / %9f',min(obj.ueq),max(obj.ueq));
    0120   end
    0121
    0122   % evaluation of AL will update all objx, xboxx, ineqx, eqx ...
    0123   % and we nned its derivative anyway...
    0124   obj.eval_alx();
    0125   obj.eval_aldx();
    0126   % just in case I use old notation... TODO remove it later and use obj.*
    0127   % directly
    0128   f=obj.objx;
    0129   constrx=[obj.xboxx;obj.ineqx];
    0130   eqx=obj.eqx;
    0131   lagrx=obj.ALx;
    0132   gradx=obj.ALdx;
    0133   gradeqx=obj.eqdx;
    0134   uineq_updt = [];
    0135   uxbox_updt = [];
    0136
    0137   obj.rNormG=norm(obj.ALdx);
    0138   rFeasM=obj.feas_ay();
    0139   obj.rFeas=max([abs(obj.eqx);max(0,obj.xboxx);max(0,obj.ineqx);rFeasM;0]);
    0140   rGap=abs(obj.objx-obj.ALx);
    0141   if (bConstrIneq)
    0142     rPmin=min([obj.pxbox;obj.pineq;obj.PYbox;obj.PA]);
    0143   else
    0144     rPmin=0;
    0145   end
    0146   if (bConstrIneq || bConstrEq)
    0147     %rCompl=max(abs(u.*constrx));
    0148     % TODO ... what if is somewhere barrier?? (and using uxbox??)
    0149     % TODO matrix complementarity??
    0150     obj.rCompl=max([0;abs(obj.uxbox.*obj.xboxx);abs(obj.uineq.*obj.ineqx);abs(obj.eqx.*obj.ueq)]);
    0151       %%% ^^^ all right? but ueq is 0 (if by default) anyway... but might be by DCF
    0152   else
    0153     obj.rCompl=0;
    0154   end
    0155   % detailed starting information
    0156   obj.print(3,Inf,'******************** Start *********************');
    0157   obj.print(3,Inf,'Objective            %27.16E',obj.objx);
    0158   obj.print(3,Inf,'Augmented Lagrangian %27.16E',obj.ALx);
    0159   obj.print(3,Inf,'|f(x) - Lagr(x)|     %27.16E',rGap);
    0160   obj.print(3,Inf,'Grad augm. lagr.     %27.16E',obj.rNormG);
    0161   obj.print(3,Inf,'Feasibility (max)    %27.16E',obj.rFeas);
    0162   obj.print(3,Inf,'Feasibility eqx      %27.16E',max(abs(obj.eqx)));
    0163   obj.print(3,Inf,'Feasibility ineq     %27.16E',max(max(0,obj.ineqx)));
    0164   obj.print(3,Inf,'Feasibility box      %27.16E',max(max(0,obj.xboxx)));
    0165   obj.print(3,Inf,'Feasibility m.ineq   %27.16E',rFeasM);
    0166   obj.print(3,Inf,'Complementarity      %27.16E',obj.rCompl);
    0167   obj.print(3,Inf,'Minimal penalty      %27.16E',rPmin);
    0168   obj.print(3,Inf,'************************************************');
    0169   obj.print(3,Inf,' ');
    0170
    0171   obj.print(2,3,'*******************************************************************************');
    0172   obj.print(2,3,'* it |     obj      | (U,G(x)) |  ||dF||  |   feas   |   pmin   |  Nwt | InIt |');
    0173   obj.print(2,3,'*******************************************************************************');
    0174   if (bConstrIneq || bConstrEq)
    0175     obj.print(2,3,'| %3d|%13.5e |%9.1e |%9.1e |%9.1e |%9.1e | %4d | %4d |',0,obj.objx,obj.rCompl,obj.rNormG,obj.rFeas,rPmin,obj.miter_last,obj.initer_last);
    0176   else
    0177     obj.print(2,3,'| %3d|%13.5e |          |%9.1e |          |          | %4d | %4d |',0,obj.objx,obj.rNormG,obj.miter_last,obj.initer_last);
    0178   end
    0179
    0180   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0181   nFlag=7;
    0182   for pbmiter=1:MAX_PBMITER
    0183
    0184     f_old = obj.objx;
    0185     rGap_old = rGap;
    0186
    0187     obj.print(3,Inf,'************* Start of outer step %3i **********',pbmiter);
    0188
    0189     % TODO tracing information?
    0190
    0191     % unconstr minimization step
    0192     if (obj.Neq==0)
    0193         %fnc.obj=@(xtmp) auglagr(xtmp,u,p,'',ps);
    0194         %fnc.obj_grad=@(xtmp) auglagr_D(xtmp,u,p,'',ps);
    0195         %fnc.obj_hess=@(xtmp) auglagr_D2(xtmp,u,p,'',ps);
    0196
    0197         [nFlagMin,rResults]=obj.unconstr_min();
    0198         lagrx=rResults(1);
    0199         obj.rNormG=rResults(2);
    0200         % or without rResults??
    0201         %lagrx=obj.ALx;
    0202         %obj.rNormG=norm(obj.ALdx);
    0203     else
    0204         %fnc.obj=@(xtmp,ueqtmp) auglagr(xtmp,u,p,ueqtmp,ps);
    0205         %fnc.obj_grad=@(xtmp,ueqtmp) auglagr_D(xtmp,u,p,ueqtmp,ps);
    0206         %fnc.obj_hess=@(xtmp,ueqtmp) auglagr_D2(xtmp,u,p,ueqtmp,ps);
    0207
    0208         [nFlagMin,rResults]=obj.eqconstr_min();
    0209         lagrx=rResults(1);
    0210         obj.rNormG=rResults(2);
    0211         % or directly?
    0212         %lagrx=obj.ALdx
    0213         %obj.rNormG=norm(obj.ALddx);
    0214     end
    0215
    0216     % these should be already evaluated by linesearch ... just to be sure
    0217     % evaluation of AL will update all objx, xboxx, ineqx, eqx ...
    0218     % and we nned its derivative anyway...
    0219     obj.eval_alx();
    0220     obj.eval_aldx();
    0221
    0222     % trial (full) update of Lagrangian multipliers
    0223     % (compute update of the multipliers as if it wasn't restricted by MU)
    0224     if (obj.Nxbox+obj.Nineq~=0)
    0225       %if (USEBARRIER)   % update only multipliers for 'non-barriered' constraints
    0226       %  indbar=[1:ps.N_BOUNDS];
    0227       %  ind=[ps.N_BOUNDS+1:ps.N_CONSTR-ps.N_EQUAL];
    0228       %else
    0229       %  indbar=[];
    0230       %  ind=[1:ps.N_CONSTR-ps.N_EQUAL];
    0231       %end
    0232
    0233       %u_updt(barind) = p(barind).*phibar_D(constrx(barind));
    0234       %u_updt(phiind) = phi2_D(constrx(phiind)./p(phiind));
    0235       uineq_updt = obj.phi2_D(obj.ineqx./obj.pineq);
    0236       uxbox_updt = obj.phi2_D(obj.xboxx./obj.pxbox);
    0237     end
    0238
    0239     % print results
    0240     % just in case I use old notation... TODO remove it later and use obj.*
    0241     % directly
    0242     f=obj.objx;
    0243     constrx=[obj.xboxx;obj.ineqx];
    0244     eqx=obj.eqx;
    0245     lagrx=obj.ALx;
    0246     gradx=obj.ALdx;
    0247     gradeqx=obj.eqdx;
    0248
    0249     obj.rNormG=norm(obj.ALdx);
    0250     rFeasM=obj.feas_ay();
    0251     obj.rFeas=max([abs(obj.eqx);max(0,obj.xboxx);max(0,obj.ineqx);rFeasM;0]);
    0252     rGap=abs(obj.objx-obj.ALx);
    0253     if (bConstrIneq)
    0254       % TODO this hasn't been change, has it?
    0255       rPmin=min([obj.pxbox;obj.pineq;obj.PYbox;obj.PA]);
    0256     end
    0257     if (bConstrIneq || bConstrEq)
    0258       % TODO ... what if is somewhere barrier?? (and using uxbox??)
    0259       % TODO matrix complementarity??
    0260       obj.rCompl=max([0;abs(obj.uxbox.*uxbox_updt.*obj.xboxx);abs(obj.uineq.*uineq_updt.*obj.ineqx);abs(obj.eqx.*obj.ueq)]);
    0261     end
    0262
    0263     if (nFlagMin>0)  % solver failed
    0264       obj.print(3,Inf,'*****!!!!!!! Result of outer step %3i !!!!!*****',pbmiter);
    0265     else
    0266       obj.print(3,Inf,'************ Result of outer step %3i **********',pbmiter);
    0267     end
    0268     obj.print(3,Inf,'Objective            %27.16E',obj.objx);
    0269     obj.print(3,Inf,'Augmented Lagrangian %27.16E',obj.ALx);
    0270     obj.print(3,Inf,'|f(x) - f(x_old)|    %27.16E',abs(f-f_old));
    0271     obj.print(3,Inf,'|f(x) - Lagr(x)|     %27.16E',rGap);
    0272     obj.print(3,Inf,'Grad augm. lagr.     %27.16E',obj.rNormG);
    0273     obj.print(3,Inf,'Feasibility (max)    %27.16E',obj.rFeas);
    0274     obj.print(3,Inf,'Feasibility eqx      %27.16E',max(abs(obj.eqx)));
    0275     obj.print(3,Inf,'Feasibility ineq     %27.16E',max(max(0,obj.ineqx)));
    0276     obj.print(3,Inf,'Feasibility box      %27.16E',max(max(0,obj.xboxx)));
    0277     obj.print(3,Inf,'Feasibility m.ineq   %27.16E',rFeasM);
    0278     obj.print(3,Inf,'Complementarity      %27.16E',obj.rCompl);
    0279     obj.print(3,Inf,'Minimal penalty      %27.16E',rPmin);
    0280     obj.print(3,Inf,'Newton steps                               %5d',obj.miter_last);
    0281     obj.print(3,Inf,'Inner steps                                %5d',obj.initer_last);
    0282     obj.print(3,Inf,'Linesearch steps                           %5d',obj.lsiter_last);
    0283     obj.print(3,Inf,'Time of the minimization step   %14g s',obj.stats_time_miter_last);
    0284     obj.print(3,Inf,'  - factorizations in the step  %14g s',obj.stats_time_fact_last);
    0285     if (nFlagMin>0)   % solver failed
    0286       obj.print(3,Inf,'*****!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*****');
    0287       warn_sign='!';
    0288     else
    0289       obj.print(3,Inf,'************************************************');
    0290       warn_sign=' ';
    0291     end
    0292     obj.print(3,Inf,' ');
    0293
    0294     if (bConstrIneq || bConstrEq)
    0295       obj.print(2,3,'| %3d|%13.5e |%9.1e |%9.1e |%9.1e |%9.1e | %4d | %4d%s|',pbmiter,obj.objx,obj.rCompl,obj.rNormG,obj.rFeas,rPmin,obj.miter_last,obj.initer_last,warn_sign);
    0296     else
    0297       obj.print(2,3,'| %3d|%13.5e |          |%9.1e |          |          | %4d | %4d%s|',pbmiter,obj.objx,obj.rNormG,obj.miter_last,obj.initer_last,warn_sign);
    0298     end
    0299
    0300     % visualisation of VTS problems??
    0301     %if (visual)
    0302     %  vis_vts3('',x);
    0303     %end
    0304
    0305     % stoping criterion
    0306     if (nFlagMin>0)
    0307       nMinFailed=nMinFailed+1;
    0308       obj.print(3,Inf,'Warning: Unconstr/eqconstr. min failed (flag %i), remaing attepmts: %i',nFlagMin,MIN_MAX_FAILED-nMinFailed);
    0309       % perhaps don't update penalty so strongly?
    0310     else
    0311       nMinFailed=0;
    0312     end
    0313     if (nMinFailed >= MIN_MAX_FAILED)
    0314       nFlag=5;
    0315       break;
    0316     end
    0317
    0318     %if (obj.rNormG<PBMALPHA)
    0319     %if (obj.rNormG<PBMALPHA && (isempty(obj.rFeas) || obj.rFeas<PBMALPHA))
    0320     if (stop_crit(obj.objx, f_old, PBMALPHA) && stop_crit(obj.objx, obj.ALx, 10*PBMALPHA) && (obj.rFeas<1e-4*abs(obj.ALx) || obj.rFeas<1e-3))
    0321       if ( (obj.rNormG<KKTALPHA && obj.rCompl<KKTALPHA && obj.rFeas<KKTALPHA) || (obj.rNormG<10*KKTALPHA && obj.rCompl<0.01*KKTALPHA && obj.rFeas<0.01*KKTALPHA) )
    0322         nFlag=1;
    0323         break;
    0324       end
    0325     end
    0326     if (rGap > GAP_NOPROGRESS)
    0327       nGapNoProgress=nGapNoProgress+1;
    0328     else
    0329       nGapNoProgress=0;
    0330     end
    0331     if (rGap > rGap_old)
    0332       nGrowingGap=nGrowingGap+1;
    0333     else
    0334       nGrowingGap=0;
    0335     end
    0336     if (nGapNoProgress>2 || nGrowingGap>40)
    0337       nFlag=6;  % no progress, stop it
    0338       break;
    0339     end
    0340     if (obj.ALx<OBJ_UNBOUNDED)
    0341       nFlag=4;  % objective seems to be unbounded, stop it
    0342       break;
    0343     end
    0344
    0345     % (restricted) multipliers update
    0346     obj.uxbox = lmlt_update(obj.uxbox,uxbox_updt,obj.xboxindphi,MU,UMIN);
    0347     obj.uineq = lmlt_update(obj.uineq,uineq_updt,obj.ineqindphi,MU,UMIN);
    0348     % (restricted) matrix multipliers update
    0349     obj.lmltm_update();
    0350     % penalty update
    0351     if (obj.Nxbox>0)
    0352       obj.pxbox(obj.xboxindbar) = obj.pxbox(obj.xboxindbar).*PENALTY_UPDT_BAR;
    0353       obj.pxbox(obj.xboxindphi) = obj.pxbox(obj.xboxindphi).*PENALTY_UPDT;
    0354     end
    0355     if (obj.Nineq>0)
    0356       %obj.pineq(obj.ineqindbar) = obj.pineq(obj.ineqindbar).*PENALTY_UPDT_BAR;
    0357       obj.pineq(obj.ineqindphi) = obj.pineq(obj.ineqindphi).*PENALTY_UPDT;
    0358     end
    0359     % matrix penalty update
    0360     obj.mpen_update();
    0361
    0362     % force recompute ALx, ALdx, ALddx (point hasn't change so
    0363     % still valid ticker but because penalties and Lagrangian
    0364     % changed, AL*x is actually not valid)
    0365     % this is a bit dirty, ideally P(:) and U(:) should have
    0366     % a counter as well and ALx would be max(xTck, PTck, UTck)
    0367     obj.ALxtck=0;
    0368     obj.ALdxtck=0;
    0369     obj.ALddxtck=0;
    0370
    0371   end    % of pbmiter cycle
    0372   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0373
    0374   % full multipliers update if not terminated by maxit
    0375   %[f,constrx,eqx] = ps.func_val(x); % not necessary, x unchanged
    0376   %if (nFlag~=4 && ps.N_CONSTR-ps.N_EQUAL~=0)
    0377   %  u = u.*u_updt;
    0378   %end
    0379   % (full) multipliers update
    0380   % TODO update it always? even if error occured
    0381   obj.uxbox(obj.xboxindphi) = obj.uxbox(obj.xboxindphi).*uxbox_updt(obj.xboxindphi);
    0382   obj.uineq(obj.ineqindphi) = obj.uineq(obj.ineqindphi).*uineq_updt(obj.ineqindphi);
    0383
    0384   % print results
    0385   obj.rNormG=norm(obj.ALdx);
    0386   obj.rFeas=max([abs(obj.eqx);max(0,obj.xboxx);max(0,obj.ineqx);0]);
    0387   rGap=abs(obj.objx-obj.ALx);
    0388   if (obj.Nxbox + obj.Nineq>0)
    0389     % TODO this hasn't been change, has it?
    0390     rPmin=min([obj.pxbox;obj.pineq]);
    0391   end
    0392   if (obj.Nxbox + obj.Nineq + obj.Neq>0)
    0393     % TODO ... what if is somewhere barrier?? (and using uxbox??)
    0394     obj.rCompl=max([abs(obj.uxbox.*obj.xboxx);abs(obj.uineq.*obj.ineqx);abs(obj.eqx.*obj.ueq)]);
    0395   end
    0396   %  ??? not necessary to recompute?
    0397   % Note: if finished with maxit -> lagrangian multipliers weren't fully updated --> printed slackness may be absolutely wrong; better to leave the earlier computed one
    0398
    0399   if (nFlag==1)
    0400     if (obj.rNormG>1.0e-4*abs(obj.ALx) && obj.rNormG>1.0e-3)
    0401       nFlag=2;
    0402     end
    0403     if(obj.rFeas>1.0e-4*abs(obj.ALx) && obj.rFeas>1.0e-3)
    0404       nFlag=3;
    0405     end
    0406   end
    0407   obj.stats_time_total = cputime - starttime;
    0408
    0409   obj.print(2,Inf,'*******************************************************************************');
    0410   if (nFlag>=1 && nFlag<length(penlab.solvermsg)-1)
    0411     obj.print(2,Inf,'PenLab %s',penlab.solvermsg{nFlag});
    0412   else
    0413     obj.print(2,Inf,'PenLab: Unknow error (%i)',nFlag);
    0414   end
    0415   obj.print(2,Inf,'*******************************************************************************');
    0416   obj.print(2,Inf,'Objective            %27.16E',obj.objx);
    0417   obj.print(3,Inf,'Augmented Lagrangian %27.16E',obj.ALx);
    0418   obj.print(2,Inf,'Relative precision   %27.16E',abs(obj.ALx-obj.objx)/max(1,obj.objx));
    0419   obj.print(2,Inf,'Compl. Slackness     %27.16E',obj.rCompl);
    0420   obj.print(2,Inf,'Grad augm. lagr.     %27.16E',obj.rNormG);
    0421   obj.print(2,Inf,'Feasibility (max)    %27.16E',obj.rFeas);
    0422   obj.print(3,Inf,'Feasibility eqx      %27.16E',max(abs(obj.eqx)));
    0423   obj.print(3,Inf,'Feasibility ineq     %27.16E',max(max(0,obj.ineqx)));
    0424   obj.print(3,Inf,'Feasibility box      %27.16E',max(max(0,obj.xboxx)));
    0425   obj.print(3,Inf,'Minimal penalty      %27.16E',rPmin);
    0426   obj.print(2,Inf,'Newton steps                               %5d',obj.miter);
    0427   obj.print(2,Inf,'Inner steps                                %5d',obj.initer);
    0428   obj.print(2,Inf,'Linesearch steps                           %5d',obj.lsiter);
    0429   obj.print(2,Inf,'Number of evaluations of');
    0430   obj.print(2,Inf,'   - function values                       %5d',obj.stats_ncall_alx);
    0431   obj.print(2,Inf,'   - gradient values                       %5d',obj.stats_ncall_aldx);
    0432   obj.print(2,Inf,'   - hessian values                        %5d',obj.stats_ncall_alddx);
    0433   obj.print(2,Inf,'Time statistics');
    0434   obj.print(2,Inf,'   - total process time         %14g s',obj.stats_time_total);
    0435   obj.print(2,Inf,'   - all minimization steps     %14g s',obj.stats_time_miters);
    0436   obj.print(2,Inf,'   - all factorizations         %14g s',obj.stats_time_fact);
    0437   obj.print(2,Inf,'   - function values evaluation %14g s',obj.stats_time_alx);
    0438   obj.print(2,Inf,'   - gradient values evaluation %14g s',obj.stats_time_aldx);
    0439   obj.print(2,Inf,'   - hessian values evaluation  %14g s',obj.stats_time_alddx);
    0440   obj.print(2,Inf,'*******************************************************************************');
    0441   obj.print(2,Inf,' ');
    0442
    0443   obj.phase=3;   % solver finished
    0444
    0445   ifail = nFlag;
    0446 end
    0447
    0448 %%%%%%%%%%%%%%% Helper routines %%%%%%%%%%%%%%%%%
    0449
    0450 % basic stopping criterion
    0451 function [b]=stop_crit(f1,f2,eps)
    0452
    0453   b = abs(f1 - f2) < eps*(1+.5*(abs(f2)+abs(f1)));
    0454   return
    0455 end
    0456
    0457 % update of lagrangian multipliers with 'mu' limitation and minimum
    0458 % value check
    0459 function [u] = lmlt_update(u,u_updt,phiind,mu,umin)
    0460   if (~isempty(phiind))
    0461     utmp = u(phiind).*u_updt(phiind);
    0462     utmp = min((1/mu)*u(phiind), utmp);
    0463     u(phiind) = max(mu*u(phiind), utmp);
    0464     % check for multiplier minima (UMIN)
    0465     xchange=abs(u)<umin & u~=0;
    0466     u(xchange)=sign(u(xchange))*umin;
    0467   end
    0468 end
    0469