BlockIP Class Reference

Main class for loading and solving problems. More...

#include <BlockIP.h>

Classes
struct	BackupLnk

Public Types
enum	TYPE_PROBLEM { LINEAR, QUADRATIC, NONLINEAR }
enum	TYPE_COMP_DY { CHOL_PWRS_PCG, FULL_CHOL, HYBRID_PCG, CHOL_THETA0_PCG }
enum	TYPE_DIRECTION { NEWTON, SECOND_ORDER, AUTOMATIC }
enum	TYPE_START_POINT { SIMPLE, QUAD_EQCONS_PROB }
enum	TYPE_REG { NO_REG, QUAD_REG, PROX_REG }
enum	OUTPUT { NONE, SCREEN, FILE, BOTH }

Public Member Functions
	BlockIP ()
	Constructor.
	BlockIP (TYPE_PROBLEM type_objective, double *&cost, double *&qcost, void(*fobj)(int n, double x[], double &fx, double Gx[], double Hx[], void *params), void *params, double *&ub, double *&rhs, int numBlocks, bool sameN, MatrixBlockIP N[], bool sameL, MatrixBlockIP L[])
	Create an instance with problem in standard form: 0 <= variables <= ub and constraints = rhs.
	BlockIP (TYPE_PROBLEM type_objective, double *&cost, double *&qcost, void(*fobj)(int n, double x[], double &fx, double Gx[], double Hx[], void *params), void *params, double *&lb, double *&ub, double *&lhs, double *&rhs, int numBlocks, bool sameN, MatrixBlockIP N[], bool sameL, MatrixBlockIP L[])
	Create an instance with general problem: lb <= variables <= ub and lhs <= constraints <= rhs.
	~BlockIP ()
	Destructor.
void	initialize ()
	Initializes class attributes.
void	free_memory ()
	Free all the memory.
void	create_problem (TYPE_PROBLEM type_objective, double *&cost, double *&qcost, void(*fobj)(int n, double x[], double &fx, double Gx[], double Hx[], void *params), void *params, double *&ub, double *&rhs, int numBlocks, bool sameN, MatrixBlockIP N[], bool sameL, MatrixBlockIP L[])
	Create a problem in standard form: 0 <= variables <= ub and constraints = rhs.
void	create_problem (TYPE_PROBLEM type_objective, double *&cost, double *&qcost, void(*fobj)(int n, double x[], double &fx, double Gx[], double Hx[], void *params), void *params, double *&lb, double *&ub, double *&lhs, double *&rhs, int numBlocks, bool sameN, MatrixBlockIP N[], bool sameL, MatrixBlockIP L[])
	Create a general problem: lb <= variables <= ub and lhs <= constraints <= rhs.
void	convert_to_standard ()
	Convert any problem to standard form, restrictions= rhs, variables >= 0, <= ub.
void	check_input_problem_and_compute_dimensions ()
	Check the input problem and compute dimensions.
void	check_sameN ()
	Check if the problems satisfies the conditions when using sameN.
int	minimize ()
void	set_defaults_minimize ()
int	get_it ()
	Return number of IP iterations of minimization.
int	get_cgit ()
	Return number of overall PCG iterations.
double	get_fobj ()
	Return objective function.
double	get_dobj ()
	Return dual objective function.
void	set_inf (double inf=INF)
	Set threshold value to be used as infinity (x>inf -> x is inf)
double	get_inf ()
	Get infinity value to be considered.
