BCLS: Bound Constrained Least Squares

Version 0.1

/examples/bcsol.c

00001 // =====================================================================
00002 // $Revision: 290 $ $Date: 2007-03-04 22:01:54 -0800 (Sun, 04 Mar 2007) $
00003 //
00004 // Driver for BCLS: a bound-constrained least-squares solver.
00005 //
00006 // 25 Aug 05: Original version.
00007 //            Michael P. Friedlander
00008 //            mpf@cs.ubc.ca
00009 //            University of British Columbia
00010 // =====================================================================
00011 
00012 #include <stdlib.h>
00013 #include <stdio.h>
00014 #include <string.h>
00015 #include <assert.h>
00016 #include <getopt.h>
00017 #include "cs.h"    // The CSparse matrix library.
00018 #include "iohb.h"  // The Harwell-Boing file reader library.
00019 #include "bcls.h"  // The BCLS solver library.
00020 
00021 // ---------------------------------------------------------------------
00022 // Global variables.
00023 // ---------------------------------------------------------------------
00024 
00025 // Name of the (required) input file for A and b.
00026 static char *in_file_Ab_name = NULL;
00027 
00028 // Name of the (optional) input file for bl, bu, c, x.
00029 static char *in_file_opt_name = NULL;
00030 
00031 // Name of the (optional) minor-output file.
00032 static char *out_file_minor = NULL;
00033 
00034 // Optional storage description.  At least RHS must be stored.
00035 static char opt_storage[5] = "lu";
00036 
00037 // Minor iteration limit.
00038 static int minor_its = -1;
00039 
00040 // Major iteration limit.
00041 static int major_its = -1;
00042 
00043 // Optimality tolerance.
00044 static double opt_tol = -1.0;
00045 
00046 // Preconditioning.  None by default.
00047 static int preconditioning = 0;
00048 
00049 // Print level.
00050 static int print_level = -1;
00051 
00052 // Type of projected search.
00053 static int proj_search = BCLS_PROJ_SEARCH_EXACT;
00054 
00055 // Method for Newton step computation.
00056 static int newton_step = BCLS_NEWTON_STEP_LSQR;
00057 
00058 // Method for Newton step computation.
00059 static int scaled_steepest = 0;
00060 
00061 // Regularization parameter.
00062 static double damp = 0.0;
00063 
00064 // If this flag is set, only check input, don't solve the problem.
00065 static int check = 0;
00066 
00067 // ---------------------------------------------------------------------
00068 // Workspace structure passed to Aprod and Usolve routines.
00069 // ---------------------------------------------------------------------
00070 typedef struct {
00071     cs *A;
00072 } worksp;
00073 
00074 // ---------------------------------------------------------------------
00075 // xmalloc
00076 // Malloc wrapper.
00077 //
00078 // Returns a pointer to the newly allocated memory, or NULL if malloc
00079 // returned an eror.
00080 // ---------------------------------------------------------------------
00081 static void *
00082 xmalloc(int n, size_t size, char * who) {
00083 
00084     register void *value = malloc(n * size);
00085     if (value == NULL) {
00086         fprintf( stderr, "Not enough memory to allocate %s.\n", who );
00087         exit(EXIT_FAILURE);
00088     }
00089     return value;
00090 }
00091 
00092 // ---------------------------------------------------------------------
00093 // dload
00094 // Load a constant into a vector.
00095 // ---------------------------------------------------------------------
00096 static void
00097 dload( const int n, const double alpha, double x[] ) {
00098 
00099     int j;
00100     for (j = 0; j < n; j++) x[j] = alpha;
00101     return;
00102 }
00103 
00104 // ---------------------------------------------------------------------
00105 // display_help
00106 // Display help at the commandline.
00107 // ---------------------------------------------------------------------
00108 static void
00109 display_help(char *my_name) {
00110 
00111     printf("Usage: %s [options...] filename\n", my_name);
00112     printf("The last argument is the filename of a Harwell-Boeing file\n"
00113            "that contains A and b (1-based indexing).\n");
00114     printf("Options:\n");
00115     printf("   -c, --check        do not solve problem, check input"
00116            " data only\n");
00117     printf("       --cgls         use CGLS for Newton step computation\n");
00118     printf("   -d, --damp         damping (regularization) parameter\n");
00119     printf("   -a, --approx       approximate projected linesearch\n");
00120     printf("   -h, --help         print this message\n");
00121     printf("   -m, --minor        minor (CGLS/LSQR) iteration limit\n");
00122     printf("   -M, --major        major iteration limit\n");
00123     printf("       --minor-out    minor (CGLS/LSQR) output file."
00124            "Use - for stdout.\n");
00125     printf("   -O  --optional     optional data file for l, u, c, x\n");
00126     printf("   -o  --opttol       optimality tolerance\n");
00127     printf("       --precon       use LU preconditioner\n");
00128     printf("   -p, --print        print level 0-6 (default is 1)\n");
00129     printf("       --scaled       scaled steepest descent"
00130            " (default is unscaled)\n");
00131     printf("   -s, --storage      storage description in optional data file"
00132            " (default is 'lu')\n");
00133     return;
00134 }
00135 
00136 // ---------------------------------------------------------------------
00137 // parse_cmdline
00138 // Parse the commandline options.
00139 // ---------------------------------------------------------------------
00140 static void
00141 parse_cmdline(int argc, char *argv[]) {
00142 
00143     int c, j;
00144     int option_index;
00145     static struct option long_options[] = {
00146         // Long options w/o short equivalents.  Only set a flag.
00147         {"cgls",     no_argument, &newton_step,   BCLS_NEWTON_STEP_CGLS},
00148         {"scaled",   no_argument, &scaled_steepest, 1},
00149         {"precon",   no_argument, &preconditioning, 1},
00150         // Long and short options.
00151         {"check",    no_argument,       0, 'c'},
00152         {"damp",     required_argument, 0, 'd'},
00153         {"approx",   no_argument,       0, 'a'},
00154         {"help",     no_argument,       0, 'h'},
00155         {"minor",    required_argument, 0, 'm'},
00156         {"major",    required_argument, 0, 'M'},
00157         {"minor-out",required_argument, 0, 'z'},
00158         {"optional", required_argument, 0, 'O'},
00159         {"opttol",   required_argument, 0, 'o'},
00160         {"print",    required_argument, 0, 'p'},
00161         {"storage",  required_argument, 0, 's'},
00162         // End of option flags.
00163         {0,0,0,0}
00164     };
00165     
00166     while (1) {
00167         
00168         c = getopt_long(argc, argv, "cd:ahm:M:O:o:p:s:",
00169                         long_options, &option_index );
00170         
00171         if ( c == -1 ) // No more options.
00172             break;
00173 
00174         switch (c) {
00175 
00176         case 0:   // Long option with no argument.
00177             break;
00178             
00179         case 'c': // -c  --check
00180             check = 1;
00181             break;
00182 
00183         case 'd': // -d  --damp
00184             damp = atof(optarg);
00185             break;
00186             
00187         case 'a': // -a --approx
00188             proj_search = BCLS_PROJ_SEARCH_APPROX;
00189             break;
00190 
00191         case 'h': // -h  --help
00192             display_help(argv[0]);
00193             exit(EXIT_SUCCESS);
00194             break;
00195 
00196         case 'm': // -m  --minor
00197             minor_its = atoi(optarg);
00198             break;
00199 
00200         case 'M': // -M  --major
00201             major_its = atoi(optarg);
00202             break;
00203 
00204         case 'O': // -O  --optional
00205             if (in_file_opt_name != NULL) {
00206                 printf("Only one optional data file allowed\n");
00207                 exit(EXIT_FAILURE);
00208             }
00209             in_file_opt_name = optarg;
00210             break;
00211 
00212         case 'z': // --minor-out
00213             out_file_minor = optarg;
00214             break;
00215 
00216         case 'o': //     --opttol
00217             opt_tol = atof(optarg);
00218             break;
00219 
00220         case 'p': // -p  --print
00221             print_level = atoi(optarg);
00222             if (print_level < 0) {
00223                 printf("print level must be 0 or positive");
00224                 exit(EXIT_FAILURE);
00225             }
00226             break;
00227 
00228         case 's': // -s  --storage
00229             strncpy( opt_storage, optarg, sizeof(opt_storage) );
00230             for (j = 0; j < strlen(opt_storage); j++) {
00231                 if (opt_storage[j] != 'l' &&
00232                     opt_storage[j] != 'u' &&
00233                     opt_storage[j] != 'c' &&
00234                     opt_storage[j] != 'x' ) {
00235                     printf("Storage descriptor '%s' not recognized.\n",
00236                            &(opt_storage[j]));
00237                     exit(EXIT_FAILURE);
00238                 }
00239             }
00240             break;
00241 
00242         case '?':
00243             // getopt_long already printed an error message.
00244             printf("Try %s --help\n", argv[0]);
00245             exit(EXIT_FAILURE);
00246             break;
00247 
00248         default:
00249             display_help(argv[0]);
00250             exit(EXIT_SUCCESS);
00251             
00252         }
00253     }
00254     
00255     // Gather the remaining non-option arguments.
00256     if (optind < argc) {
00257         while (optind < argc) {
00258             if (in_file_Ab_name != NULL) {
00259                 printf("Only one data input file allowed\n");
00260                 printf("Try %s --help\n", argv[0]);
00261                 exit(EXIT_FAILURE);
00262             }
00263             in_file_Ab_name = argv[optind];
00264             optind++;
00265         }
00266     }
00267     else {
00268         printf("No input file for A, b specified; try %s --help\n", argv[0]);
00269         exit(EXIT_FAILURE);
00270     }
00271 }
00272 
00273 // ---------------------------------------------------------------------
00274 // Aprod
00275 // Matrix-vector products.
00276 //
00277 // If     mode == BCLS_PROD_A  (0), compute y <- A *x, with x untouched;
00278 // and if mode == BCLS_PROD_At (1), compute x <- A'*y, with y untouched.
00279 // ---------------------------------------------------------------------
00280 static int
00281 Aprod( int mode, int m, int n, int nix, int ix[],
00282        double x[], double y[], void *UsrWrk ) {
00283 
00284     int     i, j, k, l;
00285     double  aij, xj, sum;
00286     worksp * Wrk = (worksp *)UsrWrk;
00287     cs *A = (cs *)Wrk->A;
00288 
00289     assert( mode == BCLS_PROD_A     ||
00290             mode == BCLS_PROD_At    ||
00291             mode == BCLS_PROD_INIT  ||
00292             mode == BCLS_PROD_TERM  );
00293     assert( nix  <= n );
00294     assert( A->m == m );
00295     assert( A->n == n );
00296 
00297     if (mode == BCLS_PROD_A) {
00298 
00299         dload( m, 0.0, y );
00300 
00301         for (l = 0; l < nix; l++) {
00302             j = ix[l];
00303             xj = x[j];
00304             if (xj == 0.0)
00305                 ; // Relax.
00306             else
00307                 for (k = A->p[j]; k < A->p[j+1]; k++) {
00308                     aij   = A->x[k];
00309                     i     = A->i[k];
00310                     y[i] += aij * xj;
00311                 }
00312         }
00313     }
00314 
00315     else if (mode == BCLS_PROD_At) {
00316         for (l = 0; l < nix; l++) {
00317             j = ix[l];
00318             sum = 0;
00319             for (k = A->p[j]; k < A->p[j+1]; k++) {
00320                 aij  = A->x[k];
00321                 i    = A->i[k];
00322                 sum += aij * y[i];
00323             }
00324             x[j] = sum;
00325         }
00326     }
00327 
00328     // Exit.
00329     return 0;
00330 }
00331 
00332 // ---------------------------------------------------------------------
00333 // Usolve
00334 // Preconditioner.
00335 //
00336 // If     mode = BCLS_PRECON_U,  solve  U v = w,
00337 // and if mode = BCLS_PRECON_Ut, solve  U'w = v.
00338 // ---------------------------------------------------------------------
00339 static int
00340 Usolve( int mode, int m, int n, int nix, int ix[],
00341         double v[], double w[], void *UsrWrk ){
00342     
00343     assert( nix  <= n );
00344     assert( mode == BCLS_PRECON_U    ||
00345             mode == BCLS_PRECON_Ut   ||
00346             mode == BCLS_PRECON_INIT ||
00347             mode == BCLS_PRECON_TERM );
00348     
00349     if ( mode == BCLS_PRECON_INIT ) {
00350         // -------------------------------------------------------------
00351         // Initialize the preconditioner.
00352         // -------------------------------------------------------------
00353 
00354         // -------------------------------------------------------------
00355         // Factorize A(:,ix).
00356         // -------------------------------------------------------------
00357 
00358 
00359     }
00360     else if ( mode == BCLS_PRECON_U ) {
00361         
00362         // Solve  U v = w.
00363 
00364     }
00365     else if ( mode == BCLS_PRECON_Ut ) {
00366 
00367         // Solve  U'w = v.
00368 
00369     }
00370     else if ( mode == BCLS_PRECON_TERM )
00371         ; // Relax.
00372     
00373     // Exit.
00374     return 0;
00375 }
00376 
00377 // ---------------------------------------------------------------------
00378 // CallBack.
00379 // Periodically called by BCLS to test if the user wants to exit.
00380 // ---------------------------------------------------------------------
00381 static int
00382 CallBack( BCLS *ls, void *UsrWrk )
00383 {
00384     int err;
00385     err = 0;  // No error.
00386     return err;
00387 }
00388 
00389 // ---------------------------------------------------------------------
00390 // pretty_printer
00391 // This is the print-routine that will be used by BCLS for its output.
00392 // ---------------------------------------------------------------------
00393 static int
00394 pretty_printer( void *io_file, char *msg ) {
00395     fprintf( io_file, msg );
00396     return 0;
00397 }
00398 
00399 // ---------------------------------------------------------------------
00400 // main
00401 // ---------------------------------------------------------------------
00402 int
00403 main(int argc, char *argv[]) {
00404 
00405     // Miscellaneous variables.
00406     int i, j, err, nRhs;
00407     double *dPtr;
00408     char line[1025], *blActiv, *buActiv, *MatType;
00409     FILE *in_file;
00410 
00411     // These variables help to define the BCLS problem
00412     BCLS   *ls;              // A BCLS problem.
00413     worksp  Wrk;             // Workspace.
00414     int     m, n, nnz;       // Problem dimensions.
00415     double *c     = NULL;    // Linear vector.
00416     double *x     = NULL;    // Solution vector.
00417     double *b     = NULL;    // RHS vector.
00418     double *bl    = NULL;    // Lower bounds.
00419     double *bu    = NULL;    // Upper bounds.
00420     double *anorm = NULL;    // Column norms of A.
00421 
00422     // -----------------------------------------------------------------
00423     // Basic initialization.
00424     // -----------------------------------------------------------------
00425     parse_cmdline( argc, argv );
00426 
00427     // -----------------------------------------------------------------
00428     // Read in the header for A and b.
00429     // -----------------------------------------------------------------
00430     if (!readHB_info( in_file_Ab_name, &m, &n, &nnz, &MatType, &nRhs ))
00431         exit(EXIT_FAILURE);
00432     free(MatType);
00433     printf( " ----------------------------------------------------\n");
00434     printf( " Reading from data file %s:\n", in_file_Ab_name );
00435     printf( " Matrix size:  m = %i  n = %i  nnz = %i  nRhs = %i\n",
00436             m, n, nnz, nRhs );
00437     printf( " ----------------------------------------------------\n");
00438 
00439     // -----------------------------------------------------------------
00440     // Allocate storage for A, b, bl, bu, x, c.
00441     // -----------------------------------------------------------------
00442 
00443     // Allocate storage for A.
00444     Wrk.A = cs_spalloc( m, n, nnz, 1, 0 );
00445 
00446     // Storage for these vectors are always needed.
00447     b  = (double *)xmalloc( m, sizeof(double), "b"  );
00448     bl = (double *)xmalloc( n, sizeof(double), "bl" );
00449     bu = (double *)xmalloc( n, sizeof(double), "bu" );
00450     x  = (double *)xmalloc( n, sizeof(double), "x"  );
00451 
00452     // Storage for c is optional.
00453     if (strchr(opt_storage, 'c'))
00454     c  = (double *)xmalloc( n, sizeof(double), "c"  );
00455 
00456     // -----------------------------------------------------------------
00457     // Load A, b.
00458     // -----------------------------------------------------------------
00459     readHB_mat_double( in_file_Ab_name, Wrk.A->p, Wrk.A->i, Wrk.A->x );
00460     readHB_aux_double( in_file_Ab_name, 'F', b );
00461 
00462     // -----------------------------------------------------------------
00463     // Load bl, bu, c, x.
00464     // -----------------------------------------------------------------
00465 
00466     // Load bl, bu, x with defaults if no values are specified.
00467     if (!strchr(opt_storage,'l') || !in_file_opt_name) dload(n,-1e20,bl);
00468     if (!strchr(opt_storage,'u') || !in_file_opt_name) dload(n, 1e20,bu);
00469     if (!strchr(opt_storage,'x') || !in_file_opt_name) dload(n,  0.0, x);
00470 
00471     // Open optional file for reading.
00472     if (in_file_opt_name) {
00473         if ( (in_file = fopen( in_file_opt_name, "r")) == NULL) {
00474             fprintf(stderr,"Error: Cannot open file: %s\n",
00475                     in_file_opt_name);
00476             exit(EXIT_FAILURE);
00477         }
00478         
00479         for (j = 0; j < strlen(opt_storage); j++) {
00480             
00481             if      (opt_storage[j] == 'l')  dPtr = bl;
00482             else if (opt_storage[j] == 'u')  dPtr = bu;
00483             else if (opt_storage[j] == 'x')  dPtr = x;  
00484             else                             dPtr = c;
00485             
00486             for (i = 0; i < n; i++ ) {
00487                 if ( fgets( line, sizeof(line), in_file ) == NULL ) {
00488                     fprintf( stderr, "Reached end of file %s before "
00489                              "reading %d entries in mode %s.\n",
00490                              in_file_opt_name, i, &(opt_storage[j]) );
00491                     exit(EXIT_FAILURE);
00492                 }
00493                 sscanf( line, "%lg", &(dPtr[i]) );
00494             }
00495         }
00496         
00497         // Done with input file.
00498         fclose( in_file );
00499     }
00500     // -----------------------------------------------------------------
00501     // Initialize a BCLS problem.  This routine MUST be called before
00502     // any other BCLS routine.
00503     // -----------------------------------------------------------------
00504     ls = bcls_create_prob( m, n );
00505 
00506     // Compute column norms.
00507     if ( scaled_steepest ) {
00508 
00509         // Allocate memory to hold the scales.
00510         anorm = xmalloc( n, sizeof(double), "anorm" );
00511 
00512         // Compute the scales and let BCLS know we have them.
00513         bcls_compute_anorm( ls, n, m, Aprod, &Wrk, anorm );
00514         bcls_set_anorm( ls, anorm );
00515     }
00516 
00517     // Instatiate a particular BCLS problem.
00518     bcls_set_problem_data( ls,    // The BCLS problem
00519                            m,     // Number of problem rows
00520                            n,     // Number of problem columns
00521                            Aprod, // The Mat-vec routine
00522                            &Wrk,  // Arbitrary data for the Mat-vec routine
00523                            damp,  // Damping parameter
00524                            x,     // Solution vector
00525                            b,     // RHS vector
00526                            c,     // Linear term (may be NULL)
00527                            bl,    // Lower-bounds vector
00528                            bu );  // Upper-bounds vector
00529 
00530     // Set the user options.
00531     bcls_set_print_hook( ls, stdout, pretty_printer );
00532     if (major_its   > 0) ls->itnMajLim     = major_its;
00533     if (minor_its   > 0) ls->itnMinLim     = minor_its;
00534     if (opt_tol     > 0) ls->optTol        = opt_tol;
00535     if (print_level > 0) ls->print_level   = print_level;
00536     if (preconditioning) bcls_set_usolve( ls, Usolve );
00537     if (out_file_minor) {
00538         if (strcmp(out_file_minor, "-") == 0)
00539             ls->minor_file = stdout;
00540         else
00541             ls->minor_file = fopen(out_file_minor, "w+");
00542     }
00543     ls->proj_search    = proj_search;
00544     ls->newton_step    = newton_step;
00545     ls->CallBack       = CallBack;
00546                          
00547     // Call the main BCLS routine.
00548     err = bcls_solve_prob( ls );
00549 
00550     // Print the solution if print level >= 4.
00551     if (print_level >= 4) {
00552         printf("\n Solution\n --------\n");
00553         printf("%4s  %18s %1s %18s %1s %18s  %18s\n",
00554                "Var","Lower","","Value","","Upper","Gradient");
00555         for (j = 0; j < n; j++) {
00556             blActiv = "";
00557             buActiv = "";
00558             if (x[j] - bl[j] < ls->epsx) blActiv = "=";
00559             if (bu[j] - x[j] < ls->epsx) buActiv = "=";
00560             printf("%4d  %18.11e %1s %18.11e %1s %18.11e  %18.11e\n",
00561                    j+1, bl[j], blActiv, x[j], buActiv, bu[j], (ls->g)[j]);
00562         }
00563     }
00564     // Deallocate the BCLS problem.
00565     err = bcls_free_prob( ls );
00566 
00567     // Deallocate memory.
00568     cs_spfree( Wrk.A );
00569     if ( b   ) free( b     );
00570     if ( bl  ) free( bl    );
00571     if ( bu  ) free( bu    );
00572     if ( x   ) free( x     );
00573     if ( c   ) free( c     );
00574     if (anorm) free( anorm );
00575 
00576     // Close files.
00577     if (out_file_minor)
00578         fclose( ls->minor_file );
00579 
00580     // Exit with no errors.  Whew!
00581     return (EXIT_SUCCESS);
00582 }

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