Version 0.1
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