void	set_m_pw_prec (int m_pw_prec=M_PW_PREC)
	Set number of terms used as preconditioner of the power series expansion of (D-C'*B^{-1}*B)^{-1}.
int	get_m_pw_prec ()
	Get number of terms used as preconditioner of the power series expansion of (D-C'*B^{-1}*B)^{-1}.
void	set_sigma (double sigma=SIGMA)
	Set reduction of the centrality parameter at each IP iteration.
double	get_sigma ()
	Get reduction of the centrality parameter at each IP iteration.
void	set_rho (double rho=RHO)
	Set reduction of the step-length for the primal and dual variables at each IP iteration.
double	get_rho ()
	Get reduction of the step-length for the primal and dual variables at each IP iteration.
void	set_optim_gap (double optim_gap=OPTIM_GAP)
	Set optimality gap tolerance.
double	get_optim_gap ()
	Get optimality gap tolerance.
void	set_optim_pfeas (double optim_pfeas=OPTIM_PFEAS)
	Set primal feasibility tolerance.
double	get_optim_pfeas ()
	Get primal feasibility tolerance.
void	set_optim_dfeas (double optim_dfeas=OPTIM_DFEAS)
	Set dual feasibility tolerance.
double	get_optim_dfeas ()
	Get dual feasibility tolerance.
void	set_output_freq (int output_freq=OUTPUT_FREQ)
	Set output information lines will be printed each output_freq IP iterations.
int	get_output_freq ()
	Get output information lines will be printed each output_freq IP iterations.
void	set_output (OUTPUT output=SCREEN)
	Set type of output.
OUTPUT	get_output ()
	Get type of output.
void	set_maxiter (int maxiter=MAXITER)
	Set maximum number of IP iterations.
int	get_maxiter ()
	Get maximum number of IP iterations.
void	set_init_pcgtol (double init_pcgtol)
	Set initial tolerance for the conjugate gradient (by default PCG_TOL_LIN if linear problem, PCG_TOL_QUAD if quadratic or -convex- nonlinear)
double	get_init_pcgtol ()
	Get initial tolerance for the conjugate gradient (by default PCG_TOL_LIN if linear problem, PCG_TOL_QUAD if quadratic or -convex- nonlinear)
void	set_min_pcgtol (double min_pcgtol=MIN_PCGTOL)
	Set minimum tolerance for the conjugate gradient.
double	get_min_pcgtol ()
	Get minimum tolerance for the conjugate gradient.
void	set_maxit_pcg (double maxit_pcg=0)
	Set maximum number of pcg iterations (if 0, then the value will be computed by the code)
double	get_maxit_pcg ()
	Get maximum number of pcg iterations (if 0, then the value will be computed by the code)
void	set_red_pcgtol (double red_pcgtol=RED_PCGTOL)
	Set reduction of pcg_tol at each IP iteration.
double	get_red_pcgtol ()
	Get reduction of pcg_tol at each IP iteration.
void	set_show_specrad (bool show_specrad=false)
	Set show_specrad at each IP iteration.
bool	get_show_specrad ()
	Get show_specrad.
void	set_type_comp_dy (TYPE_COMP_DY type_comp_dy=TYPE_COMP_DY_DEFAULT)
	Set how dy direction is computed.
TYPE_COMP_DY	get_type_comp_dy ()
	Get how dy direction is computed.
void	set_type_direction (TYPE_DIRECTION type_direction=TYPE_DIRECTION_DEFAULT)
	Set which direction is computed.
TYPE_DIRECTION	get_type_direction ()
	Get which direction is computed.
void	set_type_start_point (TYPE_START_POINT type_start_point=TYPE_START_POINT_DEFAULT)
	Set how starting point is computed.
TYPE_START_POINT	get_type_start_point ()
	Get how starting point is computed.
void	set_type_reg (TYPE_REG type_reg=TYPE_REG_DEFAULT)
	Set type of regularization.
TYPE_REG	get_type_reg ()
	Get type of regularization.
void	set_factor_reg (double factor_reg=FACTOR_REG_DEFAULT)
	Set factor of regularization.
double	get_factor_reg ()
	Get factor of regularization.
void	set_deactivateLnk (bool deactivateLnk=DEACTIVATELNK)
	Set flag on deactivation of linking.
bool	get_deactivateLnk ()
	Get flag on deactivation of linking.
int	get_k_blocks ()
	Get the number of blocks.
int	get_km ()
	Get the number of constraints without linking constraints.
int	get_kn ()
	Get the number of variables without linking constraints.
int	get_l_link ()
	Get the number of linking constraints.
int	get_n_vars ()
	Get the number of variables of the problem to be optimized.
int	get_m_cons ()
	Get the number of constraints of the problem to be optimized.
int	get_num_vars ()
	Get the number of variables of the original problem.
int	get_num_cons ()
	Get the number of constraints of the original problem.
void	set_whoperm (WHO_PERMUTES whoperm_=CHOLESKY)
	Set whopermutes constraints-related information: the Cholesky factorization of the user of it.
WHO_PERMUTES	get_whoperm ()
	Get whopermutes constraints-related information: the Cholesky factorization of the user of it.
const double *	get_x ()
	Get primal variables of problem optimized.
const double *	get_y ()
	Get dual variables of problem optimized (of equality constraints)
const double *	get_z ()
	Get dual variable of problem optimized (of x>=0) i=1..n_vars.
const double *	get_w ()
	Get dual variable of problem optimized (of x<=u) i=1..n_vars.
double	get_x (int i)
	Get value of primal variable x(i) of problem optimized i=1..n_vars.
double	get_y (int i)
	Get value of dual variable y(i) of problem optimized (of equality constraints) i=1..m_cons.
double	get_z (int i)
	Get value of dual variable z(i) (of x>=0) of problem optimized i=1..n_vars.
double	get_w (int i)
	Get value of dual variable w(i) (of x<=u) of problem optimized i=1..n_vars.
int	min_PCG_Hv (int nn, double *v, double *Hv)
int	min_PCG_Hz_eq_r (int nn, double *zz, double *rr)
int	min_PCG_Theta0z_eq_r (int nn, double *zz, double *rr)
void	set_constant_fobj (double constant)
	Set the constant to add to the objective function.
double	get_constant_fobj ()
	Get the constant to add to the objective function.
void	set_names (string *blockNames, string *varNames, string *consNames, bool copy_vectors=true)
	Set the names of the problem TODO change copy_vectors.
void	get_names (const string *&blockNames, const string *&varNames, const string *&consNames)
	Get the names of the problem.
void	convert_to_std_and_write_mps (const char *filename)
	Convert the problem to standard form (if it is not already converted) and write a mps file.
void	write_mps (const char *filename)
	Write a mps file.
void	write_mps (const char *filename, TYPE_PROBLEM type_objective, double cost[], double qcost[], double lb[], double ub[], double lhs[], double rhs[], int numBlocks, bool sameN, MatrixBlockIP N[], bool sameL, MatrixBlockIP L[], string *blockNames=NULL, string *varNames=NULL, string *consNames=NULL, double constantFObj=0)
	Write a mps file.
void	read_mps (const char *filename)
	Create a problem from a mps file.
void	write_BlockIP_format (const char *filename)
	Write the problem in BlockIP format file.
void	read_BlockIP_format (const char *filename)
	Create a problem from a BlockIP format file.
void	write_problem (const char *filename)
	Write all data related to the problem into a file.
StdForm *	get_std_form (bool keepStdFormAfterDelete=false)
	Get the StdForm related to the problem.
void	create_names ()
	Create names for blocks, variables and constraints if does not exist.
string *	get_full_var_names ()
string *	get_full_cons_names ()

Public Attributes
int	k_blocks
	Number of diagonal blocks.
int	km
	Number of constraints without linking constraints.
int	kn
	Number of variables without slacks of linking.
int	m_cons
	Number of constraints including linking constraints.
int	n_vars
	Number of variables including slacks.
int	l_link
	Number of linking constraints.
int *	ini_n
	Begin to blocks 1..k and linking for variables.
int *	ini_m
	Begin to blocks 1..k and linking for constraints.
double *	cost
	Linear cost of variables including slacks.
double *	qcost
	Quadratic cost of variables including slacks.
double *	lb
	Lower bounds, including slack bounds.
double *	ub
	Upper bounds, including slack bounds.
double *	lhs
	Lower constraint limits, including linking constraints limits.
double *	rhs
	Upper constraint limits, including linking constraints limits.
bool	sameN
	Define if the same matrix is used for each N block.
bool	sameL
	Define if the same matrix is used for each L block.
MatrixBlockIP *	N
	Diagonal blocks.
MatrixBlockIP *	L
	Linking constraints blocks.
MatrixBlockIP *	A
	Full matrix (likely made from N and L)
MatrixBlockIP *	D
	D = Theta_(k+1) + sum{i..k} L_i*Theta_i*L_i'.
MatrixBlockIP *	Theta0
void(*	fobj )(int n, double x[], double &fx, double Gx[], double Hx[], void *params)
	User function to calculate the objective function in a point.
void *	params
	User parameters to perform objective function calculations.
double	fx
double *	Gx
double *	Hx
	Objective function value, Gradient and Hessian in the actual point.
double	dobj
	dual objective
bool	inStdForm
	To control if the problem is in standard form.
string *	blockNames
	Block names of N blocks.
string *	varNames
	Variable names including slacks.
string *	consNames
	Constraint names including linking constraints.
TYPE_PROBLEM	type_objective
	Type of objective function, linear, quadratic or non-linear.
StdForm *	stdForm
	Contains all the information to perform conversions between the original and standard problem.
double	constantFObj
	Constant to add in the objective function.
int *	origNVars
	Number of original variables for each block.
int *	origNCons
	Number of original constraints for each block.

Private Member Functions
void	min_initializations ()
	Initialize parameters, tolerances, etc for BlockIP minimization algorithm.
void	min_preprocess ()
void	min_normb_normc ()
void	min_Ax (double *Ax, const double *x)
void	min_Aty (double *Aty, const double *y)
void	min_compute_s ()
void	min_compute_mu ()
void	min_compute_sigma_psi ()
void	min_update_inactivelnk ()
void	min_KKT_residuals ()
void	min_starting_point ()
void	min_objective_functions_nonlin ()
void	min_objective_functions_linquad ()
void	eval_fobj (double x[], double &fx, double Gx[], double Hx[])
void	min_write_problem_information ()
void	min_write_current_point_info ()
void	min_steplengths (double &talpha_p, double &talpha_d, double tdx[], double tdz[], double tdw[])
bool	min_optimal ()
void	min_compute_Theta ()
void	min_compute_Theta0 ()
void	min_compute_direction ()
void	min_Newton_direction ()
void	min_second_order_predictor_corrector_direction ()
void	min_update_newpoint ()
void	min_test_newpoint ()
void	min_backup_inactive (int lnk)
void	min_restore_inactive (int lnk)
void	min_compute_dy_cholpcg (double *dy, double *rhsdy, bool only_solve=false)
void	min_compute_dy_fullchol (double *dy, double *rhsdy, bool only_solve=false)
void	min_free_memory ()
void	min_check_Newton_direction_is_correct ()
void	min_check_predictor_direction_is_correct ()
void	min_check_predictor_corrector_direction_is_correct ()
void	min_solve_NThetaNt (double *v)
void	min_Ctv (double *v, double *Ctv)
void	min_Cv (double *v, double *Cv)
void	min_debug_write_variables ()
void	min_debug_variables_of_inactivelnk ()
void	min_debug_variables_without_ub ()
void	min_update_qreg ()

Private Attributes
int	input_problem
bool	initialized
	To control if the attributes have been initialized.
bool	keepStdFormAfterDelete
	To keep stdForm when BlockIP will delete.
bool	deleteMatrices
	To delete matrices when created by BlockIP (read_mps)
WHO_PERMUTES	whoperm
double	inf
	Infinity value.
double	optim_gap
	Optimality gap tolerance.
double	optim_pfeas
	Primal feasibility tolerance.
double	optim_dfeas
	Dual feasibility tolerance.
int	maxiter
	Maximum number of IP iterations.
int	output_freq
	Output information lines will be printed each output_freq IP iterations.
OUTPUT	output
	Type of output.
double	epsmach
	Stores computed machine epsilon.
double	init_pcgtol
	Initial tolerance for the conjugate gradient (by default PCG_TOL_LIN if linear problem, PCG_TOL_QUAD if quadratic or -convex- nonlinear)
double	pcgtol
	Tolerance for conjugate gradient at current iteration.
double	min_pcgtol
	Minimum tolerance for the conjugate gradient.
double	maxit_pcg
	Maximum number of pcg iterations (if 0, then the value will be computed by the code)
double	red_pcgtol
	Reduction of pcg_tol at each IP iteration.
int	m_pw_prec
	Number of terms of the power series expansion of (D-C'*B^{-1}*B)^{-1} used as preconditioner.
int	totcgit
	Overall number of CG iterations.
int	cgit
	number of CG iterations of this IPM iteration
int	it
	IPM iterations.
double	est_specrad
	Estimation of spectral radius of D^{-1}C'*B^{-1}*B (computed through Ritz values)
bool	show_specrad
bool	comp_specrad
	Whether the estimation of spectral radius has to be shown in the output.
TYPE_COMP_DY	type_comp_dy
	Whether the estimation of spectral radius has to be computed.
TYPE_DIRECTION	type_direction
	Type of direction (Newton, second order...) given by user.
TYPE_DIRECTION	this_type_direction
	Type of direction (Newton or second order) of this particular iteration.
TYPE_START_POINT	type_start_point
TYPE_REG	type_reg
	How starting point will be computed.
double	factor_reg
	Type of regularization performed.
double *	x
	initial value for regularization
double *	y
double *	z
double *	w
double *	s
	primal and dual optimization variables
double *	rb
double *	rc
double *	rxz
double *	rsw
	residuals of KKT conditions
double *	dxdz
double *	dsdw
	for rhs of corrector system, from dx,dw,ds and dw of predictor direction
double	normb
double	normc
double	normrb
double	normrc
	norms of rhs, costs, rb, and rc
double	alpha_p
double	alpha_d
	step lengths
double	rho
	Reduction of the step-length for the primal and dual variables at each IP iteration.
double *	dx
double *	dy
double *	dz
double *	dw
	Primal and dual optimization directions.
double	gap
	Duality gap.
double	mu
double	mu0
	Centrality parameter; mu0 is value of mu at starting point.
double	sigma
	Reduction of the centrality parameter at each IP iteration.
double	psi
	Reduction of dx*dz vector in the rhs of corrector system (in predictor-corrector)
double *	Theta
	Theta diagonal matrix.
double *	Ax
double *	Aty
double *	r
double *	rhsdy
double *	rhspcg
double *	Cvaux
double *	Cvaux2
	Auxiliary vectors for interior-point algorithm.
double *	r_cg
double *	z_cg
double *	p_cg
	Auxiliary vectors for PCG.
double	qreg
	Regularization term.
double *	valpha
double *	vbeta
bool	deactivateLnk
	Vectors to store information from PCG to later compute Ritz values.
int	numActiveLnk
	Number of active lnk.
int	numInactiveLnk
	Number of inactive lnk.
int	numInactiveLnkWithUb
	Number of inactive lnk with upper bound.
bool *	isActiveLnk
	linking i is active if isActiveLnk[i] is true
int *	listActiveLnk
	list of active linking, j=listActiveLnk[i], i=1..numActiveLnk, j in {1,..,l_link})
int *	listInactiveLnk
	list of inactive linking, j=listInactiveLnk[i], i=1..numInactiveLnk, j in {1,..,l_link})
int	numFreeVars
	Number of variables marked as free.
int	numUnfreeVars
	Number of variables not marked as free.
bool *	isFreeVar
	Variable i is marked as free if isFreeVar[i] is true.
int *	listFreeVars
	List of variables marked as free, j=listFreeVars[i], i=1..numFreeVars, j in {1,..,n_vars})
