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

#include <BlockIP.h>

Collaboration diagram for BlockIP:

Classes
struct	BackupLnk

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.

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.

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_iTheta_iL_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 *	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_compute_direction ()

void	min_Newton_direction ()

void	min_second_order_predictor_corrector_direction ()

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
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	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_REG	type_reg
	How starting point will be computed.

double	factor_reg
	Type of regularization performed.

double *	x
	initial value for regularization

double *	s
	primal and dual optimization variables

double *	rsw
	residuals of KKT conditions

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

double	normrc
	norms of rhs, costs, rb, and rc

double	alpha_d
	step lengths

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

double *	dw
	Primal and dual optimization directions.

double	gap
	Duality gap.

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 *	Cvaux2
	Auxiliary vectors for interior-point algorithm.

double *	p_cg
	Auxiliary vectors for PCG.

double	qreg
	Regularization term.

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	numWithUb
	Number of variables with upper bound (not infinity)

Detailed Description

Main class for loading and solving problems.

Constructor & Destructor Documentation

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.

Member Function Documentation

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

Here is the caller graph for this function:

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. !

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_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

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.

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).

Here is the caller graph for this function:

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

Here is the caller graph for this function:

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

Here is the call graph for this function:

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

Here is the call graph for this function:

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

Here is the call graph for this function:

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 (AThetaA)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 AThetaA' 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 (AThetaA)dy= rhsdy
rhsdy[1:km+l],:	on input right-hand-side of (AThetaA)dy= rhsdy
only_solve,:	if true, only numeric solve (neither symbolic nor numeric factorization are needed) since a previous call with the same AThetaA' 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_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_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.

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

Here is the call graph for this function:

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_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)):

Predictor direction

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

Compute sigma and psi

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:

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

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_qreg ( )

private

Updates regularization term

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. !

Here is the caller graph for this 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. !

Here is the call graph for this function:

Here is the caller graph for this function:

void BlockIP::set_defaults_minimize ( )

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

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::write_BlockIP_format ( const char * filename )

Write the problem in BlockIP format file.

Parameters

filename File name without extension

Here is the caller graph for this function:

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

Here is the call graph for this function:

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

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

double BlockIP::inf

private

Infinity value.

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

MatrixBlockIP* BlockIP::Theta0

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

TYPE_COMP_DY BlockIP::type_comp_dy

private

Whether the estimation of spectral radius has to be computed.

How dy direction is computed

WHO_PERMUTES BlockIP::whoperm

private

for matrix factorizations

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

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

Classes

Public Member Functions

Public Attributes

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation