BCLS: Bound Constrained Least Squares

Version 0.1

bclib.c File Reference

#include <math.h>
#include <stdarg.h>
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <assert.h>
#include <float.h>
#include <setjmp.h>
#include "bclib.h"

Include dependency graph for bclib.c:

Go to the source code of this file.

Functions

void bcls_project_step (int n, double s[], double step, double x[], double dx[], double bl[], double bu[])
 Compute a projected step.
int bcls_examine_bnds (BCLS *ls, const int n, double bl[], double bu[])
 Examine the bounds bl and bu.
void bcls_print_params (BCLS *ls)
 Print the BCLS parameters.
void bcls_examine_column_scales (BCLS *ls, double anorm[])
 Print some stats about the column scales provided by the user.
void bcls_print (BCLS *ls, int current_print_level, char *fmt,...)
 Send output to the user-supplied pretty-printer.
void bcls_fault (BCLS *ls, char *fmt,...)
 Send output to the user-supplied error-printer and then exit.
int bcls_primal_inf (int n, double x[], double bl[], double bu[], double *pInf, int *jpInf)
 Maximum primal infeasibility.
void bcls_dual_inf (int n, double x[], double z[], double bl[], double bu[], double *dInf, int *jInf)
 Maximum dual infeasibility.
void bcls_mid (int n, double x[], double bl[], double bu[])
 Project the point x into the feasible box.
void bcls_aprod (BCLS *ls, int mode, int nix, int ix[], double x[], double y[])
 Interface to the user's matrix-vector product routine.
void bcls_usolve (BCLS *ls, int mode, int nix, int ix[], double v[], double w[])
 Interface to the user's preconditioning routine.
int bcls_callback (BCLS *ls)
 Interface to the user's callback routine.
int bcls_free_vars (BCLS *ls, int n, int ix[], double x[], double g[], double bl[], double bu[])
 Determine the set of free variables.
int bcls_heap_del_min (int numElems, double x[], int ix[])
 Discard the smallest element and contract the heap.
void bcls_heap_sift (int root, int lastChild, double x[], int ix[])
 Perform the "sift" operation for the heap-sort algorithm.
void bcls_heap_build (int n, double x[], int ix[])
 Build a heap by adding one element at a time.
double bcls_vec_dnorm2 (int nix, int ix[], double x[])
 Return the 2-norm of x(ix), where ix is a set of indices.
void bcls_dload (const int n, const double alpha, double x[], const int incx)
 Load a constant into a vector.
void bcls_gather (const int nix, const int ix[], const double x[], double y[])
 Gather the components of x(ix) into y, eg, y = x(ix).
void bcls_scatter (const int nix, const int ix[], const double x[], double y[])
 Scatter the components of x into y(ix), eg, y(ix) = x.


Detailed Description

Helper routines for BCLS.

Definition in file bclib.c.


Function Documentation

void bcls_aprod ( BCLS ls,
int  mode,
int  nix,
int  ix[],
double  x[],
double  y[] 
)

Interface to the user's matrix-vector product routine.

This routine calls the user's Mat-vec routine, passed in as Aprod. We could call Aprod directly, but this interface keeps track of timings and call counts.

Parameters:
[in,out] ls The BCLS solver context.
[in] mode Determines if returning y = A x or x = A'y.
[in] nix Length of ix, x, y.
[in] ix Indices of contributing variables.
[in,out] x RHS or LHS (depending on the mode).
[in,out] y RHS or LHS (depending on the mode).

Definition at line 453 of file bclib.c.

References BCLS::Aprod, BCLS_EXIT_APROD, BCLS_PROD_INIT, BCLS_PROD_TERM, bcls_timer(), BCLS_TIMER_APROD, BCLS_TIMER_START, BCLS_TIMER_STOP, BCLS::jmp_env, BCLS::m, BCLS::n, BCLS::nAprod1, BCLS::nAprodF, BCLS::nAprodT, BCLS::stopwatch, and BCLS::UsrWrk.