int *	listUnfreeVars
	List of variables not marked as free, j=listUnfreeVars[i], i=1..numUnfreeVars, j in {1,..,n_vars})
int	numBackupLnk
BackupLnk *	backupLnk
int	numWithUb
	Number of variables with upper bound (not infinity)
int	numWithoutUb
int *	listWithoutUb
int *	listWithUb

Static Private Attributes
static constexpr double	INF = 1.0e128
static const int	M_PW_PREC = 1
static constexpr double	SIGMA = 0.1
static constexpr double	RHO = 0.995
static constexpr double	OPTIM_GAP = 1.0e-6
static constexpr double	OPTIM_GAP_SAFEGUARD = 1.0e-5
static constexpr double	OPTIM_PFEAS = 1.0e-6
static constexpr double	OPTIM_DFEAS = 1.0e-6
static const int	OUTPUT_FREQ = 1
static const int	MAXITER = 200
static constexpr double	PCG_TOL_LIN = 1.0e-2
static constexpr double	PCG_TOL_QUAD = 1.0e-3
static constexpr double	MIN_PCGTOL = 1.0e-8
static constexpr double	RED_PCGTOL = 0.95
static const TYPE_COMP_DY	TYPE_COMP_DY_DEFAULT = CHOL_PWRS_PCG
static const TYPE_DIRECTION	TYPE_DIRECTION_DEFAULT = AUTOMATIC
static const TYPE_START_POINT	TYPE_START_POINT_DEFAULT = SIMPLE
static const TYPE_REG	TYPE_REG_DEFAULT = NO_REG
static constexpr double	FACTOR_REG_DEFAULT = 1.0e-6
static const bool	DEACTIVATELNK = true

Detailed Description

Main class for loading and solving problems.

Member Enumeration Documentation

enum BlockIP::OUTPUT

Enumerator:

NONE
SCREEN
FILE
BOTH

enum BlockIP::TYPE_COMP_DY

Enumerator:

CHOL_PWRS_PCG
FULL_CHOL
HYBRID_PCG
CHOL_THETA0_PCG

enum BlockIP::TYPE_DIRECTION

Enumerator:

NEWTON
SECOND_ORDER
AUTOMATIC

enum BlockIP::TYPE_PROBLEM

Enumerator:

LINEAR
QUADRATIC
NONLINEAR

enum BlockIP::TYPE_REG

Enumerator:

NO_REG
QUAD_REG
PROX_REG

enum BlockIP::TYPE_START_POINT

Enumerator:

SIMPLE
QUAD_EQCONS_PROB

Constructor & Destructor Documentation

BlockIP::BlockIP

(

)

Constructor.

BlockIP::BlockIP	(	TYPE_PROBLEM	type_objective,
		double *&	cost,
		double *&	qcost,
		void(*)(int n, double x[], double &fx, double Gx[], double Hx[], void *params)	fobj,
		void *	params,
		double *&	ub,
		double *&	rhs,
		int	numBlocks,
		bool	sameN,
		MatrixBlockIP	N[],
		bool	sameL,
		MatrixBlockIP	L[]
	)

Create an instance with problem in standard form: 0 <= variables <= ub and constraints = rhs.

Parameters

type_objective	Type of objective function, linear, quadratic or non-linear
cost	Linear cost of variables including slacks
qcost	Quadratic cost of variables including slacks
fobj	User function to calculate the objective function in a point. If it is is not NULL, cost and qcost must be NULL
params	User parameters to perform objective function calculations. Only can be used when fobj is not NULL
ub	Upper bounds, including slack bounds
rhs	Upper constraint limits, including linking constraints limits
numBlocks	Number of diagonal blocks
sameN	Define if the same matrix is used for each N block. If true array N must have dimension 1
N	Diagonal blocks
sameL	Define if the same matrix is used for each L block If true array L must have dimension 1
L	Linking constraints blocks

BlockIP::BlockIP	(	TYPE_PROBLEM	type_objective,
		double *&	cost,
		double *&	qcost,
		void(*)(int n, double x[], double &fx, double Gx[], double Hx[], void *params)	fobj,
		void *	params,
		double *&	lb,
		double *&	ub,
		double *&	lhs,
		double *&	rhs,
		int	numBlocks,
		bool	sameN,
		MatrixBlockIP	N[],
		bool	sameL,
		MatrixBlockIP	L[]
	)

Create an instance with general problem: lb <= variables <= ub and lhs <= constraints <= rhs.

Parameters

type_objective	Type of objective function, linear, quadratic or non-linear
cost	Linear cost of variables including slacks
qcost	Quadratic cost of variables including slacks
fobj	User function to calculate the objective function in a point. If it is is not NULL, cost and qcost must be NULL
params	User parameters to perform objective function calculations. Only can be used when fobj is not NULL
lb	Lower bounds, including slack bounds
ub	Upper bounds, including slack bounds
lhs	Lower constraint limits, including linking constraints limits
rhs	Upper constraint limits, including linking constraints limits
numBlocks	Number of diagonal blocks
sameN	Define if the same matrix is used for each N block. If true array N must have dimension 1
N	Diagonal blocks
sameL	Define if the same matrix is used for each L block If true array L must have dimension 1
L	Linking constraints blocks

Note

The problem must be converted to standard form in order to be solved, it will be converted automatically when minimize function is called. After that call the problem only may be written in standard form. Before that call the problem can be written in standard form with convert_to_std_and_write_mps function or in general form with write_mps function.

If sameN is used each constraint must be the same type in each block, e.g., if constraint i is an equality in block j must be an equality in all other blocks.

If sameN is used each free variable with zero value in Hessian must be the same type in each block, e.g., if variable i is free variable with zero value in Hessian of block j, it must be a free variable with zero value in Hessians of all other blocks.

BlockIP::~BlockIP

(

)

Destructor.

Member Function Documentation

void BlockIP::check_input_problem_and_compute_dimensions

(

)

Check the input problem and compute dimensions.

void BlockIP::check_sameN

(

)

Check if the problems satisfies the conditions when using sameN.

void BlockIP::convert_to_standard

(

)

Convert any problem to standard form, restrictions= rhs, variables >= 0, <= ub.

void BlockIP::convert_to_std_and_write_mps

(

const char *

filename

)

Convert the problem to standard form (if it is not already converted) and write a mps file.

Parameters

filename

File name without extension

void BlockIP::create_names

(

)

Create names for blocks, variables and constraints if does not exist.

void BlockIP::create_problem	(	TYPE_PROBLEM	type_objective,
		double *&	cost,
		double *&	qcost,
		void(*)(int n, double x[], double &fx, double Gx[], double Hx[], void *params)	fobj,
		void *	params,
		double *&	ub,
		double *&	rhs,
		int	numBlocks,
		bool	sameN,
		MatrixBlockIP	N[],
		bool	sameL,
		MatrixBlockIP	L[]
	)

Create a problem in standard form: 0 <= variables <= ub and constraints = rhs.

