Home > source > @penlab > penlab.m

# penlab

## PURPOSE

PenLab (Pennon Laboratory for Matlab, previously PennonM)

## SYNOPSIS

This is a script file.

## DESCRIPTION

     PenLab (Pennon Laboratory for Matlab, previously PennonM)
     a NLP-SDP optimization solver

## CROSS-REFERENCE INFORMATION

This function calls:

- Y2vec Copy matrix data (with respect to the given pattern) to xall,
- clearstats Zero all statistics
- disp overloading default disp()
- eqconstr_min equality constrained minimization
- eval_alddx As eval_alx(), evaluate hessian of the Augmented Lagrangian and Jacobian
- eval_aldx As eval_alx(), evaluate gradient of the Augmented Lagrangian and Jacobian
- eval_alx Evaluate Augmented Lagrangian at the current 'xall' point with the current
- feas_ay Compute feasibility of the matrix constraints (A) and matrix variables (Y)
- init Initialize lagrangian multipliers, penalty parameters etc.
- logfile open/close log file based on option settings (in obj.allopts)
- ls_pennon rewritten line_search() @ line_search_els.c @ Pennlp v.2.3 & Pennon v.0.9
- lseq_pen dir/udir ... new directions for x/equality constraints multipliers
- penlab PenLab (Pennon Laboratory for Matlab, previously PennonM)
- phi2 Penalty function for augmented Lagrangian function.
- phi2_D Derivative od penalty function for augmented Lagrangian function.
- phi2_D2 Second derivative of penalty function for augmented Lagrangian function.
- phibar log-penalty barrier function, t can be a scalar/vector/...
- phibar_D log-penalty barrier function, t can be a scalar/vector/...
- phibar_D2 log-penalty barrier function, t can be a scalar/vector/...
- print Publisher (printing manager) for PenLab
- print_opts print all options from the object via obj.print() and add a flag if it is
- setpentype Set index sets assigning constraint numbers to the type of constraint
- solve_chol (Modified) Newton method via Cholesky factorization (with inertia control)
- solvekkt_ldl A direct solver based on LDL() function in Matlab (version 7.6, R2008a++)
- unconstr_min unconstrained minimization
- vec2Y Map 'vectorized' matrix decision variables back to cell array of matrices Y

This function is called by:

- penlab PenLab (Pennon Laboratory for Matlab, previously PennonM)

## SUBFUNCTIONS

- function obj = penlab(penm)
- function delete(obj)
- function set.opts(obj, newopts)
- function set.xall(obj, xall)
- function x=get.x(obj)
- function Y=get.Y(obj)
- function set.xinit(obj,x)
- function set.Yinit(obj,Y)
- function dname = default_probname()
- function dopts = default_opts()