Referenced by aprod_free(), aprod_free_lsqr(), bcls_newton_step(), bcls_proj_backtrack(), and bcls_solver().

Here is the call graph for this function:

int bcls_callback ( BCLS ls  ) 

Interface to the user's callback routine.

This routine calls the user's callback routine. We could call CallBack directly, this this interface keeps tracks of various odds and ends.

Parameters:
[in,out] ls The BCLS solver context.
Returns:
Returns (unchanged) the return from the user's callback function.

Definition at line 533 of file bclib.c.

References BCLS::CallBack, and BCLS::UsrWrk.

Referenced by bcls_solver().

void bcls_dload ( const int  n,
const double  alpha,
double  x[],
const int  incx 
)

Load a constant into a vector.

Loads a constant into every component of a vector x. The (usual) case incx = 1, alpha = 0.0 is treated specially.

Parameters:
[in] n The length of x.
[in] alpha The constant double to be loaded into x.
[in,out] x The vector to be loaded.
[in] incx Fill every "incx" elements.

Definition at line 736 of file bclib.c.

Referenced by bcls_cgls(), bcls_compute_anorm(), and bcls_solver().

void bcls_dual_inf ( int  n,
double  x[],
double  z[],
double  bl[],
double  bu[],
double *  dInf,
int *  jInf 
)

Maximum dual infeasibility.

Compute the dual infeasibility (i.e., the optimality) at the point x, measured as the maximum departure from complementarity (with respect to the gradient). If xj has a lower bound lj , the dual infeasibility is (xj - lj ) * zj, but to allow for infinite bounds, we use min( xj - lj, one ) * zj.

Parameters:
[in] n The length of x, z, bl, bu.
[in] x The current point.
[in] z The gradient at current point.
[in] bl The lower bounds on x.
[in] bu The upper bounds on x.
[out] dInf The maximum dual infeasibility.
[out] jInf The variable with the maximum dual infeasibility.

Definition at line 383 of file bclib.c.

Referenced by bcls_solver().

int bcls_examine_bnds ( BCLS ls,
const int  n,
double  bl[],
double  bu[] 
)

Examine the bounds bl and bu.

Print some statistics about the problem bounds, and make sure that the problem is feasible, i.e., that bl <= bu.

Parameters:
[in] ls BCLS solver context.
[in] n Length of bl and bu.
[in] bl Lower bounds on x.
[in] bu Upper bounds on x.
Returns:

Definition at line 100 of file bclib.c.

References BCLS_EXIT_INFEA, BCLS::BigNum, BCLS::epsfixed, PRINT1, and BCLS::unconstrained.

Referenced by bcls_solve_prob().

void bcls_examine_column_scales ( BCLS ls,
double  anorm[] 
)

Print some stats about the column scales provided by the user.

Parameters:
[in] ls BCLS solver context.
[in] anorm Norms of the columns of A.

Definition at line 212 of file bclib.c.

References BCLS::n, PRINT2, and BCLS::print_level.

Referenced by bcls_solve_prob().

void bcls_fault ( BCLS ls,
char *  fmt,
  ... 
)

Send output to the user-supplied error-printer and then exit.

This routine prints an error message specified by the format control string fmt and the optional parameter list. It then terminates execution of the program. See bcls_set_fault_hook.

Parameters:
[in] ls BCLS solver context.
[in] fmt String
[in] ... Any number of arguments.

Definition at line 291 of file bclib.c.

References BCLS::fault_hook, and BCLS::fault_info.

int bcls_free_vars ( BCLS ls,
int  n,
int  ix[],
double  x[],
double  g[],
double  bl[],
double  bu[] 
)

Determine the set of free variables.

Construct an index set (stored in ix[]) of the free (i.e., "inactive") variables. The j-th variable is considered free if it is further than epsx from one of its boundaries, or if its gradient points away from that bound.

