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