Parameters

type_objective	Type of objective function, linear, quadratic or non-linear
cost	Linear cost of variables including slacks
qcost	Quadratic cost of variables including slacks
fobj	User function to calculate the objective function in a point. If it is is not NULL, cost and qcost must be NULL
params	User parameters to perform objective function calculations. Only can be used when fobj is not NULL
ub	Upper bounds, including slack bounds
rhs	Upper constraint limits, including linking constraints limits
numBlocks	Number of diagonal blocks
sameN	Define if the same matrix is used for each N block. If true array N must have dimension 1
N	Diagonal blocks
sameL	Define if the same matrix is used for each L block If true array L must have dimension 1
L	Linking constraints blocks !

void BlockIP::create_problem	(	TYPE_PROBLEM	type_objective,
		double *&	cost,
		double *&	qcost,
		void(*)(int n, double x[], double &fx, double Gx[], double Hx[], void *params)	fobj,
		void *	params,
		double *&	lb,
		double *&	ub,
		double *&	lhs,
		double *&	rhs,
		int	numBlocks,
		bool	sameN,
		MatrixBlockIP	N[],
		bool	sameL,
		MatrixBlockIP	L[]
	)

Create a general problem: lb <= variables <= ub and lhs <= constraints <= rhs.

Parameters

type_objective	Type of objective function, linear, quadratic or non-linear
cost	Linear cost of variables including slacks
qcost	Quadratic cost of variables including slacks
fobj	User function to calculate the objective function in a point. If it is is not NULL, cost and qcost must be NULL
params	User parameters to perform objective function calculations. Only can be used when fobj is not NULL
lb	Lower bounds, including slack bounds
ub	Upper bounds, including slack bounds
lhs	Lower constraint limits, including linking constraints limits
rhs	Upper constraint limits, including linking constraints limits
numBlocks	Number of diagonal blocks
sameN	Define if the same matrix is used for each N block. If true array N must have dimension 1
N	Diagonal blocks
sameL	Define if the same matrix is used for each L block If true array L must have dimension 1
L	Linking constraints blocks

Note

The problem must be converted to standard form in order to be solved, it will be converted automatically when minimize function is called. After that call the problem only may be written in standard form. Before that call the problem can be written in standard form with convert_to_std_and_write_mps function or in general form with write_mps function.

If sameN is used each constraint must be the same type in each block, e.g., if constraint i is an equality in block j must be an equality in all other blocks.

If sameN is used each free variable with zero value in Hessian must be the same type in each block, e.g., if variable i is free variable with zero value in Hessian of block j, it must be a free variable with zero value in Hessians of all other blocks. !

void BlockIP::eval_fobj ( double  x[], double &  fx, double  Gx[], double  Hx[]  )

inlineprivate

void BlockIP::free_memory

(

)

Free all the memory.

int BlockIP::get_cgit ( )

inline

Return number of overall PCG iterations.

double BlockIP::get_constant_fobj ( )

inline

Get the constant to add to the objective function.

bool BlockIP::get_deactivateLnk ( )

inline

Get flag on deactivation of linking.

double BlockIP::get_dobj ( )

inline

Return dual objective function.

double BlockIP::get_factor_reg ( )

inline

Get factor of regularization.

double BlockIP::get_fobj ( )

inline

Return objective function.

string * BlockIP::get_full_cons_names

(

)

string * BlockIP::get_full_var_names

(

)

double BlockIP::get_inf ( )

inline

Get infinity value to be considered.

double BlockIP::get_init_pcgtol ( )

inline

Get initial tolerance for the conjugate gradient (by default PCG_TOL_LIN if linear problem, PCG_TOL_QUAD if quadratic or -convex- nonlinear)

int BlockIP::get_it ( )

inline

Return number of IP iterations of minimization.

int BlockIP::get_k_blocks ( )

inline

Get the number of blocks.

Returns

field k_blocks

int BlockIP::get_km ( )

inline

Get the number of constraints without linking constraints.

Returns

field km

int BlockIP::get_kn ( )

inline

Get the number of variables without linking constraints.

Returns

field kn

int BlockIP::get_l_link ( )

inline

Get the number of linking constraints.

Returns

field l_link

int BlockIP::get_m_cons ( )

inline

Get the number of constraints of the problem to be optimized.

Returns

field m_cons

int BlockIP::get_m_pw_prec ( )

inline