Parameters:
[in] ls BCLS solver context.
[in] n The length of x, b, bl, bu.
[out] ix The index set of "free" variables.
[in] x The current iterate.
[in] g The gradient at x.
[in] bl The lower bounds.
[in] bu The upper bounds.
Returns:
Returns the number of free variables. This is therefore the length of the constructed array ix[].

Todo:
The parameter eps and epsx are used in a rather arbitrary way to define if g is "zero" or if x is "close" to its bound. Need to understand how sensitive the algorithm is to this parameter.

Definition at line 569 of file bclib.c.

References BCLS::eps, and BCLS::epsx.

Referenced by bcls_solver().

void bcls_gather ( const int  nix,
const int  ix[],
const double  x[],
double  y[] 
)

Gather the components of x(ix) into y, eg, y = x(ix).

Parameters:
[in] nix The length of the index set.
[in] ix The index set.
[in] x The vector to be gathered.
[out] y The gathered vector.

Definition at line 762 of file bclib.c.

Referenced by aprod_free(), and bcls_newton_step_cgls().

void bcls_heap_build ( int  n,
double  x[],
int  ix[] 
)

Build a heap by adding one element at a time.

Parameters:
[in] n The length of x and ix.
[in,out] x The array to be heapified.
[in,out] ix The indices of the array elements.

Definition at line 686 of file bclib.c.

References bcls_heap_sift().

Here is the call graph for this function:

int bcls_heap_del_min ( int  numElems,
double  x[],
int  ix[] 
)

Discard the smallest element and contract the heap.

On entry, the numElems of the heap are stored in x[0],...,x[numElems-1], and the smallest element is x[0]. The following operations are performed:

  1. Swap the first and last elements of the heap
  2. Shorten the length of the heap by one.
  3. Restore the heap property to the contracted heap. This effectively makes x[0] the next smallest element in the list.

Parameters:
[in] numElems The number of elements in the current heap.
[in,out] x The array to be modified.
[in,out] ix The indices of the array.
Returns:
The number of elements in the heap after it has been contracted.

Definition at line 619 of file bclib.c.

References bcls_dswap, bcls_heap_sift(), and bcls_iswap.

Here is the call graph for this function:

void bcls_heap_sift ( int  root,
int  lastChild,
double  x[],
int  ix[] 
)

Perform the "sift" operation for the heap-sort algorithm.

A heap is a collection of items arranged in a binary tree. Each child node is less than or equal to its parent. If x[k] is the parent, than its children are x[2k+1] and x[2k+2].

This routine promotes ("sifts up") children that are smaller than their parents. Thus, this is a "reverse" heap, where the smallest element of the heap is the root node.

Parameters:
[in] root The root index from which to start sifting.
[in] lastChild The last child (largest node index) in the sift operation.
[in,out] x The array to be sifted.
[in,out] ix The indices of the array.

Definition at line 657 of file bclib.c.

References bcls_dswap, and bcls_iswap.

Referenced by bcls_heap_build(), and bcls_heap_del_min().

void bcls_mid ( int  n,
double  x[],
double  bl[],
double  bu[] 
)

Project the point x into the feasible box.

Project the point x into the box defined by the bounds bl and bu.

Parameters:
[in] n The length of x, bl, bu.
[in,out] x The current point.
[in] bl The lower bounds on x.
[in] bu The upper bounds on x.

Definition at line 425 of file bclib.c.

Referenced by bcls_solver().

int bcls_primal_inf ( int  n,
double  x[],
double  bl[],
double  bu[],
double *  pInf,
int *  jpInf 
)

Maximum primal infeasibility.

Compute the primal infeasibility at the point x. This is measured as the maximum violation of the bounds.

Parameters:
[in] n Length of x, bl, bu
[in] x Current point.
[in] bl Lower bounds.
[in] bu Upper bounds.
[out] pInf Maximum infeasibility.
[out] jpInf Variable that achieves the maximum infeasibility.
Returns:
  • 0 if x is feasible to within machine precision (DBL_EPSILON).
  • 1 otherwise.

