Home > utilities > setuptestpoint.m

# setuptestpoint

## PURPOSE

SETUPTESTPOINT sets a point from original Pennon for comparison

## SYNOPSIS

function [alx,aldx,alddx] = setuptestpoint(obj, filename, type)

## DESCRIPTION

     SETUPTESTPOINT sets a point from original Pennon for comparison
     prob ... initialized object with the same problem as it is stored in DCF
     filename ... name (path) to a DCF (Data Container File) generated from
       Pennon suit, so far PenSDP (for linear SDP stored in SDPA) or Pennon/Pennlp
       for Ampl problems
     type ... what's the type (origin) of the DCF? either 'pennlp' or 'pensdp'

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

- alselftest ALSELFTEST is a self test of evaluation of the Augmented Lagrangian and its

## SOURCE CODE

    0001
    0002 function [alx,aldx,alddx] = setuptestpoint(obj, filename, type)
    0003 % SETUPTESTPOINT sets a point from original Pennon for comparison
    0004 % prob ... initialized object with the same problem as it is stored in DCF
    0005 % filename ... name (path) to a DCF (Data Container File) generated from
    0006 %   Pennon suit, so far PenSDP (for linear SDP stored in SDPA) or Pennon/Pennlp
    0007 %   for Ampl problems
    0008 % type ... what's the type (origin) of the DCF? either 'pennlp' or 'pensdp'
    0009
    0010 % This file is a part of PENLAB package distributed under GPLv3 license
    0011 % Copyright (c) 2013 by  J. Fiala, M. Kocvara, M. Stingl
    0012 % Last Modified: 27 Nov 2013
    0013
    0014   % make datatransfer quiet?
    0015   global DCF_SILENT
    0016   DCF_SILENT=1;
    0017
    0018   if (strcmpi(type,'pennlp') || strcmpi(type,'nlp'))
    0019
    0020     % unfortunately, it is impossible to identify content of the dcf file
    0021     % automatically ... so it is expected that all dcf files are in the same
    0022     % format (for now, xUPgHAuf)
    0023
    0024     % read the data and check their dimensions
    0025     x=datatransfer(filename,0);     % dim Nx
    0026     pen_u=datatransfer(filename,1); % dim Nxbox+Nineq+Neq
    0027       % lagr. mlt. for all mixed together (box, NLN, LIN; NLN & LIN include eq)
    0028       % order of transformed inequalities: lb<=g(x)<=ub should be the same
    0029       % as here, i.e., lb-g(x)<=0 followed by g(x)-ub<=0
    0030     pen_p=datatransfer(filename,2); % dtto as pen_u but for penalty parameters
    0031     aldx=datatransfer(filename,3);  % gradient of augm. lagr.
    0032     pen_H=datatransfer(filename,4); % either (1,1) block of KKT or the whole one
    0033     pen_A=datatransfer(filename,5); % Jacobian of equalities if store separately
    0034     ueq=datatransfer(filename,6);   % Lagr. mlt. for equalities
    0035     alx=datatransfer(filename,7);   % value of augm. lagr. without penalization
    0036       % of violation of equalitites (i.e. augm. lagr., not merit function)
    0037
    0038     % check if the 'obj' looks compatible and copy all the data
    0039     if (length(x)~=obj.Nx || obj.NYnnz~=0)
    0040       error('DCF does not fit the number of variables.');
    0041     end
    0042     if (obj.NA~=0)
    0043       error('Problem object is not a NLP (only) problem.');
    0044     end
    0045     if (length(pen_u)~=obj.Nxbox+obj.Nineq+obj.Neq)
    0046       error('Problem object does not seem to be defined by the same NLP model');
    0047     end
    0048     if (length(pen_u)~=length(pen_p))
    0049       error('Arrays U and P in DCF are of different size? Corrupted DCF?');
    0050     end
    0051     if (length(pen_u)~=length(obj.userdata.BEQUAL_M2))
    0052       error('Mapping to separate ineq+eq is incompatible, different problem loaded?');
    0053     end
    0054     if (length(aldx)~=obj.Nx)
    0055       error('DCF looks incompatible (because of aldx).');
    0056     end
    0057     if (length(ueq)~=obj.Neq)
    0058       error('Different number of equalities.');
    0059     end
    0060     % check pen_H, pen_A ?
    0061
    0062     % get only inequality multipliers and penalties
    0063     u=pen_u(~obj.userdata.BEQUAL_M2);
    0064     p=pen_p(~obj.userdata.BEQUAL_M2);
    0065
    0066     % copy everything
    0067     obj.xall=x;
    0068     obj.uxbox=u(1:obj.Nxbox);
    0069     obj.uineq=u(obj.Nxbox+1:end);
    0070     obj.pxbox=p(1:obj.Nxbox);
    0071     obj.pineq=p(obj.Nxbox+1:end);
    0072     obj.ueq=ueq;
    0073
    0074     % transform Hessian from eqltymode=3 (e.i., separate equality Jacobian)
    0075     [pha,phb]=size(pen_H);
    0076     if (obj.Neq~=0 && isempty(pen_A) && pha==obj.Nx+obj.Neq && phb==pha)
    0077       disp('  ---Separating equalities from input data.---');
    0078       disp(' ');
    0079       pen_A = pen_H(1:obj.Nx,(obj.Nx+1):phb);
    0080       pen_H = pen_H(1:obj.Nx,1:obj.Nx);
    0081     %else
    0082     %  warning('strange pen_H??');
    0083     end
    0084
    0085     % Do NOT trust to the dense hessian, only upper half is for sure complete
    0086     if (~issparse(pen_H))
    0087       disp('  ---pen_H is dense, symmetrizing it for sure---');
    0088       disp(' ');
    0089       pen_H_ut=triu(pen_H,1);
    0090       %pen_H = pen_H.*triu(ones(obj.Nx,obj.Nx)) + pen_H_ut';
    0091       pen_H = triu(pen_H) + pen_H_ut';
    0092     end
    0093
    0094     alddx=[pen_H, pen_A; pen_A', sparse(obj.Neq,obj.Neq)];
    0095
    0096   elseif (strcmpi(type,'pensdp') || strcmpi(type,'sdp'))
    0097     % reads/decodes the saved data structure from Pennon (linear) SDP
    0098     % expect "xUPDfgHv" (default)
    0099
    0100     % don't forget that the indices in datatransfer are 0-based!
    0101     x=datatransfer(filename,0);  % dim Nx; x
    0102     u=datatransfer(filename,1);  % dim Ng; lagr. multiplers for normal inequal.
    0103     p=datatransfer(filename,2);  % dim Ng+Na; penalty params for ineq. & matr
    0104       % really for both? Yes
    0105     dense=datatransfer(filename,3); % density flag
    0106     alx=datatransfer(filename,4); % merit function (AL + penalty for equal)
    0107     aldx=datatransfer(filename,5);% dim Nx; gradient of AL
    0108     alddx=datatransfer(filename,6);% dim Nx x Nx, might be sparse; hessian of AL
    0109     % symetrize alddx if dense (perhaps even if sparse? ... Need to test)
    0110     if (~issparse(alddx))
    0111       % only upper triangle is stored
    0112       alddx = alddx+triu(alddx,1)';
    0113     end
    0114     % matrix constraint lagrangian multipliers
    0115     Na = datatransfer(filename,7); % N_MATRIX_CONSTR
    0116     umat=cell(Na,1);
    0117     for k=1:Na
    0118       % either dense ... dim x dim or sparse (dim*(dim+1)/2) x 1
    0119       % TODO careful, datatransfer expects double as the number of elements!
    0120       % and iterating index is automatically int.
    0121       umatk=datatransfer(filename,double(7+k));
    0122       if (size(umatk,1)~=size(umatk,2))
    0123         % Ak is sparse ==> unpack packed triangular format of umat
    0124         % upper triangle, column order (1,2,3,...,dim nnz)
    0125         % dim:  2*nnz=dim^2+dim -> 0=d^2+d-2nnz -> D=1+8nnz, d=(-1+-sqrt(D))/2
    0126         nnz=size(umatk,1);
    0127         dim=(-1+sqrt(1+8*nnz))/2;
    0128
    0129         % get mapping for the packed format
    0130         Mtmp=triu(ones(dim));
    0131         [irow,icol]=find(Mtmp);
    0132         utmp=sparse(irow,icol,umatk);
    0133         umat{k}=full(utmp+triu(utmp,1)');
    0134         %umat{k}(1:10,1:10)
    0135         %umatk(1:20)
    0136       else
    0137         % dense, as it is
    0138         umat{k}=umatk;
    0139       end
    0140     end
    0141
    0142     % check if the 'obj' looks compatible and copy all the data
    0143     if (length(x)~=obj.Nx || obj.NYnnz~=0)
    0144       error('DCF does not fit the number of variables.');
    0145     end
    0146     if (Na~=obj.NA || obj.NANLN~=0)
    0147       error('DCF has different number of matrix constraints.');
    0148     end
    0149     if (obj.Nxbox~=0 || obj.NgNLN~=0 || obj.Neq~=0)
    0150       error('Problem object does not seem to be from SDPA');
    0151     end
    0152     if (length(u)~=obj.Nineq)
    0153       error('DCF has a different number of linear inequalitites.');
    0154     end
    0155
    0156     obj.xall=x;
    0157     obj.uineq=u;  % same order? should be ... only linear ineq.
    0158     obj.pineq=p(1:obj.Nineq);
    0159     for k=1:Na
    0160       obj.UA{k}=umat{k};   % same order? ... should be
    0161     end
    0162     %obj.UA=umat;
    0163     % in Pennon or my Matlab test there used to be 2*p in the computation
    0164     % this version use only PA(k) ==> pass it here 2*
    0165     obj.PA=2*p(obj.Nineq+1:end);
    0166
    0167
    0168   else
    0169     error('Unknown type of DCF file.');
    0170   end
    0171