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 * Dario Andrea Bini <bini@dm.unipi.it>
9 * Giuseppe Fiorentino <fiorent@dm.unipi.it>
10 * Leonardo Robol <leonardo.robol@unipi.it>
11 */
12
13 #include <mps/mps.h>
14
15 /**
16 * @brief Utility routine that dumps the DPE coefficients of the secular equation
17 * passed as second argument. It is only used for debugging purpose.
18 *
19 * @param s The mps_context of the computation
20 * @param sec The secular equation whose DPE coefficients will be dumped
21 */
22 void
mps_secular_dump(mps_context * s,mps_secular_equation * sec)23 mps_secular_dump (mps_context * s, mps_secular_equation * sec)
24 {
25 int i;
26
27 MPS_DEBUG (s, "Dumping secular equation:");
28
29 switch (s->lastphase)
30 {
31 case float_phase:
32 for (i = 0; i < s->n; ++i)
33 {
34 MPS_DEBUG_CPLX (s, sec->afpc[i], "sec->afpc[%d]", i);
35 MPS_DEBUG_CPLX (s, sec->bfpc[i], "sec->bfpc[%d]", i);
36 }
37 break;
38
39 case dpe_phase:
40 for (i = 0; i < MPS_POLYNOMIAL (sec)->degree; i++)
41 {
42 MPS_DEBUG_CDPE (s, sec->adpc[i], "sec->adpc[%d]", i);
43 MPS_DEBUG_CDPE (s, sec->bdpc[i], "sec->bdpc[%d]", i);
44 }
45 break;
46
47 case mp_phase:
48 for (i = 0; i < s->n; ++i)
49 {
50 MPS_DEBUG_MPC (s, 20, sec->ampc[i], "sec->ampc[%d]", i);
51 MPS_DEBUG_MPC (s, 20, sec->bmpc[i], "sec->bmpc[%d]", i);
52 }
53 break;
54
55 default:
56 break;
57 }
58 }
59
60 void
mps_secular_restart(mps_context * s)61 mps_secular_restart (mps_context * s)
62 {
63 MPS_DEBUG_THIS_CALL (s);
64
65 int i;
66
67 switch (s->lastphase)
68 {
69 case float_phase:
70 for (i = 0; i < s->n; i++)
71 mpc_set_cplx (s->root[i]->mvalue, s->root[i]->fvalue);
72 break;
73
74 case dpe_phase:
75 for (i = 0; i < s->n; i++)
76 mpc_set_cdpe (s->root[i]->mvalue, s->root[i]->dvalue);
77 break;
78
79 default:
80 break;
81 }
82
83 mps_mrestart (s);
84
85 for (i = 0; i < s->n; i++)
86 {
87 mpc_get_cplx (s->root[i]->fvalue, s->root[i]->mvalue);
88 mpc_get_cdpe (s->root[i]->dvalue, s->root[i]->mvalue);
89 }
90 }
91
92 /**
93 * @brief Deflate a secular equation lowering the degree of the
94 * polynomial that represent it, if that is possible.
95 *
96 * Please note the <code>s->n</code> and <code>s->deg</code>
97 * will not be touched by this routine, so you should check
98 * that they are set according to <code>MPS_POLYNOMIAL (sec)->degree</code> if deflation
99 * takes place.
100 *
101 * @see <code>mps_context_set_degree ()</code>
102 *
103 * @param s The mps_context of the computation
104 * @param sec The secular equation that will be deflated.
105 */
106 void
mps_secular_deflate(mps_context * s,mps_secular_equation * sec)107 mps_secular_deflate (mps_context * s, mps_secular_equation * sec)
108 {
109 int i, j, k;
110
111 /* If the input is floating point check on the
112 * DPE input */
113 if (MPS_STRUCTURE_IS_FP (MPS_POLYNOMIAL (sec)->structure))
114 {
115 /* Do not deflate in floating point, since it is not working
116 * correctly right now */
117 MPS_DEBUG_WITH_INFO (s, "Floating point deflation still has some rough edges, so it's disabled");
118 return;
119 }
120
121 for (i = 0; i < MPS_POLYNOMIAL (sec)->degree; i++)
122 {
123 for (j = i + 1; j < MPS_POLYNOMIAL (sec)->degree; j++)
124 {
125 /* Otherwise, in the case of rational or integer input
126 * (that are handled in the same way) use initial_*mqpc
127 * values */
128 if (MPS_STRUCTURE_IS_INTEGER (MPS_POLYNOMIAL (sec)->structure) ||
129 MPS_STRUCTURE_IS_RATIONAL (MPS_POLYNOMIAL (sec)->structure))
130 {
131 if (mpq_equal (sec->initial_bmpqrc[i], sec->initial_bmpqrc[j])
132 && mpq_equal (sec->initial_bmpqic[i],
133 sec->initial_bmpqic[j]))
134 {
135 MPS_DEBUG_WITH_INFO (s,
136 "Coefficients b[%d] and b[%d] are equal, deflating",
137 i, j);
138 mpq_add (sec->initial_ampqrc[i], sec->initial_ampqrc[i],
139 sec->initial_ampqrc[j]);
140 mpq_add (sec->initial_ampqic[i], sec->initial_ampqic[i],
141 sec->initial_ampqic[j]);
142
143
144 /* Copy other coefficients back of one position */
145 for (k = j; k < MPS_POLYNOMIAL (sec)->degree - 1; k++)
146 {
147 mpq_set (sec->initial_ampqrc[k],
148 sec->initial_ampqrc[k + 1]);
149 mpq_set (sec->initial_ampqic[k],
150 sec->initial_ampqic[k + 1]);
151 }
152
153 MPS_POLYNOMIAL (sec)->degree--;
154 j--;
155 }
156 }
157 }
158 }
159
160 /* If the input was rational or integer, we need to reset the dpe coefficients
161 * according to it */
162 if (MPS_STRUCTURE_IS_INTEGER (MPS_POLYNOMIAL (sec)->structure) ||
163 MPS_STRUCTURE_IS_RATIONAL (MPS_POLYNOMIAL (sec)->structure))
164 {
165 mpf_t ftmp;
166 mpf_init (ftmp);
167
168 /* Set DPE coefficients */
169 for (i = 0; i < MPS_POLYNOMIAL (sec)->degree; i++)
170 {
171 mpf_set_q (ftmp, sec->initial_ampqrc[i]);
172 mpf_get_rdpe (cdpe_Re (sec->adpc[i]), ftmp);
173
174 mpf_set_q (ftmp, sec->initial_ampqic[i]);
175 mpf_get_rdpe (cdpe_Im (sec->adpc[i]), ftmp);
176
177 mpf_set_q (ftmp, sec->initial_bmpqrc[i]);
178 mpf_get_rdpe (cdpe_Re (sec->bdpc[i]), ftmp);
179
180 mpf_set_q (ftmp, sec->initial_bmpqic[i]);
181 mpf_get_rdpe (cdpe_Im (sec->bdpc[i]), ftmp);
182 }
183
184 mpf_clear (ftmp);
185 }
186
187 /* If the input was floating point update the coefficients using initial_*mpc
188 * values */
189 if (MPS_STRUCTURE_IS_FP (MPS_POLYNOMIAL (sec)->structure))
190 {
191 for (i = 0; i < MPS_POLYNOMIAL (sec)->degree; i++)
192 {
193 mpc_get_cdpe (sec->adpc[i], sec->ampc[i]);
194 mpc_get_cdpe (sec->bdpc[i], sec->bmpc[i]);
195 }
196 }
197
198 MPS_DEBUG (s, "Secular equation deflated to degree %d", MPS_POLYNOMIAL (sec)->degree);
199 }
200
201 /**
202 * @brief Raw version of mps_secular_equation_new that only
203 * allocate space for the coefficients but relies on the user
204 * to fill their values.
205 *
206 * @param s The mps_context of the computation.
207 * @param n The degree of the new secular equation to be created.
208 */
209 mps_secular_equation *
mps_secular_equation_new_raw(mps_context * s,unsigned long int n)210 mps_secular_equation_new_raw (mps_context * s, unsigned long int n)
211 {
212 int i;
213 mps_secular_equation *sec =
214 (mps_secular_equation*)mps_malloc (sizeof(mps_secular_equation));
215
216 mps_polynomial_init (s, MPS_POLYNOMIAL (sec));
217
218 MPS_POLYNOMIAL (sec)->type_name = "mps_secular_equation";
219
220 /* Hook up the overloaded methods for secular equations */
221 mps_polynomial * p = MPS_POLYNOMIAL (sec);
222 p->feval = mps_secular_poly_feval_with_error;
223 p->deval = mps_secular_poly_deval_with_error;
224 p->meval = mps_secular_poly_meval_with_error;
225 p->fstart = mps_secular_poly_fstart;
226 p->dstart = mps_secular_poly_dstart;
227 p->mstart = mps_secular_poly_mstart;
228 p->free = mps_secular_equation_free;
229 p->raise_data = mps_secular_raise_coefficient_precision;
230 p->fnewton = mps_secular_fnewton;
231 p->dnewton = mps_secular_dnewton;
232 p->mnewton = mps_secular_mnewton;
233 p->prec = 0;
234
235 /* Allocate floating point coefficients */
236 sec->afpc = cplx_valloc (n);
237 sec->bfpc = cplx_valloc (n);
238
239 /* Allocate complex dpe coefficients of the secular equation */
240 sec->adpc = cdpe_valloc (n);
241 sec->bdpc = cdpe_valloc (n);
242
243 /* Allocate multiprecision complex coefficients of the secular equation */
244 /* Prepare the double buffer */
245 struct mps_secular_equation_double_buffer *db = &sec->db;
246 db->ampc1 = mpc_valloc (n);
247 db->ampc2 = mpc_valloc (n);
248 db->bmpc1 = mpc_valloc (n);
249 db->bmpc2 = mpc_valloc (n);
250
251 sec->ampc = db->ampc1;
252 sec->bmpc = db->bmpc1;
253
254 db->active = 1;
255
256 sec->initial_ampc = mpc_valloc (n);
257 sec->initial_bmpc = mpc_valloc (n);
258 sec->initial_ampqrc = mpq_valloc (n);
259 sec->initial_bmpqrc = mpq_valloc (n);
260 sec->initial_ampqic = mpq_valloc (n);
261 sec->initial_bmpqic = mpq_valloc (n);
262
263 /* Allocate space for the moduli of the coefficients */
264 sec->aadpc = rdpe_valloc (n);
265 sec->abdpc = rdpe_valloc (n);
266 sec->aafpc = double_valloc (n);
267 sec->abfpc = double_valloc (n);
268
269 /* Init multiprecision arrays */
270 mpc_vinit2 (sec->db.ampc1, n, s->mpwp);
271 mpc_vinit2 (sec->db.bmpc1, n, s->mpwp);
272 mpc_vinit2 (sec->db.ampc2, n, s->mpwp);
273 mpc_vinit2 (sec->db.bmpc2, n, s->mpwp);
274 mpc_vinit2 (sec->initial_ampc, n, s->mpwp);
275 mpc_vinit2 (sec->initial_bmpc, n, s->mpwp);
276 mpq_vinit (sec->initial_ampqrc, n);
277 mpq_vinit (sec->initial_bmpqrc, n);
278 mpq_vinit (sec->initial_ampqic, n);
279 mpq_vinit (sec->initial_bmpqic, n);
280
281 MPS_POLYNOMIAL (sec)->degree = n;
282
283 /* Set up the mutexes for thread safety */
284 sec->ampc_mutex = mps_newv (pthread_mutex_t, MPS_POLYNOMIAL (sec)->degree);
285 sec->bmpc_mutex = mps_newv (pthread_mutex_t, MPS_POLYNOMIAL (sec)->degree);
286
287 for (i = 0; i < n; i++)
288 {
289 pthread_mutex_init (&sec->ampc_mutex[i], NULL);
290 pthread_mutex_init (&sec->bmpc_mutex[i], NULL);
291 }
292
293 pthread_mutex_init (&sec->precision_mutex, NULL);
294
295 /* Set the most general possible structure */
296 MPS_POLYNOMIAL (sec)->structure = MPS_STRUCTURE_UNKNOWN;
297
298 return sec;
299 }
300
301 /**
302 * @brief Create a new secular equation struct
303 *
304 * @param s The mps_context of the computation.
305 * @param afpc The floating point complex numerator coefficients.
306 * @param bfpc The floating point complex denominator coefficients.
307 * @param n The degree of the secular equation.
308 */
309 mps_secular_equation *
mps_secular_equation_new(mps_context * s,cplx_t * afpc,cplx_t * bfpc,unsigned long int n)310 mps_secular_equation_new (mps_context * s, cplx_t * afpc, cplx_t * bfpc,
311 unsigned long int n)
312 {
313 int i;
314
315 /* Allocate the space for the new struct */
316 mps_secular_equation *sec = mps_secular_equation_new_raw (s, n);
317
318 /* Copy the complex coefficients passed as argument */
319 for (i = 0; i < n; i++)
320 {
321 /* a_i coefficients */
322 cplx_set (sec->afpc[i], afpc[i]);
323
324 /* b_i coefficients */
325 cplx_set (sec->bfpc[i], bfpc[i]);
326 }
327
328 for (i = 0; i < MPS_POLYNOMIAL (sec)->degree; i++)
329 {
330 cdpe_init (sec->adpc[i]);
331 cdpe_set_x (sec->adpc[i], sec->afpc[i]);
332
333 mpc_set_cplx (sec->ampc[i], sec->afpc[i]);
334
335 cdpe_init (sec->bdpc[i]);
336 cdpe_set_x (sec->bdpc[i], sec->bfpc[i]);
337
338 mpc_set_cplx (sec->bmpc[i], sec->bfpc[i]);
339 }
340
341 MPS_POLYNOMIAL (sec)->structure = MPS_STRUCTURE_COMPLEX_FP;
342
343 return sec;
344 }
345
346 /**
347 * @brief Sets a coefficient of a secular equation.
348 *
349 * @param ctx The current context.
350 * @param p The secular equation whose coefficient should be set.
351 * @param i The index of the coefficient to set.
352 * @param a The numerator of the coefficient.
353 * @param b The zero of the denominator.
354 *
355 * The new coefficient will have the form a / (x - b).
356 */
357 void
mps_secular_equation_set_coefficient_f(mps_context * ctx,mps_secular_equation * p,int i,mpc_t a,mpc_t b)358 mps_secular_equation_set_coefficient_f (mps_context *ctx, mps_secular_equation *p, int i, mpc_t a, mpc_t b)
359 {
360 /* Updating data_type information */
361 if (MPS_POLYNOMIAL (p)->structure == MPS_STRUCTURE_UNKNOWN)
362 MPS_POLYNOMIAL (p)->structure = MPS_STRUCTURE_COMPLEX_FP;
363
364 long wp = mpc_get_prec (a);
365
366 if (mpc_get_prec (b) > wp)
367 wp = mpc_get_prec (b);
368
369 mpc_set_prec(p->initial_ampc[i], wp);
370 mpc_set_prec(p->initial_bmpc[i], wp);
371
372 mpc_set (p->initial_ampc[i], a);
373 mpc_set (p->initial_bmpc[i], b);
374
375 mps_secular_raise_coefficient_precision (ctx, MPS_POLYNOMIAL (p), wp);
376
377 mpc_set (p->ampc[i], a);
378 mpc_set (p->bmpc[i], b);
379
380 mpc_get_cplx (p->afpc[i], p->ampc[i]);
381 mpc_get_cplx (p->bfpc[i], p->bmpc[i]);
382
383 mpc_get_cdpe (p->adpc[i], p->ampc[i]);
384 mpc_get_cdpe (p->bdpc[i], p->bmpc[i]);
385
386 mpc_rmod (p->aadpc[i], p->ampc[i]);
387 mpc_rmod (p->abdpc[i], p->bmpc[i]);
388
389 p->aafpc[i] = rdpe_get_d (p->aadpc[i]);
390 p->abfpc[i] = rdpe_get_d (p->abdpc[i]);
391 }
392
393 /**
394 * @brief Sets a coefficient of a secular equation.
395 *
396 * @param ctx The current context.
397 * @param p The secular equation whose coefficient should be set.
398 * @param i The index of the coefficient to set.
399 * @param a The numerator of the coefficient.
400 * @param b The zero of the denominator.
401 *
402 * The new coefficient will have the form a / (x - b).
403 */
404 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)405 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)
406 {
407 /* Updating data_type information */
408 if (MPS_POLYNOMIAL (p)->structure == MPS_STRUCTURE_UNKNOWN)
409 MPS_POLYNOMIAL (p)->structure = MPS_STRUCTURE_COMPLEX_RATIONAL;
410
411 mpq_set (p->initial_ampqrc[i], ar);
412 mpq_set (p->initial_ampqic[i], ai);
413 mpq_set (p->initial_bmpqrc[i], br);
414 mpq_set (p->initial_bmpqic[i], bi);
415
416 mps_secular_raise_coefficient_precision (ctx, MPS_POLYNOMIAL (p), mpc_get_prec(p->ampc[0]));
417
418 mpc_set_q (p->ampc[i], p->initial_ampqrc[i], p->initial_ampqic[i]);
419 mpc_set_q (p->bmpc[i], p->initial_bmpqrc[i], p->initial_bmpqic[i]);
420
421 mpc_get_cplx (p->afpc[i], p->ampc[i]);
422 mpc_get_cplx (p->bfpc[i], p->bmpc[i]);
423
424 mpc_get_cdpe (p->adpc[i], p->ampc[i]);
425 mpc_get_cdpe (p->bdpc[i], p->bmpc[i]);
426
427 mpc_rmod (p->aadpc[i], p->ampc[i]);
428 mpc_rmod (p->abdpc[i], p->bmpc[i]);
429
430 p->aafpc[i] = rdpe_get_d (p->aadpc[i]);
431 p->abfpc[i] = rdpe_get_d (p->abdpc[i]);
432 }
433
434 /**
435 * @brief Free a secular equation and the data in it.
436 *
437 * @param ctx The current context.
438 * @param p The secular equation casted to a mps_polynomial
439 */
440 void
mps_secular_equation_free(mps_context * ctx,mps_polynomial * p)441 mps_secular_equation_free (mps_context *ctx, mps_polynomial * p)
442 {
443 mps_secular_equation *s = MPS_SECULAR_EQUATION (p);
444
445 /* Free internal data */
446 cplx_vfree (s->afpc);
447 cplx_vfree (s->bfpc);
448
449 cdpe_vfree (s->adpc);
450 cdpe_vfree (s->bdpc);
451
452 mpc_vclear (s->db.ampc1, MPS_POLYNOMIAL (s)->degree);
453 mpc_vclear (s->db.bmpc1, MPS_POLYNOMIAL (s)->degree);
454 mpc_vclear (s->db.ampc2, MPS_POLYNOMIAL (s)->degree);
455 mpc_vclear (s->db.bmpc2, MPS_POLYNOMIAL (s)->degree);
456
457 mpc_vfree (s->db.ampc1);
458 mpc_vfree (s->db.ampc2);
459 mpc_vfree (s->db.bmpc1);
460 mpc_vfree (s->db.bmpc2);
461
462 rdpe_vfree (s->aadpc);
463 rdpe_vfree (s->abdpc);
464 double_vfree (s->aafpc);
465 double_vfree (s->abfpc);
466
467 /* And old coefficients */
468 mpc_vclear (s->initial_ampc, MPS_POLYNOMIAL (s)->degree);
469 mpc_vclear (s->initial_bmpc, MPS_POLYNOMIAL (s)->degree);
470 mpq_vclear (s->initial_ampqrc, MPS_POLYNOMIAL (s)->degree);
471 mpq_vclear (s->initial_bmpqrc, MPS_POLYNOMIAL (s)->degree);
472 mpq_vclear (s->initial_ampqic, MPS_POLYNOMIAL (s)->degree);
473 mpq_vclear (s->initial_bmpqic, MPS_POLYNOMIAL (s)->degree);
474 mpc_vfree (s->initial_ampc);
475 mpc_vfree (s->initial_bmpc);
476 mpq_vfree (s->initial_ampqrc);
477 mpq_vfree (s->initial_bmpqrc);
478 mpq_vfree (s->initial_ampqic);
479 mpq_vfree (s->initial_bmpqic);
480
481 /* Mutexes */
482 free (s->ampc_mutex);
483 free (s->bmpc_mutex);
484
485 /* ...and then release it */
486 free (s);
487 }
488
489
490 /**
491 * @brief Evaluate secular equation in the point x.
492 *
493 * The evalutation will be done in floating point and this routine
494 * is used only for debugging purpose.
495 *
496 * @param s The mps_context of the computation.
497 * @param x The point in which the secular equation will be evaluated.
498 * @param sec_ev The result of the evalutation (output variable).
499 */
500 void
mps_secular_evaluate(mps_context * s,cplx_t x,cplx_t sec_ev)501 mps_secular_evaluate (mps_context * s, cplx_t x, cplx_t sec_ev)
502 {
503 cplx_t ctmp;
504 int i;
505 mps_secular_equation *sec = (mps_secular_equation*)s->secular_equation;
506
507 cplx_set (sec_ev, cplx_zero);
508
509 for (i = 0; i < s->n; i++)
510 {
511 /* Compute 1 / (x - b_i) */
512 cplx_sub (ctmp, x, sec->bfpc[i]);
513 cplx_inv_eq (ctmp);
514
515 /* Compute a_i / (x - b_i) */
516 cplx_mul_eq (ctmp, sec->afpc[i]);
517
518 /* Sum to the secular eqation */
519 cplx_add_eq (sec_ev, ctmp);
520 }
521
522 cplx_sub_eq (sec_ev, cplx_one);
523 }
524
525 /**
526 * @brief Secular version of <code>mps_check_data ()</code> that
527 * does nothing except to set the starting case according to
528 * <code>sec->starting_case</code>.
529 *
530 * @param s The current mps_context.
531 * @param which_case Pointer to a char that will be set to 'f' or 'd' depending
532 * on the chosen start phase.
533 */
534 void
mps_secular_check_data(mps_context * s,char * which_case)535 mps_secular_check_data (mps_context * s, char *which_case)
536 {
537 /* While we can't found a good criterion to check
538 * the possibility to start in pure floating point we
539 * use the DPE version. */
540 *which_case = (s->input_config->starting_phase == float_phase) ? 'f' : 'd';
541 }
542
543 /**
544 * @brief Raise precision of the coefficient of the secular equation
545 * (not the roots and neither the precison of the system) to <code>wp</code>.
546 *
547 * @param s The mps_context of the computation.
548 * @param p The secular equation casted to a mps_polynomial.
549 * @param wp The bits of precision to which the coefficients will be set.
550 *
551 * @see <code>mps_secular_raise_root_precision ()</code>
552 * @see <code>mps_secular_raise_precision ()</code>
553 */
554 long int
mps_secular_raise_coefficient_precision(mps_context * s,mps_polynomial * p,long int wp)555 mps_secular_raise_coefficient_precision (mps_context * s, mps_polynomial * p, long int wp)
556 {
557 MPS_DEBUG_THIS_CALL (s);
558
559 int i;
560 mps_secular_equation *sec = MPS_SECULAR_EQUATION (p);
561
562 mpc_t * raising_ampc;
563 mpc_t * raising_bmpc;
564
565 pthread_mutex_lock (&sec->precision_mutex);
566
567 if (wp < mpc_get_prec (sec->ampc[0]))
568 {
569 pthread_mutex_unlock (&sec->precision_mutex);
570 return mpc_get_prec (sec->ampc[0]);
571 }
572
573 if (sec->db.active == 1)
574 {
575 raising_ampc = sec->db.ampc2;
576 raising_bmpc = sec->db.bmpc2;
577 }
578 else
579 {
580 raising_ampc = sec->db.ampc1;
581 raising_bmpc = sec->db.bmpc1;
582 }
583
584 for (i = 0; i < p->degree; i++)
585 {
586 mpc_set_prec (raising_ampc[i], wp);
587 if (!MPS_STRUCTURE_IS_FP (p->structure))
588 {
589 mpf_set_q (mpc_Re (raising_ampc[i]), sec->initial_ampqrc[i]);
590 mpf_set_q (mpc_Im (raising_ampc[i]), sec->initial_ampqic[i]);
591 }
592 else
593 mpc_set (raising_ampc[i], sec->ampc[i]);
594
595 mpc_set_prec (raising_bmpc[i], wp);
596 if (!MPS_STRUCTURE_IS_FP (p->structure))
597 {
598 mpf_set_q (mpc_Re (raising_bmpc[i]), sec->initial_bmpqrc[i]);
599 mpf_set_q (mpc_Im (raising_bmpc[i]), sec->initial_bmpqic[i]);
600 }
601 else
602 mpc_set (raising_bmpc[i], sec->bmpc[i]);
603 }
604
605 sec->ampc = raising_ampc;
606 sec->bmpc = raising_bmpc;
607
608 sec->db.active = (sec->db.active % 2) + 1;
609
610 pthread_mutex_unlock (&sec->precision_mutex);
611
612 if (s->debug_level & MPS_DEBUG_MEMORY)
613 MPS_DEBUG_WITH_INFO (s, "Precision of the coefficients is now at %ld bits", wp);
614
615 return mpc_get_prec (sec->ampc[0]);
616 }
617
618 /**
619 * @brief Raise precision of the roots (not the coefficients nor the
620 * system) to <code>wp</code> bits.
621 *
622 * @param s The mps_context of the computation.
623 * @param wp The bits of precision to which the roots will be set.
624 *
625 * @see <code>mps_secular_raise_coefficient_precision ()</code>
626 * @see <code>mps_secular_raise_precision ()</code>
627 */
628 void
mps_secular_raise_root_precision(mps_context * s,int wp)629 mps_secular_raise_root_precision (mps_context * s, int wp)
630 {
631 MPS_DEBUG_THIS_CALL (s);
632 int i;
633
634 for (i = 0; i < s->n; i++)
635 {
636 mpc_set_prec (s->root[i]->mvalue, wp);
637 }
638 }
639
640 /**
641 * @brief Raise (or lower) the precision of the coefficients to
642 * <code>wp</code> bits. This will change precision of the coefficients
643 * via <code>mps_secular_raise_coefficient_precision ()</code>,
644 * of the roots via <code>mps_secular_raise_root_precision ()</code>
645 * and set <code>s->mpwp</code>.
646 *
647 * @param s The mps_context of the computation.
648 * @param wp The bits of precision to which all the computation will
649 * be brought.
650 */
651 void
mps_secular_raise_precision(mps_context * s,int wp)652 mps_secular_raise_precision (mps_context * s, int wp)
653 {
654 MPS_DEBUG_THIS_CALL (s);
655
656 int i;
657
658 mps_secular_raise_coefficient_precision (s, MPS_POLYNOMIAL (s->secular_equation), wp);
659 mps_secular_raise_root_precision (s, wp);
660
661 s->mpwp = wp;
662 rdpe_set_2dl (s->mp_epsilon, 1.0, -wp);
663
664 s->just_raised_precision = true;
665
666
667 for (i = 0; i < s->n; i++)
668 {
669 s->root[i]->approximated = false;
670 s->root[i]->again = true;
671 }
672 }
673
674 /**
675 * @brief Prepare data for the iteration in the new phase specified
676 * in the second parameter.
677 *
678 * Note that for now this function is only able to handle switch
679 * from floating point phases (i.e. float_phase or dpe_phase) to
680 * multiprecision, and not coming back.
681 *
682 * @param s The mps_context of the computation.
683 * @param phase The phase to switch the computation to.
684 */
685 void
mps_secular_switch_phase(mps_context * s,mps_phase phase)686 mps_secular_switch_phase (mps_context * s, mps_phase phase)
687 {
688 MPS_DEBUG_THIS_CALL (s);
689
690 s->just_raised_precision = true;
691
692 int i = 0;
693 mps_secular_equation *sec = (mps_secular_equation*)s->secular_equation;
694
695 if (phase == mp_phase)
696 {
697 /* Debug the approximations that we have now before going
698 * to the multiprecision phase */
699 if (s->debug_level & MPS_DEBUG_APPROXIMATIONS)
700 {
701 MPS_DEBUG (s, "Dumping current approximations before starting MP");
702 mps_dump (s);
703 }
704
705 mps_secular_raise_precision (s, MPS_SECULAR_STARTING_MP_PRECISION);
706 switch (s->lastphase)
707 {
708 case float_phase:
709 /* Copy the approximated roots and the
710 * secular equation coefficients */
711 for (i = 0; i < s->n; i++)
712 {
713 mpc_set_cplx (s->root[i]->mvalue, s->root[i]->fvalue);
714 mpc_set_cplx (sec->ampc[i], sec->afpc[i]);
715 mpc_set_cplx (sec->bmpc[i], sec->bfpc[i]);
716 rdpe_set_d (s->root[i]->drad, s->root[i]->frad);
717 }
718 break;
719
720 case dpe_phase:
721 /* Copy the coefficients and the approximated
722 * roots into the multiprecision values */
723 for (i = 0; i < s->n; i++)
724 {
725 mpc_set_cdpe (s->root[i]->mvalue, s->root[i]->dvalue);
726 mpc_set_cdpe (sec->ampc[i], sec->adpc[i]);
727 mpc_set_cdpe (sec->bmpc[i], sec->bdpc[i]);
728 }
729
730 default:
731 break;
732 }
733
734 /* Set lastphase to mp_phase */
735 s->lastphase = mp_phase;
736
737 /* Set epsilon */
738 rdpe_set_2dl (s->mp_epsilon, 1.0, -s->mpwp + 1);
739 }
740 else
741 {
742 fprintf (stderr, "mps_secular_switch_phase is only able to manage\n"
743 "switches from float_phase or dpe_phase to mp_phase. Aborting.");
744 exit (EXIT_FAILURE);
745 }
746 }
747
748 /**
749 * @brief Update radii of the roots according to the coefficients
750 * of the secular equation in this moment, if they are better of
751 * the radii present now.
752 *
753 * @param s The mps_context of the computation.
754 */
755 void
mps_secular_set_radii(mps_context * s)756 mps_secular_set_radii (mps_context * s)
757 {
758 MPS_DEBUG_THIS_CALL (s);
759
760 int i;
761 mps_secular_equation *sec = (mps_secular_equation*)s->secular_equation;
762
763 /* DPE and multiprecision implementation */
764 rdpe_t rad, rad_eps, rtmp;
765 cdpe_t ctmp;
766 rdpe_t * drad = rdpe_valloc (s->n);
767
768 mpc_t mtmp;
769 mpc_init2 (mtmp, mps_context_get_data_prec_max (s));
770
771 if (s->lastphase == mp_phase)
772 rdpe_set (rad_eps, s->mp_epsilon);
773 else
774 rdpe_set_d (rad_eps, DBL_EPSILON);
775
776 rdpe_mul_eq_d (rad_eps, s->n * 4);
777 rdpe_add_eq (rad_eps, rdpe_one);
778
779 /* Check if the Gerschgorin's radii are more convenient */
780 for (i = 0; i < s->n; i++)
781 {
782 mpc_get_cdpe (ctmp, sec->ampc[i]);
783 cdpe_mod (rad, ctmp);
784
785 rdpe_mul_eq (rad, rad_eps);
786
787 rdpe_mul_eq_d (rad, (double)s->n);
788
789 rdpe_set (drad[i], rad);
790
791 mpc_rmod (rtmp, s->root[i]->mvalue);
792 if (s->lastphase == mp_phase)
793 rdpe_mul_eq (rtmp, s->mp_epsilon);
794 else
795 rdpe_mul_eq_d (rtmp, DBL_EPSILON);
796 rdpe_mul_eq_d (rtmp, 4.0);
797 rdpe_add_eq (drad[i], rtmp);
798 }
799
800 switch (s->lastphase)
801 {
802 case float_phase:
803 {
804 for (i = 0; i < s->n; i++)
805 {
806 rdpe_set_d (s->root[i]->drad, s->root[i]->frad);
807 mpc_set_d (s->root[i]->mvalue, cplx_Re (s->root[i]->fvalue),
808 cplx_Im (s->root[i]->fvalue));
809 }
810
811 mps_mcluster (s, drad, 2.0 * s->n);
812 mps_fmodify (s, false);
813
814 for (i = 0; i < s->n; i++)
815 {
816 s->root[i]->frad = rdpe_get_d (s->root[i]->drad);
817 if (s->root[i]->frad == 0.0)
818 s->root[i]->frad += cplx_mod (s->root[i]->fvalue) * DBL_EPSILON;
819 }
820 }
821 break;
822
823 case dpe_phase:
824 mps_mcluster (s, drad, 2.0 * s->n);
825 mps_dmodify (s, false);
826 break;
827
828 case mp_phase:
829 mps_mcluster (s, drad, 2.0 * s->n);
830 mps_mmodify (s, false);
831 break;
832
833 default:
834 break;
835 }
836
837 rdpe_vfree (drad);
838
839 mpc_clear (mtmp);
840 }
841
842 mps_boolean
mps_secular_poly_feval_with_error(mps_context * ctx,mps_polynomial * p,cplx_t x,cplx_t value,double * error)843 mps_secular_poly_feval_with_error (mps_context * ctx, mps_polynomial * p, cplx_t x, cplx_t value, double * error)
844 {
845 cplx_t ctmp;
846 int i;
847 mps_secular_equation *sec = MPS_SECULAR_EQUATION (p);
848
849 if (!mps_secular_feval_with_error (ctx, p, x, value, error))
850 return false;
851
852 if (!cplx_eq_zero (value))
853 *error /= cplx_mod (value);
854 else
855 *error = p->degree * DBL_EPSILON;
856
857 for (i = 0; i < p->degree; i++)
858 {
859 cplx_sub (ctmp, x, sec->bfpc[i]);
860 cplx_mul_eq (value, ctmp);
861 }
862
863 cplx_mul_eq_d (value, -1.0);
864 *error *= cplx_mod (value);
865 return true;
866 }
867
868 mps_boolean
mps_secular_poly_deval_with_error(mps_context * ctx,mps_polynomial * p,cdpe_t x,cdpe_t value,rdpe_t error)869 mps_secular_poly_deval_with_error (mps_context * ctx, mps_polynomial * p, cdpe_t x, cdpe_t value, rdpe_t error)
870 {
871 cdpe_t ctmp;
872 rdpe_t rtmp;
873 int i;
874 mps_secular_equation *sec = MPS_SECULAR_EQUATION (p);
875
876 if (!mps_secular_deval_with_error (ctx, p, x, value, error))
877 return false;
878
879 cdpe_mod (rtmp, value);
880 if (!rdpe_eq_zero (rtmp))
881 rdpe_div_eq (error, rtmp);
882 else
883 rdpe_set_d (error, DBL_EPSILON * p->degree);
884
885 for (i = 0; i < p->degree; i++)
886 {
887 cdpe_sub (ctmp, x, sec->bdpc[i]);
888 cdpe_mul_eq (value, ctmp);
889 }
890
891 cdpe_mul_eq_d (value, -1.0);
892
893 cdpe_mod (rtmp, value);
894 rdpe_mul_eq (error, rtmp);
895
896 return true;
897 }
898
899 mps_boolean
mps_secular_poly_meval_with_error(mps_context * ctx,mps_polynomial * p,mpc_t x,mpc_t value,rdpe_t error)900 mps_secular_poly_meval_with_error (mps_context * ctx, mps_polynomial * p, mpc_t x, mpc_t value, rdpe_t error)
901 {
902 mpc_t ctmp;
903 rdpe_t rtmp;
904 int i;
905 mps_secular_equation *sec = MPS_SECULAR_EQUATION (p);
906
907 if (!mps_secular_meval_with_error (ctx, p, x, value, error))
908 {
909 return false;
910 }
911
912 mpc_rmod (rtmp, value);
913 if (!rdpe_eq_zero (rtmp))
914 rdpe_div_eq (error, rtmp);
915 else
916 rdpe_set_d (error, DBL_EPSILON * p->degree);
917
918 mpc_init2 (ctmp, mpc_get_prec (x));
919
920 for (i = 0; i < p->degree; i++)
921 {
922 mpc_sub (ctmp, x, sec->bmpc[i]);
923 mpc_mul_eq (value, ctmp);
924 }
925
926 mpc_set_si (ctmp, -1, 0);
927 mpc_mul_eq (value, ctmp);
928
929 mpc_clear (ctmp);
930
931 mpc_rmod (rtmp, value);
932 rdpe_mul_eq (error, rtmp);
933
934 return true;
935 }
936
937 void
mps_secular_poly_fstart(mps_context * ctx,mps_polynomial * p,mps_approximation ** approximations)938 mps_secular_poly_fstart (mps_context * ctx, mps_polynomial * p, mps_approximation ** approximations)
939 {
940 mps_secular_fstart (ctx, MPS_SECULAR_EQUATION (p), approximations);
941 }
942
943 void
mps_secular_poly_dstart(mps_context * ctx,mps_polynomial * p,mps_approximation ** approximations)944 mps_secular_poly_dstart (mps_context * ctx, mps_polynomial * p, mps_approximation ** approximations)
945 {
946 mps_secular_dstart (ctx, MPS_SECULAR_EQUATION (p), approximations);
947 }
948
949 void
mps_secular_poly_mstart(mps_context * ctx,mps_polynomial * p,mps_approximation ** approximations)950 mps_secular_poly_mstart (mps_context * ctx, mps_polynomial * p, mps_approximation ** approximations)
951 {
952 mps_secular_mstart (ctx, MPS_SECULAR_EQUATION (p), approximations);
953 }
954