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 /** 12 * @file 13 * @brief Header file for secular-related routines. 14 */ 15 16 #ifndef MPS_SECULAR_H_ 17 #define MPS_SECULAR_H_ 18 19 #include <mps/mps.h> 20 #include <float.h> 21 22 MPS_BEGIN_DECLS 23 24 #define MPS_SECULAR_EQUATION(t) (MPS_POLYNOMIAL_CAST (mps_secular_equation, t)) 25 #define MPS_IS_SECULAR_EQUATION(t) (mps_polynomial_check_type (t, "mps_secular_equation")) 26 27 #ifdef _MPS_PRIVATE 28 29 /* CONSTANTS */ 30 31 /** 32 * @brief This is the number of bits used when first passed 33 * in multiprecision. 34 */ 35 #define MPS_SECULAR_STARTING_MP_PRECISION 128 36 37 /** 38 * @brief This is the higher precision supported by GMP that is 39 * lower than the precision supported by the standard floating 40 * point machinery. It is used to set an "equivalent" precision 41 * in s->mpwp for the step of multiprecision coefficient regeneration. 42 */ 43 #define MPS_SECULAR_EQUIVALENT_FP_PRECISION (MPS_SECULAR_STARTING_MP_PRECISION / 2) 44 45 struct mps_secular_equation_double_buffer { 46 char active; 47 mpc_t *ampc1; 48 mpc_t *ampc2; 49 mpc_t *bmpc1; 50 mpc_t *bmpc2; 51 }; 52 53 /** 54 * @brief Secular equation data. 55 * 56 * A secular equation is an equation in the form 57 * \f[ 58 * \sum_{i = 1}^{n} \frac{a_i}{z - b_i} = 1 59 * \f] 60 * and this struct holds the values of the parameters \f$a_i\f$ 61 * and \f$b_i\f$. 62 */ 63 struct mps_secular_equation { 64 struct mps_polynomial __base_class__; 65 66 struct mps_secular_equation_double_buffer db; 67 68 /** 69 * @brief Vector of \f$a_i\f$ as complex floating 70 * point numbers. 71 */ 72 cplx_t *afpc; 73 74 /** 75 * @brief Same as <code>afpc</code>, but the <code>dpe</code> 76 * version. 77 */ 78 cdpe_t *adpc; 79 80 /** 81 * @brief Vector with the values of \f$b_i\f$ as complex 82 * floating point numbers. 83 */ 84 cplx_t *bfpc; 85 86 /** 87 * @brief Same as <code>bfpc</code>, but the <code>dpe</code> 88 * version. 89 */ 90 cdpe_t *bdpc; 91 92 /** 93 * @brief Same as <code>afpc</code>, but the multiprecision 94 * version. 95 */ 96 mpc_t *ampc; 97 98 /** 99 * @brief Mutexes thatn need to be locked to ensure consistent 100 * access to ampc[j] variable. 101 */ 102 pthread_mutex_t * ampc_mutex; 103 104 /** 105 * @brief Same as <code>bfpc</code>, but the multiprecision 106 * version. 107 */ 108 mpc_t *bmpc; 109 110 /** 111 * @brief Mutexes that need to be locked to ensure consistent 112 * access to bmpc[j] variable. 113 */ 114 pthread_mutex_t * bmpc_mutex; 115 116 /** 117 * @brief Moduli of the floating point a_i 118 * coefficients of the secular equation. 119 */ 120 double *aafpc; 121 122 /** 123 * @brief Moduli of the floating point b_i 124 * coefficients of the secular equation. 125 */ 126 double *abfpc; 127 128 /** 129 * @brief DPE Moduli of the CDPE of Multiprecision a_i 130 * coefficients of the secular equation. 131 */ 132 rdpe_t *aadpc; 133 134 /** 135 * @brief DPE Moduli of the CDPE of Multiprecision b_i 136 * coefficients of the secular equation. 137 */ 138 rdpe_t *abdpc; 139 140 /** 141 * @brief Initial multiprecision coefficients saved for latter 142 * regeneration in <code>mps_secular_ga_regenerate_coefficients()</code>. 143 */ 144 mpc_t *initial_ampc; 145 146 /** 147 * @brief Initial multiprecision coefficients saved for latter 148 * regeneration in <code>mps_secular_ga_regenerate_coefficients()</code>. 149 */ 150 mpc_t *initial_bmpc; 151 152 /** 153 * @brief Initial rational coefficients, if rational input is selected. 154 * This value is the real part of the \f$a_i\f$ coefficients. 155 */ 156 mpq_t *initial_ampqrc; 157 158 /** 159 * @brief Initial rational coefficients, if rational input is selected. 160 * This value is the real part of the \f$b_i\f$ coefficients. 161 */ 162 mpq_t *initial_bmpqrc; 163 164 /** 165 * @brief Initial rational coefficients, if rational input is selected. 166 * This value is the imaginary part of the \f$a_i\f$ coefficients. 167 */ 168 mpq_t *initial_ampqic; 169 170 /** 171 * @brief Initial rational coefficients, if rational input is selected. 172 * This value is the imaginary part of the \f$b_i\f$ coefficients. 173 */ 174 mpq_t *initial_bmpqic; 175 176 /** 177 * @brief This mutex is locked while changing precision. 178 */ 179 pthread_mutex_t precision_mutex; 180 }; /* End of struct mps_secular_equation {... */ 181 182 /** 183 * @brief This is a struct that represent an iteration on a root. It contains 184 * information that could be useful for mps_secular_*iterate() routine to determine 185 * some error bound and provide a method for the routine to communicate if 186 * it was able to set the radius or not (by setting the <code>radius_set</code> 187 * in the right way). 188 */ 189 struct mps_secular_iteration_data { 190 /** 191 * @brief The index of the roots on which the iterations 192 * is being carried out. 193 */ 194 long int k; 195 196 /** 197 * @brief The state of the iteration. This is a pointer 198 * to a boolean that tells if the iterator has been 199 * able to set a radius or not, because the radius 200 * that was there before was better. 201 */ 202 mps_boolean radius_set; 203 204 /** 205 * @brief Global mutex used to synchronization, but mainly 206 * while testing new MP implementations. 207 */ 208 pthread_mutex_t * gs_mutex; 209 210 /** 211 * @brief Thread local copy of the \f$a_i\f$ coefficients of the secular 212 * equation. 213 */ 214 mpc_t * local_ampc; 215 216 /** 217 * @brief Thread local copy of the \f$b_i\f$ coefficients of the secular 218 * equation. 219 */ 220 mpc_t * local_bmpc; 221 222 /** 223 * @brief Thread local copy of the floating point coefficients of the secular 224 * equation. 225 */ 226 cplx_t * local_afpc; 227 228 /** 229 * @brief Thread local copy of the floating point coefficients of the secular 230 * equation. 231 */ 232 cplx_t * local_bfpc; 233 234 /** 235 * @brief Thread local copy of the CDPE coefficients of the secular 236 * equation. 237 */ 238 cdpe_t * local_adpc; 239 240 /** 241 * @brief Thread local copy of the CDPE coefficients of the secular 242 * equation. 243 */ 244 cdpe_t * local_bdpc; 245 }; 246 247 #endif /* #ifdef _MPS_PRIVATE */ 248 249 250 /* MACROS */ 251 #define mps_secular_equation_from_status(s) (mps_secular_equation*)(s)->secular_equation 252 253 /* Routines in secular-newton.c */ 254 void mps_secular_fnewton (mps_context * st, mps_polynomial * p, mps_approximation * root, cplx_t corr); 255 void mps_secular_dnewton (mps_context * st, mps_polynomial * p, mps_approximation * root, cdpe_t corr); 256 void mps_secular_mnewton (mps_context * st, mps_polynomial * p, mps_approximation * root, mpc_t corr, long int wp); 257 258 /* Routines in secular-regeneartion.c */ 259 mps_boolean mps_secular_ga_regenerate_coefficients (mps_context * s); 260 261 /* Routines in secular.c */ 262 void mps_secular_deflate (mps_context * s, mps_secular_equation * sec); 263 264 void mps_secular_check_data (mps_context * s, char *which_case); 265 266 void mps_secular_restart (mps_context * s); 267 268 void mps_secular_switch_phase (mps_context * s, mps_phase phase); 269 270 long int mps_secular_raise_coefficient_precision (mps_context * s, mps_polynomial * p, long int wp); 271 272 void mps_secular_raise_precision (mps_context * s, int wp); 273 274 void mps_secular_raise_root_precision (mps_context * s, int wp); 275 276 /* Routines in secular-starting.c */ 277 void mps_secular_fstart (mps_context * s, mps_secular_equation * sec, mps_approximation ** approximations); 278 void mps_secular_dstart (mps_context * s, mps_secular_equation * sec, mps_approximation ** approximations); 279 void mps_secular_mstart (mps_context * s, mps_secular_equation * sec, mps_approximation ** approximations); 280 281 /* Routines in secular-iteration.c */ 282 int mps_secular_ga_fiterate (mps_context * s, int maxit, mps_boolean just_regenerated); 283 284 int mps_secular_ga_diterate (mps_context * s, int maxit, mps_boolean just_regenerated); 285 286 int mps_secular_ga_miterate (mps_context * s, int maxit, mps_boolean just_regenerated); 287 288 /* Routines in secular-ga.c */ 289 mps_boolean mps_secular_ga_check_stop (mps_context * s); 290 291 void mps_secular_ga_update_coefficients (mps_context * s); 292 293 /* Interface functions in secular.c */ 294 mps_secular_equation *mps_secular_equation_new (mps_context * s, 295 cplx_t * afpc, 296 cplx_t * bfpc, 297 unsigned long int n); 298 299 mps_secular_equation *mps_secular_equation_new_raw (mps_context * s, 300 unsigned long int n); 301 302 void mps_secular_equation_set_coefficient_f (mps_context *ctx, mps_secular_equation *p, int i, mpc_t a, mpc_t b); 303 void mps_secular_equation_set_coefficient_q (mps_context *ctx, mps_secular_equation *p, int i, mpq_t ar, mpq_t ai, mpq_t br, mpq_t bi); 304 305 void mps_secular_equation_free (mps_context * ctx, mps_polynomial * p); 306 307 void mps_secular_set_radii (mps_context * s); 308 309 mps_boolean mps_secular_poly_feval_with_error (mps_context * ctx, mps_polynomial * p, cplx_t x, cplx_t value, double * error); 310 311 mps_boolean mps_secular_poly_deval_with_error (mps_context * ctx, mps_polynomial * p, cdpe_t x, cdpe_t value, rdpe_t error); 312 313 mps_boolean mps_secular_poly_meval_with_error (mps_context * ctx, mps_polynomial * p, mpc_t x, mpc_t value, rdpe_t error); 314 315 void mps_secular_poly_fstart (mps_context * ctx, mps_polynomial * p, mps_approximation ** approximations); 316 317 void mps_secular_poly_dstart (mps_context * ctx, mps_polynomial * p, mps_approximation ** approximations); 318 319 void mps_secular_poly_mstart (mps_context * ctx, mps_polynomial * p, mps_approximation ** approximations); 320 321 mps_secular_equation * mps_secular_equation_read_from_stream (mps_context * ctx, mps_input_buffer * buffer, 322 mps_structure structure, mps_density density, 323 long int precision); 324 325 326 MPS_END_DECLS 327 328 #endif /* SECULAR_H */ 329