BCLS: Bound Constrained Least Squares

Version 0.1

bclib.c

Go to the documentation of this file.

/* bclib.c
   $Revision: 285 $ $Date: 2006-12-20 09:33:50 -0800 (Wed, 20 Dec 2006) $

   ----------------------------------------------------------------------
   This file is part of BCLS (Bound-Constrained Least Squares).

   Copyright (C) 2006 Michael P. Friedlander, Department of Computer
   Science, University of British Columbia, Canada. All rights
   reserved. E-mail: <mpf@cs.ubc.ca>.
   
   BCLS is free software; you can redistribute it and/or modify it
   under the terms of the GNU Lesser General Public License as
   published by the Free Software Foundation; either version 2.1 of the
   License, or (at your option) any later version.
   
   BCLS is distributed in the hope that it will be useful, but WITHOUT
   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
   or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
   Public License for more details.
   
   You should have received a copy of the GNU Lesser General Public
   License along with BCLS; if not, write to the Free Software
   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
   USA
   ----------------------------------------------------------------------
*/
#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"

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

    for (j = 0; j < n; j++) {

        xpj = x[j] + step*dx[j];

        if ( xpj < bl[j] )
            s[j] = bl[j] - x[j];
        
        else if ( xpj > bu[j] )
            s[j] = bu[j] - x[j];
        
        else
            s[j] = step * dx[j];
    }
    return;
}

int
bcls_examine_bnds( BCLS *ls, const int n, double bl[], double bu[] )
{
    int j, jopen;
    double blj, buj;
    const double bigL = - ls->BigNum;
    const double bigU =   ls->BigNum;
    const double epsfixed = ls->epsfixed;

    // Variable-type counts.  No. of variables that
    int neql = 0; // are "open" (have an interior)
    int nfix = 0; // are fixed
    int nneg = 0; // are unbounded below
    int npos = 0; // are unbounded above
    int nlow = 0; // have a finite lower bound
    int nupp = 0; // have a finite upper bound
    int nnrm = 0; // have finite lower and upper bounds
    int nfre = 0; // are free.

    for (j = 0; j < n; j++) {
        jopen = 0;
        blj   = bl[j];
        buj   = bu[j];
        
        if (buj < blj)                  // Infeasible interval.  Exit.
            goto infeas;
        
        else if (buj - blj <= epsfixed) // Fixed variable.
            nfix++;

        else {                          // Open interval.
            neql++;
            jopen = 1;
        }

        if     (blj <= bigL  &&  buj == 0            ) nneg++;
        if     (blj == 0     &&  buj >= bigU         ) npos++;
        if     (blj >  bigL                  && jopen) nlow++;
        if     (                 buj <  bigU && jopen) nupp++;
        if     (blj >  bigL  &&  buj <  bigU         ) nnrm++;
        if     (blj <= bigL  &&  buj >= bigU         ) nfre++;
    }
    
    ls->unconstrained = nfre == n;

    PRINT1("\n Bounds:");
    if (ls->unconstrained) {
        PRINT1(" problem is unconstrained\n");
    }
    else {
        PRINT1("\n   [-inf,0]    [0,inf]  Finite bl  Finite bu  "
               " Two bnds      Fixed       Free\n");
        PRINT1("  %9d  %9d  %9d  %9d  %9d  %9d  %9d\n",
               nneg, npos, nlow, nupp, nnrm, nfix, nfre );
    }

    return 0;

 infeas:
    PRINT1(" Infeasible problem: bl[%d] = %10.2e > %10.2e bu[%d]\n",
           j, blj, buj, j);

    return BCLS_EXIT_INFEA;
}

void
bcls_print_params( BCLS *ls )
{
    PRINT2("\n");
    PRINT2(" Parameters:\n");
    PRINT2("  Rows............. %8d\t",   ls->m);
    PRINT2("  Columns.......... %8d\n",   ls->n);
    PRINT2("  Optimality tol... %8.1e\t", ls->optTol);
    PRINT2("  Major itn limit.. %8d\n",   ls->itnMajLim);
    PRINT2("  Damp............. %8.1e\t", ls->damp);
    PRINT2("  Minor itn limit.. %8d\n",   ls->itnMinLim);
    PRINT2("  Linsearch type... %8s\t", 
           ls->proj_search == BCLS_PROJ_SEARCH_APPROX
           ? "approx" : "exact");
    PRINT2("  Linear term...... %8s\n",
           ls->c == NULL ? "no" : "yes" );
    PRINT2("  Newton step...... %8s\t",
           ls->newton_step == BCLS_NEWTON_STEP_LSQR
           ? "lsqr" : "cgls");
    PRINT2("  Gradient step.... %8s\n",
           ls->anorm == NULL
           ? "unscaled" : "scaled");
    PRINT2("  Preconditoner.... %8s\n",
           ls->Usolve == NULL
           ? "no" : "yes");

    return;
}