Get number of terms used as preconditioner of the power series expansion of (D-C'*B^{-1}*B)^{-1}.

double BlockIP::get_maxit_pcg ( )

inline

Get maximum number of pcg iterations (if 0, then the value will be computed by the code)

int BlockIP::get_maxiter ( )

inline

Get maximum number of IP iterations.

double BlockIP::get_min_pcgtol ( )

inline

Get minimum tolerance for the conjugate gradient.

int BlockIP::get_n_vars ( )

inline

Get the number of variables of the problem to be optimized.

Returns

field n_vars

void BlockIP::get_names	(	const string *&	blockNames,
		const string *&	varNames,
		const string *&	consNames
	)

Get the names of the problem.

Parameters

blockNames	Block names of N blocks, dimension k_blocks
varNames	Variable names including slacks, dimension n_vars
consNames	Constraint names including linking constraints, dimension m_cons

int BlockIP::get_num_cons ( )

inline

Get the number of constraints of the original problem.

Returns

Number of constraints of the original problem

int BlockIP::get_num_vars ( )

inline

Get the number of variables of the original problem.

Returns

Number of variables of the original problem

double BlockIP::get_optim_dfeas ( )

inline

Get dual feasibility tolerance.

double BlockIP::get_optim_gap ( )

inline

Get optimality gap tolerance.

double BlockIP::get_optim_pfeas ( )

inline

Get primal feasibility tolerance.

BlockIP::OUTPUT BlockIP::get_output ( )

inline

Get type of output.

int BlockIP::get_output_freq ( )

inline

Get output information lines will be printed each output_freq IP iterations.

double BlockIP::get_red_pcgtol ( )

inline

Get reduction of pcg_tol at each IP iteration.

double BlockIP::get_rho ( )

inline

Get reduction of the step-length for the primal and dual variables at each IP iteration.

bool BlockIP::get_show_specrad ( )

inline

Get show_specrad.

double BlockIP::get_sigma ( )

inline

Get reduction of the centrality parameter at each IP iteration.

StdForm * BlockIP::get_std_form ( bool  keepStdFormAfterDelete = false )

inline

Get the StdForm related to the problem.

Parameters

keepStdFormAfterDelete

Set if the StdForm related to the problem will be deleted or not when BlockIP will delete

Returns

The StdForm related to the problem, if the problem was created in standard form, return NULL.

BlockIP::TYPE_COMP_DY BlockIP::get_type_comp_dy ( )

inline

Get how dy direction is computed.

BlockIP::TYPE_DIRECTION BlockIP::get_type_direction ( )

inline

Get which direction is computed.

BlockIP::TYPE_REG BlockIP::get_type_reg ( )

inline

Get type of regularization.

BlockIP::TYPE_START_POINT BlockIP::get_type_start_point ( )

inline

Get how starting point is computed.

const double * BlockIP::get_w ( )

inline

Get dual variable of problem optimized (of x<=u) i=1..n_vars.

double BlockIP::get_w ( int  i )

inline

Get value of dual variable w(i) (of x<=u) of problem optimized i=1..n_vars.

WHO_PERMUTES BlockIP::get_whoperm ( )

inline

Get whopermutes constraints-related information: the Cholesky factorization of the user of it.

const double * BlockIP::get_x ( )

inline

Get primal variables of problem optimized.

double BlockIP::get_x ( int  i )

inline

Get value of primal variable x(i) of problem optimized i=1..n_vars.

const double * BlockIP::get_y ( )

inline

Get dual variables of problem optimized (of equality constraints)

double BlockIP::get_y ( int  i )

inline

Get value of dual variable y(i) of problem optimized (of equality constraints) i=1..m_cons.

const double * BlockIP::get_z ( )

inline

Get dual variable of problem optimized (of x>=0) i=1..n_vars.

double BlockIP::get_z ( int  i )

inline

Get value of dual variable z(i) (of x>=0) of problem optimized i=1..n_vars.

void BlockIP::initialize

(

)

Initializes class attributes.

void BlockIP::min_Aty ( double *  Aty, const double *  y  )

private

Computes and returns A'*y A is the constraints matrix, formed by k diagonal blocks N and last rows with k L matrices plus the Identity for slacks. Computes Ni'*yi + Li'*y_{k+1} for all block i=1..k; and I*y_{k+1} for last block k+1.

Parameters

y[1:km+l],:	dual variables. Values of duals y_{k+1} of inactive linking constraints must be 0.
Aty[kn+l],:	result of A'*y. Only values of duals y_{k+1} of active linking constraints are added (equivalently, the components of y_{k+1} of inactives are 0).

void BlockIP::min_Ax ( double *  Ax, const double *  x  )

private

Computes and returns A*x A is the constraints matrix, formed by k diagonal blocks N and last rows with k L matrices plus the Identity for slacks. Computes Ni*xi for all block i=1..k; and sum{1..k}Li*xi + x_{k+1}, for the active constraints

Parameters

x[1:kn+l],:	primal variables
Ax[km+l],:	result of A*x. Only values of active linking onstraints are filled

void BlockIP::min_backup_inactive ( int  lnk )

private

void BlockIP::min_check_Newton_direction_is_correct ( )

private

For debugging purposes, checks direction satisfies Newton system

| -H A' I -I | |dx| |rc | | A | |dy| = |rb | | Z X | |dz| |rxz| | -W S | |dw| |rsw|

If QUAD_REG then it considers H:= H+qreg for block variables

void BlockIP::min_check_predictor_corrector_direction_is_correct ( )

private

For debugging purposes, checks predictor-corrector direction satisfies system

| -H A' I -I | |dx| |rc | | A | |dy| = |rb | | Z X | |dz| |sigma*mu-XZe-Dx_p Dz_p e| | -W S | |dw| |sigma*mu-SWe-Ds_p Dw_p e|

If QUAD_REG then it considers H:= H+qreg for block variables

void BlockIP::min_check_predictor_direction_is_correct ( )

private

For debugging purposes, checks predictor direction satisfies system

| -H A' I -I | |dx| |rc | | A | |dy| = |rb | | Z X | |dz| |-XZe| | -W S | |dw| |-SWe|

If QUAD_REG then it considers H:= H+qreg for block variables

void BlockIP::min_compute_direction ( )

private

Decides whether to use standard Newton direction or second order direction (Mehrotra direction).

void BlockIP::min_compute_dy_cholpcg ( double *  dy, double *  rhsdy, bool  only_solve = false  )

private

Solves (A*Theta*A')dy = rhsdy by Cholesky for blocks and PCG for linking constraints. The structure of (A*Theta*A') dy = rhsdy is

| B C | |dy1| |rhsdy1| | C' D | |dy2| = |rhsdy2|

so it can be solved as

(D-C'*B^{-1}*C) dy2 = rhsdy2-C'*B^{-1}*rhsdy1 B dy1= rhsdy1-C*dy2

where B is made of k_blocks N_i*Theta_i*N_i'.

Note that system (A*Theta*A') is really of dimension km+numActiveLnk, then we need to implicitly consider only the active linking constraints when applying the PCG.

Parameters

dy[1:km+l],:	on input is the last dy solution computed, to be used as starting solution fof PCG; on output is the solution of (A*Theta*A)dy= rhsdy
rhsdy[1:km+l],:	right-hand-side of PCG system.
only_solve,:	if true, only numeric solve (neither symbolic nor numeric factorization are needed) since a previous call with the same A*Theta*A' but different rhsdy was made (for instance, in corrector step of predictor-corrector heuristic)

void BlockIP::min_compute_dy_fullchol ( double *  dy, double *  rhsdy, bool  only_solve = false  )

private

Solves (A*Theta*A')dy = rhsdy by Cholesky of full system. Note that system (A*Theta*A') is really of dimension km+numActiveLnk, but we reuse the initial symbolic factorization, where the dimension of the system was km+l_link. Instead of using a row of the identity matrix for the rows of inactive constraints to avoid a nonsingular system, we consider a 0, since Cholesky deals with 0 pivots (semidefinie matrices). This way, the values of dy are correctly computed.

Parameters

dy[1:km+l],:	on output is the solution of (A*Theta*A)dy= rhsdy
rhsdy[1:km+l],:	on input right-hand-side of (A*Theta*A)dy= rhsdy
only_solve,:	if true, only numeric solve (neither symbolic nor numeric factorization are needed) since a previous call with the same A*Theta*A' but different rhsdy was made (for instance, in corrector step of predictor-corrector heuristic)

void BlockIP::min_compute_mu ( )

private

Computes and returns the centrality at current point, computed as mu= (x'z+s'w)/( (kn+number_of_active_linking) + (number of variables with upper bounds)). For efficiency of the implementation, w and s of variables without upper bound are 0 and a finite (large number)

void BlockIP::min_compute_s ( )

private

Computes the slacks of upper bounds of block variables and slacks of linking constraints: s= u-x. Only slacks of upper bounds of slacks of active linking constraints filled (for inactive linking, s is 0). If some u=inf, then s=u-x~=inf which is OK for deactivating the complementarity constraints s*w=sigma*mu. There is no need in setting s=inf for vars in listWithoutUb since s=u-x will be very close (if not equal) to inf. This also applies to (non-splitted) free vars, whose upper bound is also inf.

void BlockIP::min_compute_sigma_psi ( )

private

Computes sigma (reduction of the centrality parameter at each IP iteration) and psi (reduction of dx*dz vector in the rhs of corrector system of predictor-corrector heuristic)

void BlockIP::min_compute_Theta ( )

private

void BlockIP::min_compute_Theta0 ( )

private

void BlockIP::min_Ctv ( double *  v, double *  Ctv  )

private

Given v[1:km], this routine computes C'*v= sum{i in 1..k} L_i*Theta_i*N_i'*v_i, which is a vector of dimension l_link. The result vector Ctv only contains values of active linking, components of inactive constraints are 0.

Parameters

v[1:km],:	input vector
Ctv[1:l_link],:	output vector containing C'*v

void BlockIP::min_Cv ( double *  v, double *  Cv  )

private

Given v[l], this routine computes C*v= [N_i*Theta_i*L_i'*v i=1..k], C*v is a vector [1:km]. Product L_i'*v for all components is made assuming that inactive components of v or inactive rows of L are 0 (otherwise, there is a bug in the code and the function will provide wrong results...)

Parameters

v[1:l_link],:	input vector
Cv[1:km],:	output vector containing C*v

void BlockIP::min_debug_variables_of_inactivelnk ( )

private

For debugging purposes, check variables and data of inactive slacks and constraints are 0

void BlockIP::min_debug_variables_without_ub ( )

private

For debugging purposes, check variables and data of variables without upper bounds are 0

void BlockIP::min_debug_write_variables ( )

private

For debugging purposes, write variables to file

void BlockIP::min_free_memory ( )

private

Deletes arrays locally needed for the interior-point algorithm

void BlockIP::min_initializations ( )

private

Initialize parameters, tolerances, etc for BlockIP minimization algorithm.

void BlockIP::min_KKT_residuals ( )

private

Compute the residuals of mu-KKT equations Ax=b (primal feasibility) A'*y+z-w-Gx=0 (dual feasibility) XZe= sigma*mu*e (complementarity x*z) SWe= sigma*mu*e (complementarity s*w)

The residuals are:

• rb[km+l]= b-Ax : primal infeasibility. Positions of inactive linking constraints are set to 0.
• rc[kn+l]= Gx-A'y-z+w : dual infeasibility. Positions of slacks of inactive linking constraints are set to 0. Gx is the gradient of the objective function. If it is nonlinear, is stored in Gx. If it is linear Gx=c. If it is quadratic, Gx= Qx+c. If QUAD_REG regularization then Qreg*x term has bo be added to rc (rc=Gx+qreg*x-A'y-z+w) only for block variables.
• rxz[kn+l]= sigma*mu*e- X*z : complementarity residual for Xz= mu*e, decreasing mu (see below) by sigma. Positions of inactive linking constraints are set to 0.
• rsw[kn+l]= sigma*mu*e- S*w : complemetarity residual for Sw= mu*e, decreasing mu (see below) by sigma. Positions of inactive linking constraints are set to 0.
• normrc: norm of rc weighted by (1+|c|)
• normrb: norm of rb weighted by (1+|b|)
• mu= (x'z+s'w)/(kn+length(activelnk)) : centrality parameter

Residuals rxz and rsw only needed if NEWTON direction is used. SECOND ORDER direction do not need rxz and rsw (predictor-corrector will compute later their own right-hand-sides for the complementarity equations).

void BlockIP::min_Newton_direction ( )

private

Computes the standard Newton direction by solving normal equations A*Theta*A' according to type_comp_dy:

The procedure:

For the primal convex nonlinear problem

(P) min f(x) subject to Ax=b, -x<=0, x-u<=0

the Lagrangian is

L(x,y,z,w)= f(x)+y'*(b-Ax)+z'*(-x)+w'*(x-u)

and the gradient respect to x is (Gx is gradient of f(x); Gx=c if f(x) is linear, Gx= Qx+c if f(x) is quadratic):

L(x,y,z,w)= Gx - A'*y - z + w

and the dual problem is

(D) min L(x,y,z,w) s. to L(x,y,z,w)=0 w>=0, z>=0

The condition L(x,y,z,w)=0 is the dual feasibility below.

The mu-KKT conditions are (Gx is gradienf of f(x)):

A'*y+z-w-Gx=0 (dual feasibility) Ax=b (primal feasibility) XZe= sigma*mu*e (complementarity x*z) SWe= sigma*mu*e (complementarity s*w)

The Newton system is (S= U-X, and Hx is Hessian of f(x))

| -Hx A' I -I | |dx| |rc | | A | |dy| = |rb | | Z X | |dz| |rxz| | -W S | |dw| |rsw|

which is solved by:

Theta= (Hx + S^{-1}*W + X^{-1}*Z)^{-1}

r= rc + S^{-1}*rsw - X^{-1}*rxz (A*Theta*A')dy= rb+A*Theta*r dx= Theta*(A'dy-r) dw= S^{-1}(rsw + W*dx) dz= rc + dw + Hx*dx - A'*dy

If QUAD_REG regularization is applied then Gx is (Gx+Qreg*x) and Hx is (Hx+Qreg). rc and Theta were previosuly computed considering (Gx+Qreg*x) and (Hx+Qreg). Here we have to use (Hx+Qreg) too.

void BlockIP::min_normb_normc ( )

private

Compute the L2 norm of b[] and c[]. For the linking constraints, only the active are considered.

void BlockIP::min_objective_functions_linquad ( )

private

void BlockIP::min_objective_functions_nonlin ( )

private

bool BlockIP::min_optimal ( )

private

int BlockIP::min_PCG_Hv	(	int	nn,
		double *	v,
		double *	Hv
	)

Computes Hv=(D-C'*B^(-1)*C)v Returns vector Hv of dimension l_link, assuming components of inactives are 0 (if there is no error in the code). Then we avoid using (packing and unpacking) a vector of dimension numActiveLnk.

Parameters

v[1:l_link],:	input vector
Hv[1:l_link],:	output vector containing H*v

int BlockIP::min_PCG_Hz_eq_r	(	int	nn,
		double *	zz,
		double *	rr
	)

Approximate solution of H*zz = rr, using m_pw_prec terms of the power series expansion H^{-1}= [sum{i=0..infinity} (D^{-1}*(C'B^{-1}B)C)^i]*D^{-1} Recommended values are m_pw_prec= 1 or 2, to avoid expensive computations We use vectors of dimension l_link, assuming components of inactives are 0 (if there is no error in the code). Then we avoid using (packing and unpacking) a vector of dimension numActiveLnk.

Parameters

rr[1:l_link],:	rhs of system of equations
zz[1:l_link],:	solution of the system

int BlockIP::min_PCG_Theta0z_eq_r	(	int	nn,
		double *	zz,
		double *	rr
	)

Solves Theta0*zz = rr, where Theta0 are the components of Theta associated to active linking constraints. This preconditioner may be useful when the space of linking constraints is "close" (principal angles not large) to the space of block constraints. We use vectors of dimension l_link, assuming components of inactives are 0 (if there is no error in the code). Then we avoid using (packing and unpacking) a vector of dimension numActiveLnk.

Parameters

rr[1:l_link],:	rhs of system of equations
zz[1:l_link],:	solution of the system

void BlockIP::min_preprocess ( )

private

Simple check for

• Order matrices stored in general packed format
• Convert general matrices stored in column-wise packed format to row-wise packed format
• no-zero upper bounds (variables and slacks) not allowed. Set a very small upper bound to guarantee an interior solution. Zero variables are not removed to preserve the topology of the constraints for each block then we can maintain one single symbolic factorization.
• q positive semidefinite; is set to 0 if problem is linear
• initialization of activelnk if empty
• initialization of free vars arrays if empty and all variables are unfree
• check infinity upper bounds
• check no rhs term is infinity (otherwise the algorithm will fail)

void BlockIP::min_restore_inactive ( int  lnk )

private

void BlockIP::min_second_order_predictor_corrector_direction ( )

private

Computes the second order predictor-corrector Mehrotra direction by solving normal equations A*Theta*A' according to type_comp_dy:

The procedure:

For the primal convex nonlinear problem

(P) min f(x) subject to Ax=b, -x<=0, x-u<=0

the Lagrangian is

L(x,y,z,w)= f(x)+y'*(b-Ax)+z'*(-x)+w'*(x-u)

and the gradient respect to x is (Gx is gradient of f(x); Gx=c if f(x) is linear, Gx= Qx+c if f(x) is quadratic):

L(x,y,z,w)= Gx - A'*y - z + w

and the dual problem is

(D) min L(x,y,z,w) s. to L(x,y,z,w)=0 w>=0, z>=0

The condition L(x,y,z,w)=0 is the dual feasibility below.

The mu-KKT conditions are (Gx is gradienf of f(x)):

A'*y+z-w-Gx=0 (dual feasibility) Ax=b (primal feasibility) XZe= sigma*mu*e (complementarity x*z) SWe= sigma*mu*e (complementarity s*w)

The predictor-corrector direction is computed as (S= U-X, and Hx is Hessian of f(x)):

1. Predictor direction

| -Hx A' I -I | |dx_p| |rc | | A | |dy_p| = |rb | | Z X | |dz_p| |-XZe| | -W S | |dw_p| |-SWe|

1. Compute sigma and psi

1. Predictor+Corrector direction

| -Hx A' I -I | |dx| |rc | | A | |dy| = |rb | | Z X | |dz| |sigma*mu-XZe-Dx_p Dz_p e| | -W S | |dw| |sigma*mu-SWe-Ds_p Dw_p e|

Each of the above Newton-like systems are solved as:

| -Hx A' I -I | |dx| |rc | | A | |dy| = |rb | | Z X | |dz| |rhs_xz| | -W S | |dw| |rhs_sw|

Theta= (Hx + S^{-1}*W + X^{-1}*Z)^{-1}

r= rc + S^{-1}*rhs_sw - X^{-1}*rhs_xz (A*Theta*A')dy= rb+A*Theta*r dx= Theta*(A'dy-r) dw= S^{-1}(rhs_sw + W*dx) dz= rc + dw + Hx*dx - A'*dy

If QUAD_REG regularization is applied then Gx is (Gx+Qreg*x) and Hx is (Hx+Qreg). rc and Theta were previosuly computed considering (Gx+Qreg*x) and (Hx+Qreg). Here we have to use (Hx+Qreg) too.

void BlockIP::min_solve_NThetaNt ( double *  v )

private

Solves k_blocks systems N_i*Theta_i*N_i' dy_i= rhsdy_i i=1..k_blocks where N_i can be different matrices or the same one (if sameN).

Parameters

v[1:km],:	on input contains the righ-hand-side for the k_blocks systems
v[1:km],:	on output it contain the solutions to the k_blocks systems

void BlockIP::min_starting_point ( )

private

Computation of starting point:

1. simple initialization
2. Mehrotra-like initilization (solves quadratic equality constrained problem)

void BlockIP::min_steplengths ( double &  talpha_p, double &  talpha_d, double  tdx[], double  tdz[], double  tdw[]  )

private

void BlockIP::min_test_newpoint ( )

private

void BlockIP::min_update_inactivelnk ( )

private

Detect if new linking constraints must be considered inactive This is only done if we are close enough to the optimal solution; we use (gap < 1.0), where gap was computed before as abs(pobj-dobj)/(1+abs(pobj)). This only works for linear/quadr. functions, but for nonlinear ones this update is not used (see below an explanation). Constraint i inactived if its slack x(i) is out of bounds and lagrange multiplier is close to 0: (x(i)>>0 and x(i)<<u(i) and y(i)~0) and ( c(i)==0 and q(i)==0). The term (c(i)==0 and q(i)==0) guarantees that slacks appearing in the objective function are not removed. This could even work for nonlinear functions: if component of gradient and Hessian of slack is 0, the slack could be removed (**** be careful: this can provide strange results if for some reason the first and second derivatives are 0 only in that iteration; for safety, the inactivation is not performed for nonlinear objective functions ***) We also avoid inactivation of constraints with small bounds on slacks (i.e, u<0.1), since these constraints are always quasi-active

void BlockIP::min_update_newpoint ( )

private

void BlockIP::min_update_qreg ( )

private

Updates regularization term

void BlockIP::min_write_current_point_info ( )

private

void BlockIP::min_write_problem_information ( )

private

int BlockIP::minimize

(

)

void BlockIP::read_BlockIP_format

(

const char *

filename

)

Create a problem from a BlockIP format file.

Parameters

filename

File name with extension

Note

The problem must be converted to standard form in order to be solved, it will be converted automatically when minimize function is called. After that call the problem only may be written in standard form. Before that call the problem can be written in standard form with convert_to_std_and_write_mps function or in general form with write_mps function. !

void BlockIP::read_mps

(

const char *

filename

)

Create a problem from a mps file.

Parameters

filename

File name with extension

Note

Token : separate the block name and the variable or constraint name When a constraint or variable does not have a block name it is considered a linking constraint or a slack

If some variable have as a name "Constant" with no block, this variable is considered a constant to add in the objective function

Slacks are optional, but if one is defined, all slacks must be defined, to relate an slack to one linking constraint, this slack must have a 1 coeficient in the related linking constraint.

In QUADOBJ only cannot be relation between two different variables.

The problem must be converted to standard form in order to be solved, it will be converted automatically when minimize function is called. After that call the problem only may be written in standard form. Before that call the problem can be written in standard form with convert_to_std_and_write_mps function or in general form with write_mps function. !

void BlockIP::set_constant_fobj ( double  constant )

inline

Set the constant to add to the objective function.

void BlockIP::set_deactivateLnk ( bool  deactivateLnk = DEACTIVATELNK )

inline

Set flag on deactivation of linking.

void BlockIP::set_defaults_minimize

(

)

Set to default values all the parameters that control the interior-point algorithm

void BlockIP::set_factor_reg ( double  factor_reg = FACTOR_REG_DEFAULT )

inline

Set factor of regularization.

void BlockIP::set_inf ( double  inf = INF )

inline

Set threshold value to be used as infinity (x>inf -> x is inf)

void BlockIP::set_init_pcgtol ( double  init_pcgtol )

inline

Set initial tolerance for the conjugate gradient (by default PCG_TOL_LIN if linear problem, PCG_TOL_QUAD if quadratic or -convex- nonlinear)

void BlockIP::set_m_pw_prec ( int  m_pw_prec = M_PW_PREC )

inline

Set number of terms used as preconditioner of the power series expansion of (D-C'*B^{-1}*B)^{-1}.

void BlockIP::set_maxit_pcg ( double  maxit_pcg = 0 )

inline

Set maximum number of pcg iterations (if 0, then the value will be computed by the code)

void BlockIP::set_maxiter ( int  maxiter = MAXITER )

inline

Set maximum number of IP iterations.

void BlockIP::set_min_pcgtol ( double  min_pcgtol = MIN_PCGTOL )

inline

Set minimum tolerance for the conjugate gradient.

void BlockIP::set_names	(	string *	blockNames,
		string *	varNames,
		string *	consNames,
		bool	copy_vectors = true
	)

Set the names of the problem TODO change copy_vectors.

Parameters

blockNames	Block names of N blocks
varNames	Variable names including slacks
consNames	Constraint names including linking constraints
copy_vectors	If is true the user must free the memory of arguments (arrays), If false the user must allocate the memory with new, after this call the user will not have control of the arguments (arrays) anymore and MatrixBlockIP will delete the arguments

Note

If some variable is NULL BlockIP keeps its own name

void BlockIP::set_optim_dfeas ( double  optim_dfeas = OPTIM_DFEAS )

inline

Set dual feasibility tolerance.

void BlockIP::set_optim_gap ( double  optim_gap = OPTIM_GAP )

inline

Set optimality gap tolerance.

void BlockIP::set_optim_pfeas ( double  optim_pfeas = OPTIM_PFEAS )

inline

Set primal feasibility tolerance.

void BlockIP::set_output ( OUTPUT  output = SCREEN )

inline

Set type of output.

void BlockIP::set_output_freq ( int  output_freq = OUTPUT_FREQ )

inline

Set output information lines will be printed each output_freq IP iterations.

void BlockIP::set_red_pcgtol ( double  red_pcgtol = RED_PCGTOL )

inline

Set reduction of pcg_tol at each IP iteration.

void BlockIP::set_rho ( double  rho = RHO )

inline

Set reduction of the step-length for the primal and dual variables at each IP iteration.

void BlockIP::set_show_specrad ( bool  show_specrad = false )

inline

Set show_specrad at each IP iteration.

void BlockIP::set_sigma ( double  sigma = SIGMA )

inline

Set reduction of the centrality parameter at each IP iteration.

void BlockIP::set_type_comp_dy ( TYPE_COMP_DY  type_comp_dy = TYPE_COMP_DY_DEFAULT )

inline

Set how dy direction is computed.

void BlockIP::set_type_direction ( TYPE_DIRECTION  type_direction = TYPE_DIRECTION_DEFAULT )

inline

Set which direction is computed.

void BlockIP::set_type_reg ( TYPE_REG  type_reg = TYPE_REG_DEFAULT )

inline

Set type of regularization.

void BlockIP::set_type_start_point ( TYPE_START_POINT  type_start_point = TYPE_START_POINT_DEFAULT )

inline

Set how starting point is computed.

void BlockIP::set_whoperm ( WHO_PERMUTES  whoperm_ = CHOLESKY )

inline

Set whopermutes constraints-related information: the Cholesky factorization of the user of it.

void BlockIP::write_BlockIP_format

(

const char *

filename

)

Write the problem in BlockIP format file.

Parameters

filename

File name without extension

void BlockIP::write_mps

(

const char *

filename

)

Write a mps file.

void BlockIP::write_mps	(	const char *	filename,
		TYPE_PROBLEM	type_objective,
		double	cost[],
		double	qcost[],
		double	lb[],
		double	ub[],
		double	lhs[],
		double	rhs[],
		int	numBlocks,
		bool	sameN,
		MatrixBlockIP	N[],
		bool	sameL,
		MatrixBlockIP	L[],
		string *	blockNames = NULL,
		string *	varNames = NULL,
		string *	consNames = NULL,
		double	constantFObj = 0
	)

Write a mps file.

Parameters

filename	File name without extension
type_objective	Type of objective function, linear or quadratic
cost	Linear cost of variables including slacks
qcost	Quadratic cost of variables including slacks
lb	Lower bounds, including slack bounds
ub	Upper bounds, including slack bounds
lhs	Lower constraint limits, including linking constraints limits
rhs	Upper constraint limits, including linking constraints limits
numBlocks	Number of diagonal blocks
sameN	Define if the same matrix is used for each N block. If true array N must have dimension 1
N	Diagonal blocks
sameL	Define if the same matrix is used for each L block If true array L must have dimension 1
L	Linking constraints blocks
blockNames	Block names of N blocks
varNames	Variable names including slacks
consNames	Constraint names including linking constraints

void BlockIP::write_problem

(

const char *

filename

)

Write all data related to the problem into a file.

Parameters

filename

File name where output the data

Member Data Documentation

MatrixBlockIP* BlockIP::A

Full matrix (likely made from N and L)

double BlockIP::alpha_d

private

step lengths

double BlockIP::alpha_p

private

double * BlockIP::Aty

private

double* BlockIP::Ax

private

BackupLnk* BlockIP::backupLnk

private

string* BlockIP::blockNames

Block names of N blocks.

int BlockIP::cgit

private

number of CG iterations of this IPM iteration

bool BlockIP::comp_specrad

private

Whether the estimation of spectral radius has to be shown in the output.

string* BlockIP::consNames

Constraint names including linking constraints.

double BlockIP::constantFObj

Constant to add in the objective function.

double* BlockIP::cost

Linear cost of variables including slacks.

double * BlockIP::Cvaux

private

double * BlockIP::Cvaux2

private

Auxiliary vectors for interior-point algorithm.

MatrixBlockIP* BlockIP::D

D = Theta_(k+1) + sum{i..k} L_i*Theta_i*L_i'.

const bool BlockIP::DEACTIVATELNK = true

staticprivate

bool BlockIP::deactivateLnk

private

Vectors to store information from PCG to later compute Ritz values.

Setting deactivateLnk= false the user may avoid the deactivation of linking

bool BlockIP::deleteMatrices

private

To delete matrices when created by BlockIP (read_mps)

double BlockIP::dobj

dual objective

double * BlockIP::dsdw

private

for rhs of corrector system, from dx,dw,ds and dw of predictor direction

double * BlockIP::dw

private

Primal and dual optimization directions.

double* BlockIP::dx

private

double* BlockIP::dxdz

private

double * BlockIP::dy

private

double * BlockIP::dz

private

double BlockIP::epsmach

private

Stores computed machine epsilon.

double BlockIP::est_specrad

private

Estimation of spectral radius of D^{-1}C'*B^{-1}*B (computed through Ritz values)

double BlockIP::factor_reg

private

Type of regularization performed.

constexpr double BlockIP::FACTOR_REG_DEFAULT = 1.0e-6

staticprivate

void(* BlockIP::fobj)(int n, double x[], double &fx, double Gx[], double Hx[], void *params)

User function to calculate the objective function in a point.

double BlockIP::fx

double BlockIP::gap

private

Duality gap.

double * BlockIP::Gx

double * BlockIP::Hx

Objective function value, Gradient and Hessian in the actual point.

constexpr double BlockIP::INF = 1.0e128

staticprivate

double BlockIP::inf

private

Infinity value.

For the default values see the constant terms in BlockIP.h

int* BlockIP::ini_m

Begin to blocks 1..k and linking for constraints.

int* BlockIP::ini_n

Begin to blocks 1..k and linking for variables.

double BlockIP::init_pcgtol

private

Initial tolerance for the conjugate gradient (by default PCG_TOL_LIN if linear problem, PCG_TOL_QUAD if quadratic or -convex- nonlinear)

bool BlockIP::initialized

private

To control if the attributes have been initialized.

int BlockIP::input_problem

private

bool BlockIP::inStdForm

To control if the problem is in standard form.

bool* BlockIP::isActiveLnk

private

linking i is active if isActiveLnk[i] is true

bool* BlockIP::isFreeVar

private

Variable i is marked as free if isFreeVar[i] is true.

int BlockIP::it

private

IPM iterations.

int BlockIP::k_blocks

Number of diagonal blocks.

bool BlockIP::keepStdFormAfterDelete

private

To keep stdForm when BlockIP will delete.

int BlockIP::km

Number of constraints without linking constraints.

int BlockIP::kn

Number of variables without slacks of linking.

MatrixBlockIP* BlockIP::L

Linking constraints blocks.

int BlockIP::l_link

Number of linking constraints.

double* BlockIP::lb

Lower bounds, including slack bounds.

double* BlockIP::lhs

Lower constraint limits, including linking constraints limits.

int* BlockIP::listActiveLnk

private

list of active linking, j=listActiveLnk[i], i=1..numActiveLnk, j in {1,..,l_link})

int* BlockIP::listFreeVars

private

List of variables marked as free, j=listFreeVars[i], i=1..numFreeVars, j in {1,..,n_vars})

int* BlockIP::listInactiveLnk

private

list of inactive linking, j=listInactiveLnk[i], i=1..numInactiveLnk, j in {1,..,l_link})

int* BlockIP::listUnfreeVars

private

List of variables not marked as free, j=listUnfreeVars[i], i=1..numUnfreeVars, j in {1,..,n_vars})

int* BlockIP::listWithoutUb

private

int* BlockIP::listWithUb

private

int BlockIP::m_cons

Number of constraints including linking constraints.

const int BlockIP::M_PW_PREC = 1

staticprivate

int BlockIP::m_pw_prec

private

Number of terms of the power series expansion of (D-C'*B^{-1}*B)^{-1} used as preconditioner.

double BlockIP::maxit_pcg

private

Maximum number of pcg iterations (if 0, then the value will be computed by the code)

const int BlockIP::MAXITER = 200

staticprivate

int BlockIP::maxiter

private

Maximum number of IP iterations.

constexpr double BlockIP::MIN_PCGTOL = 1.0e-8

staticprivate

double BlockIP::min_pcgtol

private

Minimum tolerance for the conjugate gradient.

double BlockIP::mu

private

double BlockIP::mu0

private

Centrality parameter; mu0 is value of mu at starting point.

MatrixBlockIP* BlockIP::N

Diagonal blocks.

int BlockIP::n_vars

Number of variables including slacks.

double BlockIP::normb

private

double BlockIP::normc

private

double BlockIP::normrb

private

double BlockIP::normrc

private

norms of rhs, costs, rb, and rc

int BlockIP::numActiveLnk

private

Number of active lnk.

int BlockIP::numBackupLnk

private

int BlockIP::numFreeVars

private

Number of variables marked as free.

int BlockIP::numInactiveLnk

private

Number of inactive lnk.

int BlockIP::numInactiveLnkWithUb

private

Number of inactive lnk with upper bound.

int BlockIP::numUnfreeVars

private

Number of variables not marked as free.

int BlockIP::numWithoutUb

private

int BlockIP::numWithUb

private

Number of variables with upper bound (not infinity)

constexpr double BlockIP::OPTIM_DFEAS = 1.0e-6

staticprivate

double BlockIP::optim_dfeas

private

Dual feasibility tolerance.

constexpr double BlockIP::OPTIM_GAP = 1.0e-6

staticprivate

double BlockIP::optim_gap

private

Optimality gap tolerance.

constexpr double BlockIP::OPTIM_GAP_SAFEGUARD = 1.0e-5

staticprivate

constexpr double BlockIP::OPTIM_PFEAS = 1.0e-6

staticprivate

double BlockIP::optim_pfeas

private

Primal feasibility tolerance.

int* BlockIP::origNCons

Number of original constraints for each block.

int* BlockIP::origNVars

Number of original variables for each block.

OUTPUT BlockIP::output

private

Type of output.

const int BlockIP::OUTPUT_FREQ = 1

staticprivate

int BlockIP::output_freq

private

Output information lines will be printed each output_freq IP iterations.

double * BlockIP::p_cg

private

Auxiliary vectors for PCG.

void* BlockIP::params

User parameters to perform objective function calculations.

constexpr double BlockIP::PCG_TOL_LIN = 1.0e-2

staticprivate

constexpr double BlockIP::PCG_TOL_QUAD = 1.0e-3

staticprivate

double BlockIP::pcgtol

private

Tolerance for conjugate gradient at current iteration.

double BlockIP::psi

private

Reduction of dx*dz vector in the rhs of corrector system (in predictor-corrector)

double* BlockIP::qcost

Quadratic cost of variables including slacks.

double BlockIP::qreg

private

Regularization term.

double * BlockIP::r

private

double* BlockIP::r_cg

private

double* BlockIP::rb

private

double * BlockIP::rc

private

constexpr double BlockIP::RED_PCGTOL = 0.95

staticprivate

double BlockIP::red_pcgtol

private

Reduction of pcg_tol at each IP iteration.

constexpr double BlockIP::RHO = 0.995

staticprivate

double BlockIP::rho

private

Reduction of the step-length for the primal and dual variables at each IP iteration.

double* BlockIP::rhs

Upper constraint limits, including linking constraints limits.

double * BlockIP::rhsdy

private

double * BlockIP::rhspcg

private

double * BlockIP::rsw

private

residuals of KKT conditions

double * BlockIP::rxz

private

double * BlockIP::s

private

primal and dual optimization variables

bool BlockIP::sameL

Define if the same matrix is used for each L block.

bool BlockIP::sameN

Define if the same matrix is used for each N block.

bool BlockIP::show_specrad

private

constexpr double BlockIP::SIGMA = 0.1

staticprivate

double BlockIP::sigma

private

Reduction of the centrality parameter at each IP iteration.

StdForm* BlockIP::stdForm

Contains all the information to perform conversions between the original and standard problem.

double* BlockIP::Theta

private

Theta diagonal matrix.

MatrixBlockIP* BlockIP::Theta0

Stores Theta0 preconditioner (Theta_(k+1)) of active linking constraints)

TYPE_DIRECTION BlockIP::this_type_direction

private

Type of direction (Newton or second order) of this particular iteration.

int BlockIP::totcgit

private

Overall number of CG iterations.

TYPE_COMP_DY BlockIP::type_comp_dy

private

Whether the estimation of spectral radius has to be computed.

How dy direction is computed

const TYPE_COMP_DY BlockIP::TYPE_COMP_DY_DEFAULT = CHOL_PWRS_PCG

staticprivate

TYPE_DIRECTION BlockIP::type_direction

private

Type of direction (Newton, second order...) given by user.

const TYPE_DIRECTION BlockIP::TYPE_DIRECTION_DEFAULT = AUTOMATIC

staticprivate

TYPE_PROBLEM BlockIP::type_objective

Type of objective function, linear, quadratic or non-linear.

TYPE_REG BlockIP::type_reg

private

How starting point will be computed.

const TYPE_REG BlockIP::TYPE_REG_DEFAULT = NO_REG

staticprivate

TYPE_START_POINT BlockIP::type_start_point

private

const TYPE_START_POINT BlockIP::TYPE_START_POINT_DEFAULT = SIMPLE

staticprivate

double* BlockIP::ub

Upper bounds, including slack bounds.

double* BlockIP::valpha

private

string* BlockIP::varNames

Variable names including slacks.

double * BlockIP::vbeta

private

double * BlockIP::w

private

WHO_PERMUTES BlockIP::whoperm

private

for matrix factorizations

double* BlockIP::x

private

initial value for regularization

double * BlockIP::y

private

double * BlockIP::z

private

double * BlockIP::z_cg

private

---

The documentation for this class was generated from the following files:

• /home/jcastro2/intpoint/BlockIP/BlockIP/src_blockip-v2/BlockIP.h
• /home/jcastro2/intpoint/BlockIP/BlockIP/src_blockip-v2/BlockIP.C
• /home/jcastro2/intpoint/BlockIP/BlockIP/src_blockip-v2/BlockIPminimize.C