Definition at line 334 of file bclib.c.

Referenced by bcls_solve_prob().

void bcls_print ( BCLS ls,
int  current_print_level,
char *  fmt,
  ... 
)

Send output to the user-supplied pretty-printer.

This routine prints an informative message specified by the format control string fmt and the optional parameter list. See bcls_set_print_hook.

Parameters:
[in] ls BCLS solver context.
[in] current_print_level Print level of the current call.
[in] fmt String
[in] ... Any number of arguments.

Definition at line 253 of file bclib.c.

References BCLS::print_hook, BCLS::print_info, and BCLS::print_level.

void bcls_print_params ( BCLS ls  ) 

Print the BCLS parameters.

Print the user-settable BCLS parameters. Print level must be >= 2 for anything to print.

Parameters:
[in] ls BCLS solver context.

Definition at line 175 of file bclib.c.

References BCLS::anorm, BCLS_NEWTON_STEP_LSQR, BCLS_PROJ_SEARCH_APPROX, BCLS::c, BCLS::damp, BCLS::itnMajLim, BCLS::itnMinLim, BCLS::m, BCLS::n, BCLS::newton_step, BCLS::optTol, PRINT2, BCLS::proj_search, and BCLS::Usolve.

Referenced by bcls_solve_prob().

void bcls_project_step ( int  n,
double  s[],
double  step,
double  x[],
double  dx[],
double  bl[],
double  bu[] 
)

Compute a projected step.

Given a current point x and a step dx, compute the projected step s = P[x + step*dx] - x.

Parameters:
[in] n Length of x, s, dx, bl, bu.
[out] s The projected step.
[in] step Step length.
[in] x The current point.
[in] dx The step that needs to be projected.
[in] bl Lower bounds on x.
[in] bu Upper bounds on x.

Definition at line 60 of file bclib.c.

Referenced by bcls_proj_backtrack().

void bcls_scatter ( const int  nix,
const int  ix[],
const double  x[],
double  y[] 
)

Scatter the components of x into y(ix), eg, y(ix) = x.

Parameters:
[in] nix The length of the index set.
[in] ix The index set.
[in] x The vector to be scattered.
[out] y The scattered vector.

Definition at line 782 of file bclib.c.

Referenced by aprod_free(), and bcls_solver().

void bcls_usolve ( BCLS ls,
int  mode,
int  nix,
int  ix[],
double  v[],
double  w[] 
)

Interface to the user's preconditioning routine.

This routine calls the user's preconditioning routine, passed in as Usolve. We could call Usolve directly, this this interface keeps tracks of timings and call counts.

Parameters:
[in,out] ls The BCLS solver context.
[in] mode Determines if solving U v = w or U'w = v.
[in] nix Length of ix, v, w.
[in] ix Indices of contributing variables.
[in,out] v The solution or RHS (depending on the mode).
[in,out] w The solution or RHS (depending on the mode).

Definition at line 497 of file bclib.c.

References BCLS_EXIT_USOLVE, BCLS_PRECON_U, BCLS_PRECON_Ut, bcls_timer(), BCLS_TIMER_START, BCLS_TIMER_STOP, BCLS_TIMER_USOLVE, BCLS::jmp_env, BCLS::m, BCLS::n, BCLS::nUsolve, BCLS::stopwatch, BCLS::Usolve, and BCLS::UsrWrk.

Referenced by aprod_free_lsqr(), and bcls_newton_step().

Here is the call graph for this function:

double bcls_vec_dnorm2 ( int  nix,
int  ix[],
double  x[] 
)

Return the 2-norm of x(ix), where ix is a set of indices.

Parameters:
[in] nix Length of ix.
[in] ix Set of indices to access in x.
[in] x Input vector.
Returns:
The 2-norm of the subvector.

Definition at line 708 of file bclib.c.

Generated on Sun Mar 4 22:50:03 2007 by Doxygen 1.5.1