void
bcls_examine_column_scales( BCLS *ls, double anorm[] )
{
    int j;
    int    minj;     // Index of the column with the minimum col scale.
    int    maxj;     // ............................ maximum col scale.
    double minnorm;  // The value of the minimum             col scale.
    double maxnorm;  // ........................             col scale.

    if ( anorm == NULL ) return;
    if ( ls->print_level < 2 ) return;

    minnorm = 1.0e20;
    maxnorm = 0.0;

    for (j = 0; j < ls->n; j++) {
        if ( minnorm > anorm[j] ) { minnorm = anorm[j]; minj = j; };
        if ( maxnorm < anorm[j] ) { maxnorm = anorm[j]; maxj = j; };
    }
    PRINT2( "\n" );
    PRINT2( " Column scales:  %6s   %7s\n", "column", "norm" );
    PRINT2( "       minimum   %6d   %7.1e\n", minj, minnorm  );
    PRINT2( "       maximum   %6d   %7.1e\n", maxj, maxnorm  );

    return;
}

void
bcls_print( BCLS *ls, int current_print_level, char *fmt, ...)
{ 
    va_list arg;
    char msg[4095+1];

    // Exit immediately if the current print level is too low.
    if ( ls->print_level < current_print_level )
        return;

    // Format the message.
    va_start(arg, fmt);
    vsprintf(msg, fmt, arg);
    assert(strlen(msg) <= 4095);
    va_end(arg);
    
    // Pass the message to the user-defined hook routine.
    if (ls->print_hook != NULL && ls->print_hook(ls->print_info, msg) == 0)
        return;

    // Othwerise, send the message to the standard output */
      fprintf(stdout, "%s", msg);

      return;
}

void
bcls_fault( BCLS *ls, char *fmt, ...)
{ 
    va_list arg;
    char msg[4095+1];

    // Format the message.
    va_start(arg, fmt);
    vsprintf(msg, fmt, arg);
    assert(strlen(msg) <= 4095);
    va_end(arg);
    
    // Pass the message to the user-defined hook routine.
    if (ls->fault_hook != NULL &&
        ls->fault_hook(ls->fault_info, msg) == 0) goto skip;

    // Othwerise, send the message to the standard output */
      fprintf(stderr, "%s\n", msg);

 skip:
      // return to the calling program.
      exit(EXIT_FAILURE);     
}

int
bcls_primal_inf( int n, double x[], double bl[], double bu[],
                 double *pInf, int *jpInf )
{
    int j;
    double pj = 0;
    double blj, buj, xj;

    *pInf  = 0.0;
    *jpInf = -1;
    
    for (j = 0; j < n; j++) {
        blj = bl[j];
        buj = bu[j];
        xj  =  x[j];
        if ( xj > buj )
            pj = xj - buj;
        else if ( blj > xj )
            pj = blj - xj;
        
        if (*pInf  < pj ) {
            *pInf  = pj;
            *jpInf = j;
        }
        
    }

    return *pInf > DBL_EPSILON;
}

void
bcls_dual_inf( int n, double x[], double z[], double bl[], double bu[],
               double *dInf, int *jInf )
{
    int    j;
    double dj, blj, buj;

    *dInf = 0.0;
    *jInf = 0;
   
    for (j = 0; j < n; j++) {
        blj = bl[j];
        buj = bu[j];
        if (blj < buj) {
            dj  = z[j];
            if      (dj > 0)
                dj =   dj * fmin( x[j] - blj, 1.0 );
            else if (dj < 0)
                dj = - dj * fmin( buj - x[j], 1.0 );
     
            if (*dInf  < dj) {
                *dInf  = dj;
                *jInf  = j;
            }
        }
    }

    return;
}

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

    for (j = 0; j < n; j++ ) {
        if      ( x[j] < bl[j] ) x[j] = bl[j];
        else if ( x[j] > bu[j] ) x[j] = bu[j];
    }
    return;
}