## SOURCE CODE

    0001 % PenLab (Pennon Laboratory for Matlab, previously PennonM)
    0002 % a NLP-SDP optimization solver
    0003 classdef penlab < handle
    0004
    0005   properties   % fully public
    0006     % Problem (object) name with optional comments
    0007     probname = 'No problem loaded';
    0008     comment = '';
    0009
    0010     % user data structure to be passed to callbacks
    0011     userdata = [];
    0012
    0013     % (accepted) user's option settings
    0014     opts = struct();
    0015
    0016     % starting point given by user
    0017     xinit = [];
    0018     Yinit = [];
    0019
    0020   end
    0021
    0022   properties (SetAccess = private)
    0023     % phase of the object, 0 ~ EMPTY, 1 ~ INIT, 2 ~ SOLVING, 3 ~ FINISHED
    0024     phase = 0;
    0025
    0026     % number of variables (in a vector)
    0027     Nx = 0;
    0028     % number of matrix variables
    0029     NY = 0;
    0030     % number of variables in all matrix variables together
    0031     NYnnz = 0;
    0032     % how are variables mapped into matricies Y, just to inform user
    0033     Ymap = [];
    0034
    0035     % number of nonlinear constraint functions
    0036     NgNLN = 0;
    0037     % number of linear constraint functions
    0038     NgLIN = 0;
    0039     % number of nonlinear matrix constraint functions
    0040     NANLN = 0;
    0041     % number of linear matrix constraint functions
    0042     NALIN = 0;
    0043
    0044     % current point (or result), only updated via get.x(), get.Y() as
    0045     % a projection of xall
    0046     x = [];
    0047     Y = [];
    0048
    0049     %% Statistics
    0050     % no of objfun calls ...
    0051     % no of outer iterations (completed)
    0052     % no of inner iterations (completed)
    0053     stats_ncall_alx = 0;
    0054     stats_time_alx = 0;
    0055     stats_ncall_aldx = 0;
    0056     stats_time_aldx = 0;
    0057     stats_ncall_alddx = 0;
    0058     stats_time_alddx = 0;
    0059     stats_time_fact_last = 0;
    0060     stats_time_fact = 0;
    0061     stats_time_miter_last = 0;
    0062     stats_time_miters = 0;
    0063     stats_time_total = 0;
    0064     miter = 0;
    0065     miter_last = 0;
    0066     initer = 0;
    0067     initer_last = 0;
    0068     lsiter = 0;
    0069     lsiter_last = 0;
    0070     % add major iteration counter, inner iteration counter, LS steps
    0071     % time factor, total
    0072
    0073     % status/meassures while solving and when solver finished
    0074     solveflag = [];
    0075     rFeas = [];
    0076     rCompl = [];
    0077     rNormG = [];
    0078
    0079     % or set it up as 'empty routine'?
    0080     % objfun = @(x,Y,userdata) deal([],userdata);
    0081     objfun = [];
    0082     objgrad = [];
    0083     objhess = [];
    0084     confun = [];
    0085     congrad = [];
    0086     conhess = [];
    0087     lagrhess = [];
    0088     mconfun = [];
    0089     mcongrad = [];
    0090     mconhess = [];
    0091     mconlagrhess = [];
    0092
    0093     % all option settings (default + user's)
    0094     % allopts is linked with opts via set.opts() and shouldn't be changed
    0095     % directly
    0096     allopts = penlab.defopts(1);
    0097
    0098     %%%%%%%%%%%%%%%%% MAKE ME PRIVATE LATER %%%%%%%%%%%%%%%%%%%%%%
    0099     % real number of constraints
    0100
    0101     % number of box constraints on X
    0102     Nxbox = 0;
    0103     % mapping: Nxbox (internal numebring) -> Nx bounds
    0104     xboxmap = [];
    0105     xboxmlt = [];
    0106     xboxshift = [];
    0107     % number of matrix box constraints (no of inequal on Y matrix variables)
    0108     NYbox = 0;
    0109     % vectorized variable (xall(Nx+1:Nx+NYnnz) into Y map
    0110     vec2Ymap = [];
    0111     % mapping: no "box" matrix inequality constraint -> k  of appropriate Y{k}
    0112     Yboxmap = [];
    0113     Yboxmlt = [];
    0114     Yboxshift = [];
    0115     % box constraints on the elements of the matrix variables <== among xbox*
    0116     %NYxbox = 0;
    0117     %Yxboxmap = [];
    0118     %Yxboxmlt = [];
    0119     %Yxboxshift = [];
    0120     % equality on the elements...? <== at the moment 2 inequalities
    0121     %NYxboxeq = 0;
    0122     %Yxboxmapeq = [];
    0123     %Yxboxshifteq = [];
    0124
    0125     % function inequalities, nonlin & lin merged (but in this order)
    0126     %NineqNLN = 0;
    0127     %NineqLIN = 0;
    0128     Nineq = 0;
    0129     ineqmap = [];
    0130     ineqmlt = [];
    0131     ineqshift = [];
    0132     % function equalities, nonlin & lin
    0133     Neq = 0;
    0134     eqmap = [];
    0135     eqshift = [];
    0136     % matrix inequalitites, nonlin & lin merged (kept in order)
    0137     NA = 0;
    0138     Amap = [];
    0139     Amlt = [];
    0140     Ashift = [];
    0141     % list of dependent variables on A(k) k=1..NANLN+NALIN, i.e., the original
    0142     % numbering of matrix constraints
    0143     Adep=[];
    0144
    0145     % type of the penalty used for function and matrix inequalitites
    0146     xboxtype = [];
    0147     ineqtype = [];
    0148     Yboxtype = [];
    0149     Atype = [];
    0150
    0151     % type changed into list of indicies (use setpentype())
    0152     xboxindbar = [];
    0153     xboxindphi = [];
    0154     %ineqindbar = [];
    0155     ineqindphi = [];
    0156     Yboxindbar = [];
    0157     Yboxindphi = [];
    0158     %Aindbar = [];
    0159     Aindphi = [];
    0160
    0161
    0162   end
    0163
    0164   properties %(SetAccess = private, GetAccess = private)
    0165
    0166     % time ticker (with every change of xall adds 1)
    0167     ticker = 1;
    0168     % (users) x and Y copy (projection) of xall
    0169     xYtck = 0;
    0170
    0171
    0172     % all variables, current internal point
    0173     xall = [];
    0174
    0175     % Lagrangian multipliers for box inequalities, function ineq. & eq. and new
    0176     % candidates
    0177     uxbox = [];
    0178     uxboxnew = [];
    0179     uineq = [];
    0180     uineqnew = [];
    0181     ueq = [];
    0182     % Lagrangian multipliers for matrix inequalities: box constr. of matrix
    0183     % variables and matrix constraints inequalities
    0184     UYbox = [];
    0185     UYboxnew = [];
    0186     UA = [];
    0187     UAnew = [];
    0188
    0189     % Penalty parameters for all... (except equalities obviously)
    0190     pxbox = [];
    0191     pineq = [];
    0192     PYbox = [];
    0193     PA = [];
    0194
    0195     %%% Various values of transformed user data
    0196     % TODO add ticker!!!
    0197     objx = 0;
    0198     xboxx = [];
    0199     ineqx = [];
    0200     eqx = [];
    0201     % ticker and Jacobian of equalities
    0202     eqdxtck = 0;
    0203     eqdx = [];
    0204
    0205     % ticker of the augmented lagrangian and its value at 'xall' point
    0206     % and value of all function equalitites
    0207     ALxtck = 0;
    0208     ALx = 0;
    0209
    0210     % ticker and the value of the 1st derivative of ALx
    0211     ALdxtck = 0;
    0212     ALdx = [];
    0213
    0214     % ticker and the value of the Hessian matrix of the Augmented Lagrangian
    0215     ALddxtck = 0;
    0216     ALddx = [];
    0217
    0218     % stream for log file output
    0219     fid = -1;
    0220
    0221     % solve_chol & solvekkt_ldl specific
    0222     chol_lmlast = 0;     % the last nonzero diagonal perturbation
    0223     chol_factperm = [];  % nontrivial permutation used in the last attempt
    0224     % linesearch specific
    0225     ls_small_cnt = 0;
    0226     ls_rnu = [];
    0227   end
    0228
    0229   properties (Constant)
    0230     % solver name & version
    0231     solvername = 'PenLab 1.04 (20140125)';
    0232
    0233     % final solver message based on solverflag
    0234     solvermsg = {'converged: optimal solution', ...
    0235        'converged: suboptimal solution', ...
    0236        'finished: solution primal infeasible', ...
    0237        'finished: objective function seems to be unbounded', ...
    0238        'didn''t converge: unconstrained minimization failed', ...
    0239        'didn''t converge: no progress, iteration process terminated', ...
    0240        'didn''t converge: outer iterations reached', ...
    0241        'didn''t converge: unknown error', ...
    0242        'hasn''t finished yet'};
    0243
    0244     % default options
    0245     % struct array: name, {default_value, type, restriction}
    0246     % where  name is the name of the option used in allopts or opts
    0247     %        default_value   (==> defopts(1) ~ default allopts structure)
    0248     %        restriction ... depending on the type, restriction to the accepted values
    0249     %        type: 'O' other, no restriction, no checking (e.g. usr_prn)
    0250     %              'S' string, no restriction
    0251     %              'I', 'R' integer or real within range restriction(1)~restriction(2)
    0252     %              'M' multiple choice - one of the elements in array restriction
    0253     defopts = struct(...
    0254       'outlev', {2, 'I', [0, Inf]}, ...
    0255       'outlev_file', {5, 'I', [0, Inf]}, ...
    0256       'out_filename', {'penm_log.txt', 'S', []}, ...
    0257       'user_prn', {[], 'O', []}, ...
    0258       'maxotiter', {100, 'I', [0, Inf]}, ...
    0259       'maxiniter', {100, 'I', [0, Inf]}, ...
    0260       ... % from pennon.m
    0261       'penalty_update', {0.5, 'R', [0, 1]}, ...       % PENALTY_UPDT
    0262       'penalty_update_bar', {0.3, 'R', [0, 1]}, ...   % PENALTY_UPDT_BAR
    0263       'mpenalty_update', {0.5, 'R', [0, 1]}, ...       %
    0264       'mpenalty_min', {1e-6, 'R', [0, 1]}, ...       %
    0265       'mpenalty_border', {1e-6, 'R', [0, 1]}, ...       %
    0266       'max_outer_iter', {100, 'I', [0, Inf]}, ...       % MAX_PBMITER
    0267       'outer_stop_limit', {1e-6, 'R', [1e-20, 1]}, ...    % PBMALPHA
    0268       'kkt_stop_limit', {1e-4, 'R', [1e-20, 1]}, ...      % KKTALPHA
    0269       'mlt_update', {0.3, 'R', [0, 1]}, ...            % MU
    0270       'mmlt_update', {0.1, 'R', [0, 1]}, ...            % MU2
    0271       'uinit', {1, 'R', [-Inf, Inf]}, ...                  % UINIT
    0272       'uinit_box', {1, 'R', [-Inf, Inf]}, ...              % UINIT_BOX
    0273       'uinit_eq', {0, 'R', [-Inf, Inf]}, ...               % UINIT_EQ
    0274       'umin', {1e-10, 'R', [0, 1]}, ...               % UMIN
    0275       'pinit', {1, 'R', [0, 1]}, ...                  % PINIT
    0276       'pinit_bar', {1, 'R', [0, 1]}, ...              % PINIT_BAR
    0277       'usebarrier', {0, 'M', [0, 1]}, ...             % USEBARRIER
    0278       'xinit_mod', {0, 'M', [0, 1]}, ...             % modification of xinit
    0279       ... % from unconstr_min.m
    0280       'max_inner_iter', {100, 'I', [0, Inf]}, ...     % MAX_MITER
    0281       'inner_stop_limit', {1e-2, 'R', [0, 1]}, ...    % ALPHA
    0282       'unc_dir_stop_limit', {1e-2, 'R', [0, 1]}, ...  % TOL_DIR
    0283       'unc_solver', {0, 'M', [0, 1, 2, 3, 4, 5]}, ... % solver
    0284       'unc_linesearch', {3, 'M', [0, 1, 2, 3]}, ...   % linesearch
    0285       ... % from eqconstr_min.m
    0286       'eq_dir_stop_limit', {1e-2, 'R', [0, 1]}, ...   % TOL_DIR
    0287       'eq_solver', {0, 'M', [0, 1]}, ...              % solver
    0288       'eq_linesearch', {3, 'M', [0, 1, 2, 3]}, ...    % linesearch
    0289       'eq_solver_warn_max', {4, 'I', [0, 10]}, ...    % solver_warn_max
    0290       'ls_short_max', {3, 'I', [0, 10]}, ...          % ls_short_max
    0291       'min_recover_strategy', {0, 'M', [0, 1]}, ...   % recover_strategy
    0292       'min_recover_max', {3, 'I', [0, 10]}, ...       % recover_max
    0293       ... % from phi2.m
    0294       'phi_R', {-0.5, 'R', [-1, 1]}, ...              % R_default
    0295       ... % linesearchs - ls_armijo.m
    0296       'max_ls_iter', {20, 'I', [0, Inf]}, ...   % max tries before LS fails
    0297       'max_lseq_iter', {20, 'I', [0, Inf]}, ...   % same for LS for equality constrained problems
    0298       'armijo_eps', {1e-2, 'R', [0, 1]}, ...  % when is armijo step satisfactory? P(alp) - P(0) <= eps*alp*P'(0)
    0299       ... % solve_simple_chol.m, solkvekkt_ldl.m, solvekkt_lu.m
    0300       'ldl_pivot', {1e-5, 'R', [0, 1]}, ...  % pivot tolerance for delayed pivots in LDL
    0301       'pert_update', {2., 'R', [0, 100]}, ...  % known aka LMUPDATE, multiplier of the lambda-perturbation factor
    0302       'pert_min', {1e-6, 'R', [0, 1]}, ...   % LMLOW, minimal (starting) perturbation
    0303       'pert_try_max', {50, 'I', [0, Inf]}, ... % max number of attempts to successfully perturbate a matrix
    0304       'pert_faster', {1, 'M', [0, 1]}, ...   % use the last known negative curvature vector to determine perturbation
    0305       ... % solve_simple_chol.m
    0306       'chol_ordering', {1, 'M', [0, 1]}, ... % use symamd for sparse matrices before Cholesky factor.?
    0307       ... % solvers
    0308       'luk3_diag', {1., 'R', [0, Inf]} ...  % diagonal of the (2,2)-block in Luksan 3 preconditioner
    0309     );
    0310
    0311
    0312
    0313   end
    0314
    0315   methods
    0316     Y = vec2Y(obj, vec);
    0317     vec = Y2vec(obj, Y, Ydefault);
    0318     [errmsg] = print(obj, minlev, maxlev, msg, varargin);
    0319     [errmsg] = print_opts(obj, minlev, maxlev);
    0320     []=setpentype(obj,lbxbar,ubxbar,lbYbar,ubYbar);
    0321     [] = init(obj, forceupdate);
    0322     [ret] = phi2(obj,t);
    0323     [ret] = phi2_D(obj,t);
    0324     [ret] = phi2_D2(obj,t);
    0325     [ret] = phibar(obj,t);
    0326     [ret] = phibar_D(obj,t);
    0327     [ret] = phibar_D2(obj,t);
    0328     [status] = eval_alx(obj);
    0329     [status] = eval_aldx(obj);
    0330     [status] = eval_alddx(obj);
    0331     [] = clearstats(obj);
    0332     [dir,nFlag]=solve_chol(obj,matrix,rhs);
    0333     [x_dir,ueq_dir,nFlag,pert]=solvekkt_ldl(obj,H,A,rhs1,rhs2);
    0334     [rRelStep, nFlag]=ls_pennon(obj, dir);
    0335     [rRelStep, nFlag]=lseq_pen(obj, dir, udir, miter);
    0336     [nFlag,rResults]=unconstr_min(obj);
    0337     [nFlag,rResults]=eqconstr_min(obj);
    0338     [] = logfile(obj,task);
    0339     disp(obj);
    0340     [mfeas] = feas_ay(obj);
    0341
    0342     % constructor, you need to provide 'penm' structure
    0343     function obj = penlab(penm)
    0344       if (nargin > 0)
    0345         if (isstruct(penm))
    0346
    0347           % fill in nonexistent (optional) fields
    0348           penm=defaultsfiller(penm);
    0349
    0350           % open the filestream if needed
    0351           obj.logfile(1);
    0352
    0353           % first options to set up printing
    0354           %obj.opts=penm.opts;
    0355
    0356           obj.probname=penm.probname;
    0357           obj.comment=penm.comment;
    0358
    0359           obj.userdata=penm.userdata;
    0360
    0361
    0362           if (penm.Nx<0)
    0363             error('ERR: wrong Nx');
    0364           else
    0365             obj.Nx=penm.Nx;
    0366             % equalities not allowed, unconstrained are OK
    0367             [obj.xboxmap, obj.xboxmlt, obj.xboxshift]= ...
    0368                boundschecker(penm.Nx,penm.lbx,penm.ubx,false,true,[]);
    0369             obj.Nxbox = length(obj.xboxmap);
    0370           end
    0371
    0372           if (penm.NY<0)
    0373             error('ERR: wrong NY')
    0374           else
    0375             obj.NY=penm.NY;
    0376             obj.NYnnz=0;
    0377             filterout=[];
    0378
    0379             for k=1:penm.NY
    0380               [mapper, obj.NYnnz]=createY1map(penm.Y{k}, obj.NYnnz);
    0381               obj.vec2Ymap{k}=mapper;
    0382               if (isempty(mapper.xmap))
    0383                 % ignore/exclude from matrix inequality constraints
    0384                 filterout = [filterout, k];
    0385               end
    0386             end
    0387
    0388             % box equalities on elements of Y (in my notation 'Yx' stuff)
    0389             % what to do with different bounds on the "same" elements
    0390             % (e.g. 5<=x_45 but 3<=x_54) ? --> only lower diag is considered
    0391             % get vectorized lower & upper values (+/-Inf for defaults)
    0392             lbYxvec=obj.Y2vec(penm.lbYx,-Inf);
    0393             ubYxvec=obj.Y2vec(penm.ubYx,Inf);
    0394             % map inequalities after 'Nx' box inequalitites, equal on elemets
    0395             % is ok but put as 2 inequal, unconstr is allowed
    0396             [ineqmap, ineqmlt, ineqshift]= ...
    0397                boundschecker(obj.NYnnz,lbYxvec,ubYxvec,-1,true,[]);
    0398             ineqmap=ineqmap+obj.Nx;
    0399             obj.xboxmap = [obj.xboxmap, ineqmap];
    0400             obj.xboxmlt = [obj.xboxmlt, ineqmlt];
    0401             obj.xboxshift = [obj.xboxshift, ineqshift];
    0402             obj.Nxbox = obj.Nxbox + length(ineqmap);
    0403
    0404             % equalities not allowed, empty matrices filtered out
    0405             [obj.Yboxmap, obj.Yboxmlt, obj.Yboxshift]= ...
    0406                boundschecker(penm.NY,penm.lbY,penm.ubY,false,true,filterout);
    0407             obj.NYbox = length(obj.Yboxmap);
    0408
    0409             % how xall (matrix variable) part maps into Y...
    0410             obj.Ymap = obj.vec2Y([obj.Nx+1:obj.Nx+obj.NYnnz]);
    0411           end
    0412
    0413           % any decision variables at all?
    0414           if (obj.Nx+obj.NYnnz==0)
    0415             error('ERR: no decision variables set!');
    0416           end
    0417
    0418           % starting point (unless user offers a better one later)
    0419           obj.xall=zeros(obj.Nx+obj.NYnnz,1);
    0420           if (isfield(penm,'xinit'))
    0421             obj.xinit=penm.xinit;
    0422           end
    0423           if (isfield(penm,'Yinit'))
    0424             obj.Yinit=penm.Yinit;
    0425           end
    0426
    0427           if (penm.NgNLN<0 || penm.NgLIN<0)
    0428             error('Ng*<0 error');
    0429           else
    0430             obj.NgNLN=penm.NgNLN;
    0431             obj.NgLIN=penm.NgLIN;
    0432
    0433             % equalities are allowed, unconstrained not
    0434             [obj.ineqmap, obj.ineqmlt, obj.ineqshift, obj.eqmap, obj.eqshift]= ...
    0435                boundschecker(penm.NgNLN+penm.NgLIN,penm.lbg,penm.ubg,true,false,[]);
    0436             obj.Nineq = length(obj.ineqmap);
    0437             obj.Neq = length(obj.eqmap);
    0438
    0439           end
    0440
    0441           if (penm.NANLN<0 || penm.NALIN<0)
    0442             error('NA*<0 error');
    0443           else
    0444             obj.NANLN=penm.NANLN;
    0445             obj.NALIN=penm.NALIN;
    0446             % neither equalities inor unconstr are allowed
    0447             [obj.Amap, obj.Amlt, obj.Ashift]= ...
    0448                boundschecker(penm.NANLN+penm.NALIN,penm.lbA,penm.ubA,false,false,[]);
    0449             obj.NA = length(obj.Amap);
    0450
    0451
    0452           end
    0453
    0454           % Are callbacks present?
    0455           if (~isfield(penm,'objfun') || ~isfield(penm,'objgrad'))
    0456             % could consider it as feasible point problem but at the moment
    0457             % it's an error
    0458             error('objfun and/or objgrad not defined');
    0459           end
    0460           obj.objfun=penm.objfun;
    0461           obj.objgrad=penm.objgrad;
    0462           if (obj.NgNLN+obj.NgLIN>0)
    0463             if (~isfield(penm,'confun') || ~isfield(penm,'congrad'))
    0464               error('confun and/or congrad not defined and constraint are present');
    0465             end
    0466             obj.confun=penm.confun;
    0467             obj.congrad=penm.congrad;
    0468           end
    0469           if (isfield(penm,'lagrhess'))
    0470             % hessian of the Lagrangian will be used rather hess of obj & constr
    0471             obj.lagrhess=penm.lagrhess;
    0472           else
    0473             if (isfield(penm,'objhess'))
    0474               obj.objhess=penm.objhess;
    0475             else
    0476               error('Neither lagrhess nor objhess is defined.');
    0477             end
    0478             if (obj.NgNLN>0 && ~isfield(penm,'conhess'))
    0479               error('There are nonlinear function constraints and neither lagrhess nor conhess is defined.');
    0480             elseif (obj.NgNLN>0)
    0481               obj.conhess=penm.conhess;
    0482             end
    0483           end
    0484           if (obj.NALIN+obj.NANLN>0)
    0485             if (~isfield(penm,'mconfun') || ~isfield(penm,'mcongrad'))
    0486               error('mconfun and/or mcongrad not defined and matrix constraints present');
    0487             end
    0488             obj.mconfun=penm.mconfun;
    0489             obj.mcongrad=penm.mcongrad;
    0490           end
    0491           if (obj.NANLN>0)
    0492             if (isfield(penm,'mconlagrhess'))
    0493               obj.mconlagrhess=penm.mconlagrhess;
    0494             elseif (isfield(penm,'mconhess'))
    0495               obj.mconhess=penm.mconhess;
    0496             else
    0497               error('Neither mconlagrhess nor mconhess is defined and nonlinear matrix constraints present');
    0498             end
    0499           end
    0500
    0501           % generate Adep - is the current point good for it?
    0502            if (isfield(penm,'Adep'))
    0503                obj.Adep=penm.Adep;
    0504            else
    0505               obj.Adep=cell(obj.NANLN+obj.NALIN,1);
    0506               for kuser=[1:obj.NANLN+obj.NALIN]
    0507                   xh=rand(size(obj.x));
    0508                   if length(obj.Y)>0
    0509                   for iY=1:length(obj.Y), Yh{iY}=rand(size(obj.Y{iY}));end
    0510                   else Yh{1}=[];end
    0511                   [Akx,obj.userdata] = obj.mconfun(xh, Yh, kuser, obj.userdata);
    0512                   list=[];
    0513                   for i=[1:obj.Nx+obj.NYnnz]
    0514                       [Akdx, obj.userdata] = obj.mcongrad(xh,Yh,kuser,i,obj.userdata);
    0515                       if (~isempty(Akdx)&& nnz(Akdx)>0)
    0516                           list=[list,i];
    0517                       end
    0518                   end
    0519                   obj.Adep{kuser}=list;
    0520               end
    0521            end
    0522
    0523           % generate: *type
    0524
    0525
    0526           % initialize penalty type lists for each constraints (bar/phi?)
    0527           obj.setpentype(penm.lbxbar,penm.ubxbar,penm.lbYbar,penm.ubYbar);
    0528
    0529         else
    0530           error('ERR: provide a penm struct!');
    0531         end
    0532       else
    0533         disp('WARNING: Empty object! Call with "penm" structure!');
    0534       end
    0535       obj.phase = 1;
    0536     end
    0537
    0538     % destructor
    0539     function delete(obj)
    0540       % close log file
    0541       obj.logfile(0);
    0542     end
    0543
    0544     % whenever user changes 'opts', check their validity and copy through
    0545     % to allopts; if after INIT phase (=not called from the constructor),
    0546     % report the change to the printer as well (TODO)
    0547     function set.opts(obj, newopts)
    0548       %disp('changing options?');
    0549       % to accept the option it
    0550       %  - must be already present in 'allopts' (==> must be in penlab.defopts())
    0551       %  - consider only new values
    0552       %  - only if they are allowed (in bounds etc.)
    0553
    0554       %newopts
    0555
    0556       % isstruct?? & nonempty?
    0557       if (isempty(newopts) || ~isstruct(newopts) || isempty(fieldnames(newopts)))
    0558         % restoring all to defautls
    0559         % if wasn't empty & not init... to avoid the message while init
    0560         obj.allopts = penlab.defopts(1);
    0561         obj.opts = struct();
    0562         obj.logfile(1);
    0563         disp('Restoring all option settings to their defaults.');
    0564       else
    0565
    0566         names=fieldnames(newopts);
    0567         discard=find(~isfield(penlab.defopts(1),names));
    0568         if (~isempty(discard))
    0569           str = sprintf(' %s',names{discard});
    0570           fprintf('Warning: these options don''t exist and are ignored:%s\n',str);
    0571         end
    0572
    0573         %keep=find(isfield(obj.allopts,names) && (~isfield(obj.opts,names) || getfield(obj.opts,names)~=getfield(newopts,names)));
    0574         keep=find(isfield(penlab.defopts(1),names));
    0575         reopenlog=0;
    0576         for i=keep'    % need a row vector
    0577           defvalue=getfield(penlab.defopts(1),names{i});
    0578           value=getfield(newopts,names{i});
    0579           if (isfield(obj.opts, names{i}))
    0580             origvalue=getfield(obj.opts,names{i});
    0581           else
    0582             origvalue=[];
    0583           end
    0584           oldvalue=getfield(obj.allopts,names{i});
    0585           newvalue=oldvalue;
    0586
    0587           % delete, new or change obj.opts.names{i} ??
    0588           if (isempty(value) && ~isempty(origvalue))
    0589             % remove this one from user options and restore defaults
    0590             obj.opts = rmfield(obj.opts, names{i});
    0591             obj.allopts = setfield(obj.allopts, names{i}, defvalue);
    0592             newvalue=defvalue;
    0593             disp(sprintf('option %s set back to defautls',names{i}));
    0594
    0595           elseif (isempty(origvalue) || origvalue~=value)
    0596             % add or change
    0597             % check validity of the new value
    0598             obj.opts=setfield(obj.opts, names{i}, value);
    0599             obj.allopts=setfield(obj.allopts, names{i}, value);
    0600             newvalue=value;
    0601             %disp(sprintf('new option %s set to %g',names{i},value));
    0602               % universal format to print?? value can be anything
    0603           end
    0604
    0605           % need to reopen the log file?
    0606           if (strcmp(names{i},'out_filename') || ...
    0607             strcmp(names{i},'outlev_file') && oldvalue*newvalue==0)
    0608             reopenlog=1;
    0609           end
    0610         end
    0611         if (reopenlog)
    0612           obj.logfile(1);
    0613         end
    0614       end
    0615     end
    0616
    0617     % Keep an eye on changes of xall (internal variables) and update ticker
    0618     % to set the rest as invalid; this is the only place where ticker should
    0619     % be updated
    0620     function set.xall(obj, xall)
    0621       %disp('...Tick, xall updated');
    0622       % TODO if ticker is maxint, reset all tickers
    0623       obj.ticker=obj.ticker+1;
    0624       obj.xall=xall;
    0625     end
    0626
    0627     % Update x & Y from xall if not already up-to-date
    0628     function x=get.x(obj)
    0629       if (obj.xYtck<obj.ticker)
    0630         %disp('...update of x,Y initialized by x');
    0631         obj.x=obj.xall(1:obj.Nx);
    0632         obj.Y=obj.vec2Y(obj.xall(obj.Nx+1:end));
    0633         obj.xYtck=obj.ticker;
    0634       end
    0635       x=obj.x;
    0636     end
    0637
    0638     % Update x & Y from xall if not already up-to-date
    0639     function Y=get.Y(obj)
    0640       if (obj.xYtck<obj.ticker)
    0641         %disp('...update of x,Y initialized by Y');
    0642         obj.x=obj.xall(1:obj.Nx);
    0643         obj.Y=obj.vec2Y(obj.xall(obj.Nx+1:end));
    0644         obj.xYtck=obj.ticker;
    0645       end
    0646       Y=obj.Y;
    0647     end
    0648
    0649     % changing starting x by user is not allowed during solving
    0650     function set.xinit(obj,x)
    0651       if (obj.phase==2)
    0652         error('Changing x during computation is not allowed.');
    0653       end
    0654
    0655       %%%disp('...changing x starting point');
    0656       if (length(x)~=obj.Nx)
    0657         error('wrong dimension of x');
    0658       end
    0659       obj.xall(1:obj.Nx) = x;
    0660       if (obj.phase==3)
    0661         % throw away all results, TODO prehaps need to do something more?
    0662         obj.phase=1;
    0663       end
    0664     end
    0665
    0666     % changing starting Y (not allowed during solving)
    0667     function set.Yinit(obj,Y)
    0668       if (obj.phase==2)
    0669         error('Changing Y during computation is not allowed.');
    0670       end
    0671
    0672       %%%disp('...changing Y starting point');
    0673       % if Y is not complete, use current xall to fill in gaps
    0674       obj.xall(obj.Nx+1:end)=obj.Y2vec(Y);
    0675       if (obj.phase==3)
    0676         % throw away all results, TODO prehaps need to do something more?
    0677         obj.phase=1;
    0678       end
    0679     end
    0680
    0681
    0682   end
    0683
    0684   methods (Static)
    0685
    0686     % default Problem name
    0687     function dname = default_probname()
    0688       dname = ['Pennon NLP-SDP problem ', datestr(now) ];
    0689     end
    0690
    0691     % default option settings (only values can be changed, new options
    0692     % need to have a default value first)
    0693     function dopts = default_opts()
    0694       dopts = [];
    0695       dopts.outlev = 2;
    0696       dopts.outlev_file = 5;
    0697       dopts.out_filename = 'penm_log.txt';
    0698       dopts.user_prn = [];
    0699       dopts.maxotiter = 100;
    0700       dopts.maxiniter = 100;
    0701       % ...
    0702       % from pennon.m
    0703       dopts.penalty_update = 0.5;       % PENALTY_UPDT
    0704       dopts.penalty_update_bar = 0.3;   % PENALTY_UPDT_BAR
    0705       dopts.max_outer_iter = 100;       % MAX_PBMITER
    0706       dopts.outer_stop_limit = 1e-6;    % PBMALPHA
    0707       dopts.kkt_stop_limit = 1e-4;      % KKTALPHA
    0708       dopts.mlt_update =0.3;            % MU
    0709       dopts.uinit = 1;                  % UINIT
    0710       dopts.uinit_box = 1;              % UINIT_BOX
    0711       dopts.uinit_eq = 0;               % UINIT_EQ
    0712       dopts.umin = 1e-10;               % UMIN
    0713       dopts.pinit = 1;                  % PINIT
    0714       dopts.pinit_bar = 1;              % PINIT_BAR
    0715       dopts.usebarrier = 0;             % USEBARRIER
    0716       % from unconstr_min.m
    0717       dopts.max_inner_iter = 100;       % MAX_MITER
    0718       dopts.inner_stop_limit = 1e-2;    % ALPHA
    0719       dopts.unc_dir_stop_limit = 1e-2;  % TOL_DIR
    0720       dopts.unc_solver = 0;             % solver
    0721       dopts.unc_linesearch = 3;         % linesearch
    0722       % from eqconstr_min.m
    0723       dopts.eq_dir_stop_limit = 1e-2;   % TOL_DIR
    0724       dopts.eq_solver = 0;              % solver
    0725       dopts.eq_linesearch = 3;          % linesearch
    0726       dopts.eq_solver_warn_max = 4;     % solver_warn_max
    0727       dopts.ls_short_max = 3;           % ls_short_max
    0728       dopts.min_recover_strategy = 0;   % recover_strategy
    0729       dopts.min_recover_max = 3;        % recover_max
    0730       % from phi2.m
    0731       dopts.phi_R = -0.5;               % R_default
    0732
    0733       % linesearchs - ls_armijo.m
    0734       dopts.max_ls_iter = 20;   % max tries before LS fails
    0735       dopts.max_lseq_iter = 20;   % same for LS for equality constrained problems
    0736       dopts.armijo_eps = 1e-2;  % when is armijo step satisfactory? P(alp) - P(0) <= eps*alp*P'(0)
    0737
    0738       % solve_simple_chol.m, solkvekkt_ldl.m, solvekkt_lu.m
    0739       dopts.pert_update = 2.;  % known aka LMUPDATE, multiplier of the lambda-perturbation factor
    0740       dopts.pert_min = 1e-6;   % LMLOW, minimal (starting) perturbation
    0741       dopts.pert_try_max = 50; % max number of attempts to successfully perturbate a matrix
    0742       dopts.pert_faster = 1;   % use the last known negative curvature vector to determine perturbation
    0743
    0744       % solve_simple_chol.m
    0745       dopts.chol_ordering = 1; % use symamd for sparse matrices before Cholesky factor.?
    0746
    0747       % solvers
    0748       dopts.luk3_diag = 1;  % diagonal of the (2,2)-block in Luksan 3 preconditioner
    0749       %%% different stuff from previous 'popt' structure %%%
    0750       % equality constrained/unconstrained minimization
    0751       % - eqconstr_min.m, unconstr_min.m
    0752       %dopts.eq_dir_max_prec = 1e-6; % upper limit to the precision of a solution, do not demand a better (absolute) precision than this one; set 0 to turn it off
    0753
    0754       % tracing of one inner loop
    0755       %dopts.trace = 0;   % turn on tracing?
    0756       %dopts.trace_outer_iter = 1; % which outer iteration to trace? (1~the whole first iter ~ from the very beginning)
    0757       %dopts.trace_filename='trace_point.dcf'; % name of the dcf file if used
    0758
    0759     end
    0760
    0761   end
    0762
    0763
    0764 end
    0765