1 /* 2 * This file is part of MPSolve 3.2.1 3 * 4 * Copyright (C) 2001-2020, Dipartimento di Matematica "L. Tonelli", Pisa. 5 * License: http://www.gnu.org/licenses/gpl.html GPL version 3 or higher 6 * 7 * Authors: 8 * Leonardo Robol <leonardo.robol@unipi.it> 9 */ 10 11 #ifndef MPS_CONTEXT_H 12 #define MPS_CONTEXT_H 13 14 #include <mps/mps.h> 15 #include <pthread.h> 16 17 #ifdef __cplusplus 18 extern "C" 19 { 20 #endif 21 22 /** 23 * @file 24 * @brief This file contains the definition of <code>mps_context</code> and 25 * most of its fields. 26 */ 27 28 /** 29 * @brief Routine that performs the computation loop to solve the polynomial 30 * or the secular equation 31 */ 32 typedef void (*mps_mpsolve_ptr)(mps_context *status); 33 34 /** 35 * @brief Pointer to the callback for the async version of mpsolve 36 */ 37 typedef void* (*mps_callback)(mps_context * status, void * user_data); 38 39 40 /* 41 * Macros for casting user functions 42 */ 43 #define MPS_FNEWTON_PTR(x) (mps_fnewton_ptr) & (x) 44 #define MPS_DNEWTON_PTR(x) (mps_dnewton_ptr) & (x) 45 #define MPS_MNEWTON_PTR(x) (mps_mnewton_ptr) & (x) 46 #define MPS_CHECK_DATA_PTR(x) (mps_check_data_ptr) & (x) 47 #define MPS_FSTART_PTR(x) (mps_fstart_ptr) & (x) 48 #define MPS_DSTART_PTR(x) (mps_dstart_ptr) & (x) 49 #define MPS_MPSOLVE_PTR(x) (mps_mpsolve_ptr) & (x) 50 51 #ifdef _MPS_PRIVATE 52 /** 53 * @brief this struct holds the state of the mps computation 54 */ 55 struct mps_context { 56 /** 57 * @brief true if an error has occurred during the computation. 58 */ 59 mps_boolean error_state; 60 61 /** 62 * @brief The text describing the last error occurred. 63 */ 64 char * last_error; 65 66 /** 67 * @brief This value is non NULL if mpsolve is launched via mps_mpsolve_async() 68 * and in that case holds a pointer to the thread pool used to manage 69 * asynchronous callbacks. 70 * 71 * It will be automatically freed by mps_free_data(). 72 */ 73 mps_thread_pool * self_thread_pool; 74 75 /** 76 * @brief true if we are trying to resume previously interrupted. 77 * 78 * Not yet implemented. 79 */ 80 mps_boolean resume; 81 82 /** 83 * @brief True if check of radius should be performed at the end 84 * of the algorithm. 85 * 86 * Only works for algorithm MPS_ALGORITHM_SECULAR_GA. 87 */ 88 mps_boolean chkrad; 89 90 /** 91 * @brief Callback called when the async version of mps_mpsolve(), i.e. 92 * terminate the computation. 93 */ 94 mps_callback callback; 95 96 /** 97 * @brief Pointer to user_data passed to the callback. 98 */ 99 void * user_data; 100 101 /** 102 * @brief The operation running now. Can be used to debug what's happening 103 * event if mpsolve was launched without debug enabled. 104 */ 105 mps_operation operation; 106 107 /** 108 * @brief This value is true if the data for the computation has been allocated 109 * by calling mps_allocate_data(). It is used by mps_free_data() to know what has 110 * to be freed. 111 */ 112 mps_boolean initialized; 113 114 /** 115 * @brief Bytes containing the flags of debug enabled. 116 */ 117 mps_debug_level debug_level; 118 119 /** 120 * @brief True if the computation has reached the maximum allowed precision. 121 */ 122 mps_boolean over_max; 123 124 /** 125 * @brief Configuration of the input of MPSolve 126 */ 127 mps_input_configuration * input_config; 128 129 /** 130 * @brief Output configuration for MPSolve. 131 */ 132 mps_output_configuration * output_config; 133 134 /** 135 * @brief Newton isolation of the cluster. 136 */ 137 int newtis; 138 139 /** 140 * @brief Old value for the newton isolation of the cluster. 141 */ 142 int newtis_old; 143 144 /* 145 * INPUT / OUTPUT STREAMS 146 */ 147 148 /** 149 * @brief <code>true</code> if log are needed. They will 150 * be written to <code>logstr</code> 151 * 152 * @see logstr 153 */ 154 mps_boolean DOLOG; 155 156 /** 157 * @brief <code>true</code> if warning are needed. 158 */ 159 mps_boolean DOWARN; 160 161 /** 162 * @brief <code>true</code> if root sorting is desired. It will 163 * be performed with routines in <code>mps_sort.c</code>. 164 */ 165 mps_boolean DOSORT; 166 167 /** 168 * @brief Default input stream. 169 */ 170 FILE *instr; 171 172 /** 173 * @brief Default output stream. 174 */ 175 FILE *outstr; 176 177 /** 178 * @brief Default log stream 179 */ 180 FILE *logstr; 181 182 /** 183 * @brief Stream used to resume an interrupted computation or to load 184 * the approximations from a custom file. 185 */ 186 FILE *rtstr; 187 188 /* 189 * CONSTANT, PARAMETERS 190 */ 191 192 /** 193 * @brief number of max packets of iterations 194 */ 195 int max_pack; 196 197 /** 198 * @brief number of max iterations per packet 199 */ 200 int max_it; 201 202 /** 203 * @brief Number of max newton iterations for gravity center 204 * computations. 205 */ 206 int max_newt_it; 207 208 /** 209 * @brief Maximum allowed number of bits for mp numbers: used in 210 * high precision shift. 211 */ 212 long int mpwp_max; 213 214 /** 215 * @brief Maximum precision reached during the computation. 216 */ 217 mps_long_int_mt data_prec_max; 218 219 /** 220 * @brief Precision operation give best results when done one 221 * thread at a time :) 222 */ 223 pthread_mutex_t precision_mutex; 224 225 /** 226 * @brief True if this is the first iteration after the precision has been 227 * raised. 228 */ 229 mps_boolean just_raised_precision; 230 231 /** 232 * @brief mps_boolean value that determine if we should 233 * use a random seed for starting points 234 */ 235 mps_boolean random_seed; 236 237 /* 238 * POLYNOMIAL DATA: SHARED VARIABLES 239 */ 240 241 /** 242 * @brief degree of zero-deflated polynomial. 243 */ 244 int n; 245 246 /** 247 * @brief input degree and allocation size. 248 */ 249 int deg; 250 251 /* Solution related variables */ 252 253 /** 254 * @brief Last computing phase. 255 */ 256 mps_phase lastphase; 257 258 /** 259 * @brief Selected starting case, can be 'd' for DPE 260 * or 'f' for floating point 261 */ 262 mps_phase starting_case; 263 264 /** 265 * @brief Set to true if the approximation are the best that 266 * can be obtained with the current precision 267 */ 268 mps_boolean best_approx; 269 270 /** 271 * @brief shift in the angle in the positioning of the 272 * starting approximation for the last cluster. It will 273 * be used to determine the new sigma to maximize distance 274 * between starting points. 275 */ 276 double last_sigma; 277 278 /** 279 * @brief Vector containing count of in, out and uncertaing roots. 280 */ 281 int count[3]; 282 283 /** 284 * @brief Number of zero roots. 285 */ 286 int zero_roots; 287 288 /** 289 * @brief Output index order 290 */ 291 int *order; 292 293 /** 294 * @brief Vector of points to the 295 * current root approximations. 296 */ 297 mps_approximation ** root; 298 299 /** 300 * @brief <code>true</code> if the float phase should be skipped, 301 * passing directly do dpe phase. 302 */ 303 mps_boolean skip_float; 304 305 /** 306 * @brief Input precision of the coefficients. 307 */ 308 rdpe_t eps_in; 309 310 /** 311 * @brief Output precision of the roots. 312 */ 313 rdpe_t eps_out; 314 315 /** 316 * @brief Logarithm of the max modulus of the coefficients. 317 */ 318 double lmax_coeff; 319 320 /** 321 * @brief Bits of working precision that mpsolve is using. 322 */ 323 long int mpwp; 324 325 /** 326 * @brief Current multiprecision epsilon. 327 */ 328 rdpe_t mp_epsilon; 329 330 /** 331 * @brief Log of the lower bound to the minumum distance of the roots 332 */ 333 double sep; 334 335 /** 336 * @brief Clusterization object that represent the clusterization 337 * detected on the roots. 338 * 339 * This value is updated with the <code>mps_*cluster</code> 340 * routines. 341 * 342 * @see mps_fcluster(), mps_dcluster(), mps_mcluster() 343 */ 344 mps_clusterization * clusterization; 345 346 /** 347 * @brief Standard complex coefficients of the polynomial. 348 * 349 * This is used as a temporary vector while shifting the polynomial 350 * with a new gravity center in <code>mps_fshift()</code>. 351 */ 352 cplx_t *fppc1; 353 354 /** 355 * @brief <code>dpe</code> complex coefficients of the polynomial. 356 * 357 * This is used as a temporary vector while shifting the polynomial 358 * with a new gravity center in <code>mps_dshift()</code>. 359 */ 360 cdpe_t *dpc1; 361 362 /** 363 * @brief <code>dpe</code> complex coefficients of the polynomial. 364 * 365 * This is used as a temporary vector while shifting the polynomial 366 * with a new gravity center in <code>mps_dshift()</code>. 367 */ 368 cdpe_t *dpc2; 369 370 /** 371 * @brief Multiprecision complex coefficients of the polynomial. 372 * 373 * This is used as a temporary vector while shifting the polynomial 374 * with a new gravity center in <code>mps_mshift()</code>. 375 */ 376 mpc_t *mfpc1; 377 378 /** 379 * @brief Multiprecision complex coefficients of the 380 * first derivative of the polynomial. 381 * 382 * This is used as a temporary vector while shifting the polynomial 383 * with a new gravity center in <code>mps_mshift()</code>. 384 */ 385 mpc_t *mfppc1; 386 387 /** 388 * @brief Vector representing sparsity of the polynomial in the 389 * same way that <code>spar</code> does. 390 * 391 * It is used as a temporary vector. 392 * 393 * @see spar 394 */ 395 mps_boolean *spar1; 396 397 /** 398 * @brief Old value of <code>punt</code> (temporary vector). 399 * 400 * @see punt 401 */ 402 int *oldpunt; 403 404 /** 405 * @brief Vector containing the moduli of the coefficients 406 * of the polynomial as floating point numbers. 407 * 408 * It is used in the computation of the newton polygonal in 409 * <code>mps_fcompute_starting_radii()</code>. 410 * 411 * @see mps_fcompute_starting_radii() 412 */ 413 double *fap1; 414 415 /** 416 * @brief Vector containing the logarithm of the moduli of 417 * the coefficients of the polynomial as floating 418 * point numbers. 419 * 420 * It is used in the computation of the newton polygonal in 421 * <code>mps_fcompute_starting_radii()</code>. 422 * 423 * @see mps_fcompute_starting_radii() 424 * @see fap1 425 */ 426 double *fap2; 427 428 /** 429 * @brief Vector containing the moduli of the coefficients 430 * of the polynomial as <code>dpe</code> numbers. 431 * 432 * It is used in the computation of the newton polygonal in 433 * <code>mps_dcompute_starting_radii()</code>. 434 * 435 * @see mps_dcompute_starting_radii() 436 */ 437 rdpe_t *dap1; 438 439 /** 440 * @brief Temporary vector containing the old value of 441 * <code>again</code>. 442 * 443 * @see again 444 */ 445 mps_boolean *again_old; 446 447 /* SECTION -- Algorihtm selection */ 448 449 /** 450 * @brief This is used in the program to switch behavious based 451 * on the algorithm that is been used now. 452 */ 453 mps_algorithm algorithm; 454 455 /** 456 * @brief Strategy used to dispose the starting approximations. 457 */ 458 mps_starting_strategy starting_strategy; 459 460 /** 461 * @brief Routine that performs the loop needed to coordinate 462 * root finding. It has to be called to do the hard work. 463 */ 464 void (*mpsolve_ptr)(mps_context *status); 465 466 /** 467 * @brief This is the polynomial that is currently being solved in MPSolve. 468 */ 469 mps_polynomial * active_poly; 470 471 /** 472 * @brief Pointer to the secular equation used in computations. 473 */ 474 mps_secular_equation * secular_equation; 475 476 /** 477 * @brief Number of threads to be spawned. 478 */ 479 int n_threads; 480 481 /** 482 * @brief The thread pool used for the concurrent part of MPSolve. 483 */ 484 mps_thread_pool * pool; 485 486 /** 487 * @brief Auxiliary memory used in regeneation to avoid thread-safeness 488 * issues. 489 */ 490 mpc_t * bmpc; 491 492 /** 493 * @brief True if Jacobi-style iterations must be used in the secular 494 * algorithm. 495 */ 496 mps_boolean jacobi_iterations; 497 498 /** 499 * @brief Char to be intersted after the with statement in the output piped to gnuplot. 500 */ 501 const char * gnuplot_format; 502 503 /* DEBUG SECTIONS */ 504 505 unsigned long int regeneration_time; 506 unsigned long int mp_iteration_time; 507 unsigned long int dpe_iteration_time; 508 unsigned long int fp_iteration_time; 509 510 mps_boolean exit_required; 511 512 long int minimum_gmp_precision; 513 514 /** 515 * @brief In case this field is set to true MPSolve will avoid a multiprecision 516 * phase, and exit instead. 517 * 518 * Note that this may imply that not all the required digits/isolation condition 519 * may have been computed. 520 */ 521 mps_boolean avoid_multiprecision; 522 523 /** 524 * @brief This flags enables the "crude" only approximation mode of MPSolve. 525 * 526 * If this mode is activated MPSolve will only perform a basic Aberth iteration 527 * in floating point and then exit. Note that the output result will still be 528 * guaranteed but in general it will not be possible to reach arbitrary precision 529 * and the results may be quite far from the roots for bad conditioned polynomials. 530 */ 531 mps_boolean crude_approximation_mode; 532 533 /** 534 * @brief This is a pointer to the regeneration driver that performs the standard regeneration 535 * step. MPSolve provides a default implementation of this can be overloaded by 536 * the user. 537 */ 538 mps_regeneration_driver *regeneration_driver; 539 540 }; /* End of typedef struct { ... */ 541 542 #endif /* #ifdef _MPS_PRIVATE */ 543 544 /* Allocator, deallocator, constructors.. */ 545 mps_context * mps_context_new (void); 546 void mps_context_free (mps_context * s); 547 548 /* Accessor functions (setters) */ 549 void mps_context_abort (mps_context * s); 550 int mps_context_set_poly_d (mps_context * s, cplx_t * coeff, 551 long unsigned int n); 552 void mps_context_set_input_poly (mps_context * s, mps_polynomial * p); 553 int mps_context_set_poly_i (mps_context * s, int *coeff, long unsigned int n); 554 void mps_context_select_algorithm (mps_context * s, mps_algorithm algorithm); 555 void mps_context_set_degree (mps_context * s, int n); 556 557 #ifdef _MPS_PRIVATE 558 void mps_context_allocate_poly_inplace (mps_context * s, int n); 559 #endif 560 561 /* Accessor functions */ 562 long int mps_context_get_data_prec_max (mps_context * s); 563 long int mps_context_get_minimum_precision (mps_context * s); 564 int mps_context_get_degree (mps_context * s); 565 int mps_context_get_roots_d (mps_context * s, cplx_t ** roots, double **radius); 566 int mps_context_get_roots_m (mps_context * s, mpc_t ** roots, rdpe_t ** radius); 567 int mps_context_get_zero_roots (mps_context * s); 568 mps_root_status mps_context_get_root_status (mps_context * ctx, int i); 569 mps_boolean mps_context_get_over_max (mps_context * s); 570 mps_polynomial * mps_context_get_active_poly (mps_context * ctx); 571 mps_approximation** mps_context_get_approximations (mps_context * ctx); 572 573 /* I/O options and flags */ 574 void mps_context_set_input_prec (mps_context * s, long int prec); 575 void mps_context_set_output_prec (mps_context * s, long int prec); 576 void mps_context_set_output_format (mps_context * s, mps_output_format format); 577 void mps_context_set_output_goal (mps_context * s, mps_output_goal goal); 578 void mps_context_set_starting_phase (mps_context * s, mps_phase phase); 579 void mps_context_set_log_stream (mps_context * s, FILE * logstr); 580 void mps_context_set_jacobi_iterations (mps_context * s, mps_boolean jacobi_iterations); 581 void mps_context_select_starting_strategy (mps_context * s, mps_starting_strategy strategy); 582 void mps_context_set_avoid_multiprecision (mps_context * s, mps_boolean avoid_multiprecision); 583 void mps_context_set_crude_approximation_mode (mps_context * s, mps_boolean crude_approximation_mode); 584 void mps_context_set_regeneration_driver (mps_context * s, mps_regeneration_driver * rd); 585 586 /* Debugging */ 587 void mps_context_set_debug_level (mps_context * s, mps_debug_level level); 588 void mps_context_add_debug_domain (mps_context * s, mps_debug_level level); 589 590 /* Get input and output config pointers */ 591 mps_input_configuration * mps_context_get_input_config (mps_context * s); 592 mps_output_configuration * mps_context_get_output_config (mps_context * s); 593 594 /* Error handling */ 595 mps_boolean mps_context_has_errors (mps_context * s); 596 char * mps_context_error_msg (mps_context * s); 597 598 599 #ifdef __cplusplus 600 } 601 #endif 602 603 #endif 604