void
bcls_aprod( BCLS *ls, int mode, int nix,
            int ix[], double x[], double y[] )
{    
    // Increase the mat-vec counters, unless this is an initialization
    // or de-initialization call.
    if (mode == BCLS_PROD_INIT ||
        mode == BCLS_PROD_TERM )
        ;// Relax.
    else {
        if      (nix == 1    ) ls->nAprod1++;
        else if (nix  < ls->n) ls->nAprodF++;
                               ls->nAprodT++;
    }

    // Call the user's routine.
    bcls_timer( &(ls->stopwatch[BCLS_TIMER_APROD]), BCLS_TIMER_START );

    int err = ls->Aprod( mode, ls->m, ls->n, nix, ix, x, y, ls->UsrWrk );

    bcls_timer( &(ls->stopwatch[BCLS_TIMER_APROD]), BCLS_TIMER_STOP );

    // Return to top of BCLS's stack if error is signaled.
    if (err)
        longjmp(ls->jmp_env, BCLS_EXIT_APROD);
    return;
}

void
bcls_usolve( BCLS *ls, int mode, int nix,
             int ix[], double v[], double w[] )
{    
    // Increase the counters only for actual solves.
    if (mode == BCLS_PRECON_U  ||  mode == BCLS_PRECON_Ut )
        ls->nUsolve++;

    // Call the user's routine.
    bcls_timer( &(ls->stopwatch[BCLS_TIMER_USOLVE]), BCLS_TIMER_START );

    int err = ls->Usolve( mode, ls->m, ls->n, nix, ix, v, w, ls->UsrWrk );

    bcls_timer( &(ls->stopwatch[BCLS_TIMER_USOLVE]), BCLS_TIMER_STOP );

    // Return to top of BCLS's stack if error is signaled.
    if (err)
        longjmp(ls->jmp_env, BCLS_EXIT_USOLVE);
    return;
}

int
bcls_callback( BCLS *ls )
{    
    if ( ls->CallBack )
        return ls->CallBack( ls, ls->UsrWrk );
    else
        return 0;
}

int
bcls_free_vars( BCLS *ls, int n, int ix[], double x[], 
                double g[], double bl[], double bu[] )
{
    int j;
    int nFree = 0;
    int low, upp, inc, dec;
    const double epsx = ls->epsx;
    const double eps  = ls->eps;

    for (j = 0; j < n; j++) {

        low  = x[j]  - bl[j]  <  epsx;
        upp  = bu[j] -  x[j]  <  epsx;

        inc  = g[j]  <  eps;
        dec  = g[j]  >  eps;

        if      ( low )
            ; // Relax -- lower bound binding.
        else if ( upp )
            ; // Relax -- upper bound binding.
        else {
              // Var is free to move.
            ix[nFree] = j;
            nFree++;
        }
    }
    return nFree;
}

int
bcls_heap_del_min( int numElems, double x[], int ix[] )
{
    assert(numElems > 0);
    
    int lastChild = numElems - 1;

    // Swap the smallest element with the lastChild.
    bcls_dswap(  x[0],  x[lastChild] );
    bcls_iswap( ix[0], ix[lastChild] );

    // Contract the heap size, thereby discarding the smallest element.
    lastChild--;
    
    // Restore the heap property of the contracted heap.
    bcls_heap_sift( 0, lastChild, x, ix );

    return numElems - 1;
}

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

    for (; (child = (root * 2) + 1) <= lastChild; root = child) {

        if (child < lastChild)
            if ( x[child] > x[child+1] )
                child++;
        
        if ( x[child] >= x[root] )
            break;

        bcls_iswap( ix[root], ix[child] );
        bcls_dswap(  x[root],  x[child] );
    }
    return;
}

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

    for (i = n/2; i >= 0; i--)
        bcls_heap_sift( i, n-1, x, ix );

    return;
}

double
bcls_vec_dnorm2( int nix, int ix[], double x[] )
{    
    int j, k;
    double norm = 0.0;
    
    for (j = 0; j < nix; j++) {
        k     = ix[j];
        norm +=  x[k] * x[k];
    }
    norm = sqrt(norm);
    
    return norm;
}

void
bcls_dload( const int n, const double alpha, double x[], const int incx )
{    
    int i, ix;

    if (incx <= 0) return;

    ix = 0;

    for (i = 0; i < n; i++) {
        x[ix]  = alpha;
        ix    += incx;
    }
    return;
}

void
bcls_gather( const int nix, const int ix[], const double x[], double y[] )
{
    int j, k;
    for (j = 0; j < nix; j++) {
        k = ix[ j ];
        y[ j ] = x[ k ];
    }
}

void
bcls_scatter( const int nix, const int ix[], const double x[], double y[] )
{
    int j, k;
    for (j = 0; j < nix; j++) {
        k = ix[ j ];
        y[ k ] = x[ j ];
    }
}

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