BCLS: Bound Constrained Least Squares

Version 0.1

bcls.c

Go to the documentation of this file.
00001 /* bcls.c
00002    $Revision: 282 $ $Date: 2006-12-17 17:38:00 -0800 (Sun, 17 Dec 2006) $
00003 
00004    ---------------------------------------------------------------------
00005    This file is part of BCLS (Bound-Constrained Least Squares).
00006 
00007    Copyright (C) 2006 Michael P. Friedlander, Department of Computer
00008    Science, University of British Columbia, Canada. All rights
00009    reserved. E-mail: <mpf@cs.ubc.ca>.
00010 
00011    BCLS is free software; you can redistribute it and/or modify it
00012    under the terms of the GNU Lesser General Public License as
00013    published by the Free Software Foundation; either version 2.1 of the
00014    License, or (at your option) any later version.
00015 
00016    BCLS is distributed in the hope that it will be useful, but WITHOUT
00017    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
00018    or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
00019    Public License for more details.
00020 
00021    You should have received a copy of the GNU Lesser General Public
00022    License along with BCLS; if not, write to the Free Software
00023    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
00024    USA
00025    ---------------------------------------------------------------------
00026 */
00032 #ifdef __APPLE__
00033   #include <vecLib/vecLib.h>
00034 #else
00035   #include <cblas.h>
00036 #endif
00037 #include <stdlib.h>
00038 #include <stdio.h>
00039 #include <float.h>
00040 #include <math.h>
00041 #include <assert.h>
00042 #include <setjmp.h>
00043 
00044 #include "bcls.h"
00045 #include "bclib.h"
00046 #include "bcsolver.h"
00047 #include "bcversion.h"
00048 
00063 static void *
00064 xmalloc (int len, size_t size, char * who)
00065 {
00066     register void *value = malloc(len * size);
00067     if (value == NULL)
00068         fprintf( stderr, "Not enough memory to allocate %s.\n", who );
00069     return value;
00070 }
00071 
00092 BCLS *
00093 bcls_create_prob( int mmax, int nmax )
00094 {
00095     int mnmax = imax( nmax, mmax );
00096 
00097     assert( mmax > 0 && nmax > 0 );
00098     
00099     // Allocate the problem.
00100     BCLS *ls = (BCLS *)xmalloc(1, sizeof(BCLS), "ls" );
00101 
00102     // Check if the problem has been successfully allocated.
00103     if (ls == NULL) {
00104         fprintf( stderr, "XXX Could not allocate a BCLS problem.\n");
00105         return ls;
00106     }
00107 
00108     // -----------------------------------------------------------------
00109     // Allocate workspace vectors large enough to accomdate the problem.
00110     // -----------------------------------------------------------------
00111 
00112     // Record the maximum alloted workspace.
00113     ls->mmax = mmax;
00114     ls->nmax = nmax;
00115 
00116     // Residual: r(m+n).
00117     ls->r = (double *)xmalloc( mmax+nmax, sizeof(double), "r" );
00118     if (ls->r == NULL) goto error;
00119 
00120     // Gradient: g(n).
00121     ls->g = (double *)xmalloc( nmax, sizeof(double), "g" );
00122     if (ls->g == NULL) goto error;
00123 
00124     // Search direction, full space: dx(n).
00125     ls->dx = (double *)xmalloc( nmax, sizeof(double), "dx" );
00126     if (ls->dx == NULL) goto error;
00127 
00128     // Search direction, subspace: dxFree(n).
00129     ls->dxFree = (double *)xmalloc( nmax, sizeof(double), "dxFree" );
00130     if (ls->dxFree == NULL) goto error;
00131 
00132     // Step to each breakpoint: aBreak(n).
00133     ls->aBreak = (double *)xmalloc( nmax, sizeof(double), "aBreak" );
00134     if (ls->aBreak == NULL) goto error;
00135 
00136     // Indices of each breakpoint: iBreak(n).
00137     ls->iBreak = (int *)xmalloc( nmax, sizeof(int), "iBreak" );
00138     if (ls->iBreak == NULL) goto error;
00139 
00140     // Variable indices: ix(n).
00141     ls->ix = (int *)xmalloc( nmax, sizeof(int), "ix" );
00142     if (ls->ix == NULL) goto error;
00143 
00144     // Workspace: wrk_v( max(n, m) ).
00145     ls->wrk_u = (double *)xmalloc( mnmax, sizeof(double), "wrk_u" );
00146     if (ls->wrk_u == NULL) goto error;
00147 
00148     // Workspace: wrk_v( max(n, m) ).
00149     ls->wrk_v = (double *)xmalloc( mnmax, sizeof(double), "wrk_v" );
00150     if (ls->wrk_v == NULL) goto error;
00151 
00152     // Workspace: wrk_w( max(n, m) ).
00153     ls->wrk_w = (double *)xmalloc( mnmax, sizeof(double), "wrk_w" );
00154     if (ls->wrk_w == NULL) goto error;
00155 
00156     // -----------------------------------------------------------------
00157     // Initialize this new problem instance.
00158     // -----------------------------------------------------------------
00159     bcls_init_prob( ls );
00160 
00161     // -----------------------------------------------------------------
00162     // Exits.
00163     // -----------------------------------------------------------------
00164 
00165     // Successfull exit.
00166     return ls;
00167 
00168  error:
00169     // Unsuccessful exit.
00170     bcls_free_prob( ls );
00171     return ls;
00172 }
00173 
00185 void
00186 bcls_init_prob( BCLS *ls )
00187 {
00188     // Some quick error checking.
00189     assert( ls->mmax > 0 && ls->nmax > 0 );
00190 
00191     // Check if the problem has been successfully allocated.
00192     if (ls == NULL) {
00193         fprintf( stderr, "XXX The BCLS problem is NULL.\n");
00194         return;
00195     }
00196     
00197     // Initialize the problem structure.
00198     ls->print_info    =  NULL;
00199     ls->print_hook    =  NULL;
00200     ls->fault_info    =  NULL;
00201     ls->fault_hook    =  NULL;
00202     ls->Aprod         =  NULL;
00203     ls->Usolve        =  NULL;
00204     ls->CallBack      =  NULL;
00205     ls->UsrWrk        =  NULL;
00206     ls->anorm         =  NULL;
00207     ls->print_level   =  1;
00208     ls->proj_search   =  BCLS_PROJ_SEARCH_EXACT;
00209     ls->newton_step   =  BCLS_NEWTON_STEP_LSQR;
00210     ls->minor_file    =  NULL;
00211     ls->itnMaj        =  0;
00212     ls->itnMajLim     =  5  * (ls->nmax);
00213     ls->itnMin        =  0;
00214     ls->itnMinLim     =  10 * (ls->nmax);
00215     ls->nAprodT       =  0;
00216     ls->nAprodF       =  0;
00217     ls->nAprod1       =  0;
00218     ls->nUsolve       =  0;
00219     ls->m             =  0;
00220     ls->n             =  0;
00221     ls->unconstrained =  0;
00222     ls->damp          =  0.0;
00223     ls->damp_min      =  1.0e-4;
00224     ls->exit          =  BCLS_EXIT_UNDEF;
00225     ls->soln_rNorm    = -1.0;
00226     ls->soln_dInf     =  0.0;
00227     ls->soln_jInf     = -1;
00228     ls->soln_stat     =  BCLS_SOLN_UNDEF;
00229     ls->optTol        =  1.0e-6;
00230     ls->conlim        =  1.0 / ( 10.0 * sqrt( DBL_EPSILON ) );
00231     ls->mu            =  1.0e-2;
00232     ls->backtrack     =  1.0e-1;
00233     ls->backtrack_limit = 10;
00234 
00235     // Initialize timers.
00236     bcls_timer( &(ls->stopwatch[BCLS_TIMER_TOTAL] ),  BCLS_TIMER_INIT );
00237     ls->stopwatch[BCLS_TIMER_TOTAL].name = "Total time";
00238 
00239     bcls_timer( &(ls->stopwatch[BCLS_TIMER_APROD] ),  BCLS_TIMER_INIT );
00240     ls->stopwatch[BCLS_TIMER_APROD].name = "Total time for Aprod";
00241 
00242     bcls_timer( &(ls->stopwatch[BCLS_TIMER_USOLVE] ),  BCLS_TIMER_INIT );
00243     ls->stopwatch[BCLS_TIMER_USOLVE].name = "Total time for Usolve";
00244 
00245     bcls_timer( &(ls->stopwatch[BCLS_TIMER_LSQR] ),  BCLS_TIMER_INIT );
00246     ls->stopwatch[BCLS_TIMER_LSQR].name = "Total time for LSQR";
00247 
00248     // Initialize constants.
00249     ls->eps         = DBL_EPSILON;
00250     ls->eps2        = pow( DBL_EPSILON, 0.50 );
00251     ls->eps3        = pow( DBL_EPSILON, 0.75 );
00252     ls->epsx        = ls->eps3;
00253     ls->epsfixed    = DBL_EPSILON;
00254     ls->BigNum      = BCLS_INFINITY;
00255 
00256     return;
00257 }
00258 
00268 int
00269 bcls_free_prob( BCLS *ls )
00270 {
00271     // Report an error if it's already NULL.
00272     if (ls == NULL) return 1;
00273 
00274     // Deallocate the workspace.
00275     assert( ls->r       != NULL ); free( ls->r       );
00276     assert( ls->g       != NULL ); free( ls->g       );
00277     assert( ls->dx      != NULL ); free( ls->dx      );
00278     assert( ls->dxFree  != NULL ); free( ls->dxFree  );
00279     assert( ls->aBreak  != NULL ); free( ls->aBreak  );
00280     assert( ls->iBreak  != NULL ); free( ls->iBreak  );
00281     assert( ls->ix      != NULL ); free( ls->ix      );
00282     assert( ls->wrk_u   != NULL ); free( ls->wrk_u   );
00283     assert( ls->wrk_v   != NULL ); free( ls->wrk_v   );
00284     assert( ls->wrk_w   != NULL ); free( ls->wrk_w   );
00285     
00286     // Deallocate the BCLS problem.
00287     free( ls );
00288     
00289     return 0;
00290 }
00291 
00317 void
00318 bcls_set_print_hook( BCLS *ls,
00319                      void *info,
00320                      int (*hook)(void *info, char *msg))
00321 {
00322     ls->print_info = info;
00323     ls->print_hook = hook;
00324     return;
00325 }
00326 
00352 void
00353 bcls_set_fault_hook( BCLS *ls,
00354                      void *info,
00355                      int (*hook)(void *info, char *msg))
00356 {
00357     ls->fault_info = info;
00358     ls->fault_hook = hook;
00359     return;
00360 }
00361 
00373 void
00374 bcls_set_anorm( BCLS *ls,
00375                 double anorm[] )
00376 {
00377     ls->anorm = anorm;
00378     return;
00379 }
00380 
00402 void
00403 bcls_set_usolve( BCLS *ls,
00404                  int (*Usolve)( int mode, int m, int n, int nix,
00405                                 int ix[], double v[], double w[],
00406                                 void *UsrWrk ) )
00407 {
00408     assert( Usolve != NULL );
00409     ls->Usolve = Usolve;
00410     return;
00411 }
00412 
00413 
00436 int
00437 bcls_compute_anorm( BCLS *ls, int n, int m,
00438                     int (*Aprod)
00439                     ( int mode, int m, int n, int nix,
00440                       int ix[], double x[], double y[], void *UsrWrk ),
00441                     void *UsrWrk,
00442                     double anorm[] )
00443 {    
00444     // Make sure that the problem instance has been created.
00445     assert( ls    != NULL );
00446     assert( anorm != NULL );
00447 
00448     int j;
00449     int err;
00450     const int mode = 1;         // Always:  aj = A*ej.
00451     const int nix  = 1;         // Only one column is ever needed.
00452     int    *ix = ls->ix;
00453     double *e  = ls->wrk_u;
00454     double *aj = ls->wrk_v;
00455 
00456     bcls_dload( n, 1.0, e, 1 );
00457     
00458     for (j = 0; j < n; j++) {
00459 
00460         // Multiply A times the jth unit vector.
00461         ix[0] = j;
00462         
00463         // aj <- A * ej.
00464         err = Aprod( mode, m, n, nix, ix, e, aj, UsrWrk );
00465 
00466         // Exit if Aprod returned an error.
00467         if (err) break;
00468 
00469         // Compute the norm of aj and store it in anorm.
00470         anorm[j] = cblas_dnrm2( m, aj, 1 );
00471 
00472         // Make sure that the column norm is too small.
00473         anorm[j] = fmax( BCLS_MIN_COLUMN_NORM, anorm[j] );
00474     }
00475 
00476     return err;
00477 }
00478 
00500 void
00501 bcls_set_problem_data( BCLS *ls, int m, int n,
00502                        int (*Aprod)( int mode, int m, int n, int nix,
00503                                      int ix[], double x[], double y[],
00504                                      void *UsrWrk ),
00505                        void *UsrWrk, double damp,
00506                        double x[], double b[], double c[],
00507                        double bl[], double bu[] )
00508 {
00509     assert( m <= ls->mmax );
00510     assert( n <= ls->nmax );
00511 
00512     assert( Aprod  != NULL ); ls->Aprod  = Aprod;
00513     assert( m      >= 0    ); ls->m      = m;
00514     assert( n      >= 0    ); ls->n      = n;
00515                               ls->UsrWrk = UsrWrk;
00516     assert( damp   >= 0.0  ); ls->damp   = damp;
00517     assert( x      != NULL ); ls->x      = x;
00518     assert( b      != NULL ); ls->b      = b;
00519                               ls->c      = c;
00520     assert( bl     != NULL ); ls->bl     = bl;
00521     assert( bu     != NULL ); ls->bu     = bu;
00522 
00523     return;
00524 }
00525 
00535 char *
00536 bcls_exit_msg( int flag )
00537 {
00538     char *msg;
00539 
00540     if      (flag==BCLS_EXIT_CNVGD)  msg="Optimal solution found";
00541     else if (flag==BCLS_EXIT_MAJOR)  msg="Too many major iterations";
00542     else if (flag==BCLS_EXIT_MINOR)  msg="Too many minor iterations";
00543     else if (flag==BCLS_EXIT_UNBND)  msg="Found direction of infinite descent";
00544     else if (flag==BCLS_EXIT_INFEA)  msg="Bounds are inconsistent";
00545     else if (flag==BCLS_EXIT_APROD)  msg="Aprod requested immediate exit";
00546     else if (flag==BCLS_EXIT_USOLVE) msg="Usolve requested immediate exit";
00547     else if (flag==BCLS_EXIT_CALLBK) msg="CallBack requested immediate exit";
00548     else                             msg="Undefined exit";
00549 
00550     return msg;
00551 }
00552 
00567 int
00568 bcls_solve_prob( BCLS *ls )
00569 {
00570     int err;
00571     int jpInf;
00572     const int preconditioning = ls->Usolve != NULL;
00573     double pInf, timeTot;
00574     double bNorm;  // Norm of the RHS (or 1.0, which ever is bigger).
00575     BCLS_timer watch;
00576 
00577     // Print the banner.
00578     PRINT1(" ----------------------------------------------------\n");
00579     PRINT1(" BCLS -- Bound-Constrained Least Squares, Version %s\n",
00580            bcls_version_info() );
00581     PRINT1(" Compiled on %s\n",  bcls_compilation_info() );
00582     PRINT1(" ----------------------------------------------------\n");
00583 
00584     // -----------------------------------------------------------------
00585     // Print parameters and diagnostic information.
00586     // -----------------------------------------------------------------
00587     bcls_print_params( ls );
00588 
00589     // Examine the bounds and print a summary about them.
00590     // Also check if the problem is feasible!
00591     err = bcls_examine_bnds( ls, ls->n, ls->bl, ls->bu );
00592     if (err) {
00593         ls->exit = err;
00594         goto direct_exit;
00595     }
00596 
00597     // Set the return for longjmp.  First call to setjmp returns 0.
00598     err = setjmp(ls->jmp_env);
00599     if (err) {
00600         ls->exit = err;
00601         goto direct_exit;
00602     }
00603 
00604     // Compute and print some statistics about the column scales.
00605     bcls_examine_column_scales( ls, ls->anorm );
00606 
00607     // Print a log header.
00608     PRINT1("\n %5s  %5s    %9s  %9s  %11s  %9s %7s %7s %7s\n",
00609            "Major","Minor","Residual","xNorm","Optimal",
00610            ls->newton_step == BCLS_NEWTON_STEP_LSQR
00611            ? "LSQR" : "CGLS",
00612            "nSteps", "Free", "OptFac");
00613 
00614     // -----------------------------------------------------------------
00615     // Call the BCLS solver.
00616     // -----------------------------------------------------------------
00617     bcls_timer( &(ls->stopwatch[BCLS_TIMER_TOTAL]), BCLS_TIMER_START );
00618 
00619     bcls_solver( ls, ls->m, ls->n, &bNorm,
00620                  ls->x, ls->b, ls->c, ls->bl, ls->bu, ls->r, ls->g,
00621                  ls->dx, ls->dxFree, ls->ix,
00622                  ls->aBreak, ls->iBreak, ls->anorm );
00623 
00624     bcls_timer( &(ls->stopwatch[BCLS_TIMER_TOTAL]), BCLS_TIMER_STOP );
00625 
00626 
00627     // =================================================================
00628     // Exits.
00629     // =================================================================
00630  direct_exit:
00631     // -----------------------------------------------------------------
00632     // Check feasibility and print a solution summary (iterations, etc.)
00633     // -----------------------------------------------------------------
00634     bcls_primal_inf( ls->n, ls->x, ls->bl, ls->bu, &pInf, &jpInf );
00635 
00636     PRINT1("\n");
00637     PRINT1(" BCLS Exit %4d -- %s\n",ls->exit, bcls_exit_msg(ls->exit));
00638     PRINT1("\n");
00639     PRINT1(" No. of iterations      %8d %7s", ls->itnMin, "");
00640     PRINT1(" Objective value       %17.9e\n", ls->soln_obj  );
00641     PRINT1(" No. of major iterations%8d %7s", ls->itnMaj, "");
00642     PRINT1(" Optimality           (%7d) %8.1e\n",
00643            ls->soln_jInf,ls->soln_dInf);
00644     PRINT1(" No. of calls to Aprod  %8d", ls->nAprodT);
00645 
00646     if (ls->proj_search == BCLS_PROJ_SEARCH_EXACT)
00647         PRINT1(" (%5d)", ls->nAprod1 );
00648     else
00649         PRINT1("  %5s ", "");
00650     PRINT1(" Norm of RHS            %16.1e\n", bNorm);
00651     if (pInf > ls->eps)
00652         PRINT1(" Feasibility          (%7d) %7.1e\n", jpInf, pInf);
00653     if (preconditioning)
00654         PRINT1(" No. of calls to Usolve %8d\n", ls->nUsolve);
00655     PRINT1("\n");
00656 
00657     // -----------------------------------------------------------------
00658     // Print timing statistics.
00659     // -----------------------------------------------------------------
00660     timeTot = fmax(ls->eps,ls->stopwatch[BCLS_TIMER_TOTAL].total);
00661 
00662     watch   = ls->stopwatch[BCLS_TIMER_TOTAL];
00663     timeTot = fmax(1e-6, watch.total); // Safeguard against timeTot = 0.
00664 
00665     PRINT1(" %-25s %5.1f (%4.2f) secs\n",
00666            watch.name, watch.total, 1.0 );
00667 
00668     watch = ls->stopwatch[BCLS_TIMER_LSQR];
00669     PRINT1(" %-25s %5.1f (%4.2f) secs\n",
00670            watch.name, watch.total, watch.total / timeTot );
00671 
00672     watch = ls->stopwatch[BCLS_TIMER_APROD];
00673     PRINT1(" %-25s %5.1f (%4.2f) secs\n",
00674            watch.name, watch.total, watch.total / timeTot );
00675 
00676     // Only print time for Usolve is user provided this routine.
00677     if (preconditioning) {
00678         watch = ls->stopwatch[BCLS_TIMER_USOLVE];
00679         PRINT1(" %-25s %5.1f (%4.2f) secs\n",
00680                watch.name, watch.total, watch.total / timeTot );
00681     }
00682 
00683     return (ls->exit + err);
00684 }

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