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
14 #include <mps/mps.h>
15 #include <math.h>
16 #include <string.h>
17
18 /* Routine called to perform the floating point newton iterations. */
19 static void *
__mps_secular_ga_fiterate_worker(void * data_ptr)20 __mps_secular_ga_fiterate_worker (void* data_ptr)
21 {
22 mps_thread_worker_data *data = (mps_thread_worker_data*)data_ptr;
23 mps_context *s = data->s;
24 int i;
25 cplx_t corr, abcorr;
26 double modcorr;
27 mps_thread_job job;
28
29 while (true && !s->exit_required)
30 {
31 job = mps_thread_job_queue_next (s, data->queue);
32 i = job.i;
33
34 if (job.iter == MPS_THREAD_JOB_EXCEP || *data->nzeros >= s->n)
35 goto cleanup;
36
37 pthread_mutex_lock (&data->roots_mutex[i]);
38
39 if (job.iter == MPS_THREAD_JOB_EXCEP || *data->nzeros >= s->n)
40 {
41 pthread_mutex_unlock (&data->roots_mutex[i]);
42 goto cleanup;
43 }
44
45 if (s->root[i]->again && !s->root[i]->approximated)
46 {
47 /* Increment the number of performed iterations */
48 #if defined(__GCC__)
49 __sync_add_and_fetch (data->it, 1);
50 #else
51 pthread_mutex_lock (data->gs_mutex);
52 (*data->it)++;
53 pthread_mutex_unlock (data->gs_mutex);
54 #endif
55 cdpe_set_x (s->root[i]->dvalue, s->root[i]->fvalue);
56
57 mps_secular_fnewton (s, MPS_POLYNOMIAL (s->secular_equation), s->root[i], corr);
58
59 if (s->root[i]->status == MPS_ROOT_STATUS_NOT_FLOAT)
60 {
61 *data->excep = true;
62 pthread_mutex_unlock (&data->roots_mutex[i]);
63 break;
64 }
65
66 /* Apply Aberth correction */
67 mps_faberth_wl (s, i, abcorr, data->aberth_mutex);
68
69 if (isnan (cplx_Re (abcorr)) || isnan (cplx_Im (abcorr)))
70 {
71 s->root[i]->again = false;
72 pthread_mutex_unlock (&data->roots_mutex[i]);
73 continue;
74 }
75
76 cplx_mul_eq (abcorr, corr);
77 cplx_sub (abcorr, cplx_one, abcorr);
78 cplx_div (abcorr, corr, abcorr);
79
80 if (cplx_check_fpe (abcorr))
81 {
82 s->root[i]->again = false;
83 pthread_mutex_unlock (&data->roots_mutex[i]);
84 continue;
85 }
86
87 if (!s->root[i]->again || s->root[i]->approximated)
88 {
89 if (s->debug_level & MPS_DEBUG_APPROXIMATIONS)
90 MPS_DEBUG (s, "Root %d again was set to false on iteration %d by thread %d", i, *data->it, data->thread);
91
92 #if defined(__GCC__)
93 __sync_add_and_fetch (data->nzeros, 1);
94 #else
95 pthread_mutex_lock (data->gs_mutex);
96 (*data->nzeros)++;
97 pthread_mutex_unlock (data->gs_mutex);
98 #endif
99 }
100 else
101 {
102 pthread_mutex_lock (&data->aberth_mutex[i]);
103 cplx_sub_eq (s->root[i]->fvalue, abcorr);
104 pthread_mutex_unlock (&data->aberth_mutex[i]);
105
106 /* Correct the radius */
107 modcorr = cplx_mod (abcorr);
108 s->root[i]->frad += modcorr;
109 }
110 }
111
112 pthread_mutex_unlock (&data->roots_mutex[i]);
113 }
114
115 cleanup:
116
117 return NULL;
118 }
119
120 /**
121 * @brief Routine that performs a block of iteration
122 * in floating point on the secular equation.
123 *
124 * @param s the pointer to the mps_context struct.
125 * @param maxit Maximum number of iteration to perform.
126 * @param just_regenerated true if this is the first iteration after a coefficient
127 * regeneration. If just_regenerated is true and the iteration packet is completed
128 * in less than 2 * (n - computed_roots) iterations that best_approx is set to true
129 * in s->secular_equation so a raise in the precision will be triggered.
130 * @return The number of approximated roots after the iteration.
131 */
132 MPS_PRIVATE int
mps_secular_ga_fiterate(mps_context * s,int maxit,mps_boolean just_regenerated)133 mps_secular_ga_fiterate (mps_context * s, int maxit, mps_boolean just_regenerated)
134 {
135 int computed_roots = 0;
136 int approximated_roots = 0;
137 int root_neighborhood_roots = 0;
138 int i;
139 int nit = 0;
140 int it_threshold = 0;
141 mps_boolean excep = false;
142
143 #ifndef DISABLE_DEBUG
144 clock_t *my_clock = mps_start_timer ();
145 #endif
146
147 s->operation = MPS_OPERATION_ABERTH_FP_ITERATIONS;
148
149 mps_thread_worker_data *data;
150 pthread_mutex_t *aberth_mutex =
151 (pthread_mutex_t*)mps_malloc (sizeof(pthread_mutex_t) * s->n);
152 pthread_mutex_t *roots_mutex =
153 (pthread_mutex_t*)mps_malloc (sizeof(pthread_mutex_t) * s->n);
154
155 pthread_mutex_t gs_mutex = PTHREAD_MUTEX_INITIALIZER;
156
157 for (i = 0; i < s->n; i++)
158 {
159 pthread_mutex_init (roots_mutex + i, NULL);
160 pthread_mutex_init (aberth_mutex + i, NULL);
161 }
162
163 data = mps_newv (mps_thread_worker_data, s->n_threads);
164
165 MPS_DEBUG_THIS_CALL (s);
166
167 /* Mark the approximated roots as ready for output */
168 for (i = 0; i < s->n; i++)
169 {
170 /* Set again to false if the root is already approximated. If a root is approximated but
171 * it has less digits than the current precision don't stop the iterations on that component. */
172 if (s->root[i]->status == MPS_ROOT_STATUS_ISOLATED ||
173 s->root[i]->status == MPS_ROOT_STATUS_APPROXIMATED)
174 {
175 if (s->debug_level & MPS_DEBUG_APPROXIMATIONS)
176 {
177 MPS_DEBUG_WITH_INFO (s, "Setting again[%d] to false since the root is ready for output (or isolated)", i);
178 }
179 s->root[i]->again = false;
180
181 if (s->root[i]->status == MPS_ROOT_STATUS_APPROXIMATED)
182 s->root[i]->approximated = true;
183 }
184
185 if (!s->root[i]->again || s->root[i]->approximated)
186 computed_roots++;
187 }
188
189 MPS_DEBUG_WITH_INFO (s, "%d roots %s already approximated at the start of the packet",
190 computed_roots,
191 (computed_roots == 1) ? "is" : "are");
192
193 it_threshold = s->n - computed_roots;
194
195 mps_thread_job_queue *queue = mps_thread_job_queue_new (s);
196
197 for (i = 0; i < s->n_threads; i++)
198 {
199 data[i].it = &nit;
200 data[i].nzeros = &computed_roots;
201 data[i].s = s;
202 data[i].thread = i;
203 data[i].n_threads = s->n_threads;
204 data[i].aberth_mutex = aberth_mutex;
205 data[i].roots_mutex = roots_mutex;
206 data[i].queue = queue;
207 data[i].gs_mutex = &gs_mutex;
208 data[i].excep = &excep;
209
210 mps_thread_pool_assign (s, s->pool,
211 __mps_secular_ga_fiterate_worker, data + i);
212 }
213
214 mps_thread_pool_wait (s, s->pool);
215
216 /* Check if the roots are improvable in floating point */
217 MPS_DEBUG_WITH_INFO (s, "Performed %d iterations with floating point arithmetic",
218 nit);
219
220 if (s->debug_level & MPS_DEBUG_APPROXIMATIONS)
221 mps_dump (s);
222
223 /* Check if we need to get higher precision for the roots */
224 s->best_approx = true;
225 for (i = 0; i < s->n; i++)
226 {
227 if (!s->root[i]->approximated)
228 s->best_approx = false;
229 if (s->root[i]->approximated)
230 approximated_roots++;
231 if (!s->root[i]->again)
232 root_neighborhood_roots++;
233 }
234
235 if (just_regenerated && (nit <= it_threshold))
236 s->best_approx = true;
237
238 MPS_DEBUG_WITH_INFO (s, "%d roots are approximated with the current precision", approximated_roots);
239 MPS_DEBUG_WITH_INFO (s, "%d roots are in the root neighborhood", root_neighborhood_roots);
240 MPS_DEBUG_WITH_INFO (s, "%d roots have reached a stop condition", computed_roots);
241
242 if (excep)
243 {
244 MPS_DEBUG_WITH_INFO (s, "Switching to DPE arithmetic since there are roots not representable in standard floating point");
245 for (i = 0; i < s->n; i++)
246 {
247 cdpe_set_x (s->root[i]->dvalue, s->root[i]->fvalue);
248 rdpe_set_d (s->root[i]->drad, s->root[i]->frad);
249 s->root[i]->status = MPS_ROOT_STATUS_CLUSTERED;
250 }
251 s->lastphase = dpe_phase;
252 }
253
254 /* Compute the inclusion radii with Gerschgorin so we can compute
255 * clusterizations for the roots. */
256 /* mps_fradii (s, fradii); */
257 /* mps_fcluster (s, fradii, 2.0 * s->n); */
258 /* mps_fmodify (s, false); */
259
260 /* These lines are used to debug the again vector, but are not useful
261 * at the moment being */
262 if (s->debug_level & MPS_DEBUG_APPROXIMATIONS)
263 {
264 __MPS_DEBUG (s, "Again vector = ");
265 for (i = 0; i < s->n; i++)
266 {
267 fprintf (s->logstr, "%d ", s->root[i]->again);
268 }
269 fprintf (s->logstr, "\n");
270 }
271
272 if (s->debug_level & MPS_DEBUG_APPROXIMATIONS)
273 mps_dump (s);
274
275 /* Count time taken */
276 #ifndef DISABLE_DEBUG
277 s->fp_iteration_time += mps_stop_timer (my_clock);
278 #endif
279
280 mps_thread_job_queue_free (queue);
281 free (data);
282 free (roots_mutex);
283 free (aberth_mutex);
284
285 /* Return the number of approximated roots */
286 return computed_roots;
287 }
288
289 /* Routine called to perform the floating point newton iterations with DPE */
290 static void *
__mps_secular_ga_diterate_worker(void * data_ptr)291 __mps_secular_ga_diterate_worker (void* data_ptr)
292 {
293 mps_thread_worker_data *data = (mps_thread_worker_data*)data_ptr;
294 mps_context *s = data->s;
295 int i;
296 cdpe_t corr, abcorr, droot;
297 rdpe_t modcorr;
298 mps_thread_job job;
299
300 while (true && !s->exit_required)
301 {
302 job = mps_thread_job_queue_next (s, data->queue);
303 i = job.i;
304
305 if (job.iter == MPS_THREAD_JOB_EXCEP)
306 {
307 return NULL;
308 }
309
310 pthread_mutex_lock (&data->roots_mutex[i]);
311
312 if (s->root[i]->again && !s->root[i]->approximated)
313 {
314 /* Lock this roots to make sure that we are the only one working on it */
315 cdpe_set (droot, s->root[i]->dvalue);
316
317 (*data->it)++;
318
319 mps_secular_dnewton (s, MPS_POLYNOMIAL (s->secular_equation), s->root[i], corr);
320
321 /* Apply Aberth correction */
322 mps_daberth_wl (s, i, abcorr, data->aberth_mutex);
323 cdpe_mul_eq (abcorr, corr);
324 cdpe_sub (abcorr, cdpe_one, abcorr);
325 cdpe_div (abcorr, corr, abcorr);
326
327 cdpe_sub_eq (droot, abcorr);
328
329 /* Correct the radius */
330 if (s->root[i]->again)
331 {
332 cdpe_mod (modcorr, abcorr);
333 rdpe_add_eq (s->root[i]->drad, modcorr);
334 }
335
336 if (!s->root[i]->again || s->root[i]->approximated)
337 {
338 if (s->debug_level & MPS_DEBUG_APPROXIMATIONS)
339 MPS_DEBUG (s, "Root %d again was set to false on iteration %d by thread %d", i, *data->it, data->thread);
340 (*data->nzeros)++;
341 }
342 else
343 cdpe_set (s->root[i]->dvalue, droot);
344 }
345
346 pthread_mutex_unlock (&data->roots_mutex[i]);
347 }
348
349 return NULL;
350 }
351
352 /**
353 * @brief Routine that performs a block of iteration
354 * in floating point on the secular equation using
355 * CDPE
356 *
357 * @param s the pointer to the mps_context struct.
358 * @param maxit Maximum number of iteration to perform.
359 * @return The number of approximated roots after the iteration.
360 * @param just_regenerated true if this is the first iteration after a coefficient
361 * regeneration. If just_regenerated is true and the iteration packet is completed
362 * in less than 2 * (n - computed_roots) iterations that best_approx is set to true
363 * in s->secular_equation so a raise in the precision will be triggered.
364 */
365 MPS_PRIVATE int
mps_secular_ga_diterate(mps_context * s,int maxit,mps_boolean just_regenerated)366 mps_secular_ga_diterate (mps_context * s, int maxit, mps_boolean just_regenerated)
367 {
368 int computed_roots = 0;
369 int root_neighborhood_roots = 0;
370 int approximated_roots = 0;
371 int i;
372 int nit = 0;
373 int it_threshold = 0;
374
375 s->operation = MPS_OPERATION_ABERTH_DPE_ITERATIONS;
376
377 #ifndef DISABLE_DEBUG
378 clock_t *my_clock = mps_start_timer ();
379 #endif
380
381 mps_thread_worker_data *data;
382 pthread_mutex_t *aberth_mutex =
383 (pthread_mutex_t*)mps_malloc (sizeof(pthread_mutex_t) * s->n);
384 pthread_mutex_t *roots_mutex =
385 (pthread_mutex_t*)mps_malloc (sizeof(pthread_mutex_t) * s->n);
386
387 for (i = 0; i < s->n; i++)
388 {
389 pthread_mutex_init (roots_mutex + i, NULL);
390 pthread_mutex_init (aberth_mutex + i, NULL);
391 }
392
393 data = mps_newv (mps_thread_worker_data, s->n_threads);
394
395 MPS_DEBUG_THIS_CALL (s);
396
397 s->best_approx = false;
398
399 /* Mark the approximated roots as ready for output */
400 for (i = 0; i < s->n; i++)
401 {
402 /* Set again to false if the root is already approximated. If a root is approximated but
403 * it has less digits than the current precision don't stop the iterations on that component. */
404 if (s->root[i]->status == MPS_ROOT_STATUS_ISOLATED ||
405 s->root[i]->status == MPS_ROOT_STATUS_APPROXIMATED)
406 {
407 if (s->debug_level & MPS_DEBUG_APPROXIMATIONS)
408 {
409 MPS_DEBUG_WITH_INFO (s, "Setting again[%d] to false since the root is ready for output (or isolated)", i);
410 }
411 s->root[i]->again = false;
412 s->root[i]->approximated = true;
413 }
414
415 if (!s->root[i]->again || s->root[i]->approximated)
416 computed_roots++;
417 }
418
419 it_threshold = (s->n - computed_roots);
420
421 MPS_DEBUG_WITH_INFO (s, "%d roots %s already approximated at the start of the packet",
422 computed_roots,
423 (computed_roots == 1) ? "is" : "are");
424
425 mps_thread_job_queue *queue = mps_thread_job_queue_new (s);
426
427 for (i = 0; i < s->n_threads; i++)
428 {
429 data[i].it = &nit;
430 data[i].nzeros = &computed_roots;
431 data[i].s = s;
432 data[i].thread = i;
433 data[i].n_threads = s->n_threads;
434 data[i].aberth_mutex = aberth_mutex;
435 data[i].roots_mutex = roots_mutex;
436 data[i].queue = queue;
437
438 mps_thread_pool_assign (s, s->pool, __mps_secular_ga_diterate_worker, data + i);
439 }
440
441 mps_thread_pool_wait (s, s->pool);
442
443 /* Check if the roots are improvable in floating point */
444 MPS_DEBUG_WITH_INFO (s, "Performed %d iterations with CDPE arithmetic",
445 nit);
446
447 if (s->debug_level & MPS_DEBUG_APPROXIMATIONS)
448 mps_dump (s);
449
450 /* Check if we need to get higher precision for the roots */
451 s->best_approx = true;
452 for (i = 0; i < s->n; i++)
453 {
454 if (!s->root[i]->approximated)
455 s->best_approx = false;
456 if (s->root[i]->approximated)
457 approximated_roots++;
458 if (!s->root[i]->again)
459 root_neighborhood_roots++;
460 }
461
462 if (just_regenerated && (nit <= it_threshold))
463 s->best_approx = true;
464
465 MPS_DEBUG_WITH_INFO (s, "%d roots are approximated with the current precision", approximated_roots);
466 MPS_DEBUG_WITH_INFO (s, "%d roots are in the root neighborhood", root_neighborhood_roots);
467 MPS_DEBUG_WITH_INFO (s, "%d roots have reached a stop condition", computed_roots);
468
469 /* Clock the routine */
470 #ifndef DISABLE_DEBUG
471 s->dpe_iteration_time += mps_stop_timer (my_clock);
472 #endif
473
474 mps_thread_job_queue_free (queue);
475 free (aberth_mutex);
476 free (roots_mutex);
477
478 free (data);
479
480 /* Return the number of approximated roots */
481 return computed_roots;
482 }
483
484 /* Routine called to perform the floating point newton iterations with MP */
485 static void *
__mps_secular_ga_miterate_worker(void * data_ptr)486 __mps_secular_ga_miterate_worker (void* data_ptr)
487 {
488 mps_thread_worker_data *data = (mps_thread_worker_data*)data_ptr;
489 mps_context *s = data->s;
490 /* mps_secular_equation *sec = s->secular_equation; */
491 int i;
492 mpc_t corr, abcorr;
493 mpc_t mroot;
494 rdpe_t modcorr;
495 mps_thread_job job;
496
497 mps_cluster * cluster = NULL;
498
499 mpc_init2 (corr, s->mpwp);
500 mpc_init2 (abcorr, s->mpwp);
501 mpc_init2 (mroot, s->mpwp);
502
503 /* Get a copy of the MP coefficients that is local to this thread */
504 while (true && !s->exit_required)
505 {
506 job = mps_thread_job_queue_next (s, data->queue);
507 i = job.i;
508
509 if (job.iter == MPS_THREAD_JOB_EXCEP || *data->nzeros >= s->n)
510 goto cleanup;
511
512 pthread_mutex_lock (&data->roots_mutex[i]);
513
514 if (job.iter == MPS_THREAD_JOB_EXCEP || *data->nzeros >= s->n)
515 {
516 pthread_mutex_unlock (&data->roots_mutex[i]);
517 goto cleanup;
518 }
519
520 /* printf ("Thread %d iterating on root %d\n", data->thread, i); */
521
522 cluster = job.cluster_item->cluster;
523
524 if (s->root[i]->again && !s->root[i]->approximated)
525 {
526 /* Lock this roots to make sure that we are the only one working on it */
527 pthread_mutex_lock (&data->aberth_mutex[i]);
528 mpc_set (mroot, s->root[i]->mvalue);
529 pthread_mutex_unlock (&data->aberth_mutex[i]);
530
531 /* Check if, while we were waiting, excep condition has been reached,
532 * or all the zeros has been approximated. */
533 if ((*data->nzeros) >= s->n)
534 {
535 pthread_mutex_unlock (&data->roots_mutex[i]);
536 goto cleanup;
537 }
538
539 /* pthread_mutex_lock (data->gs_mutex); */
540 (*data->it)++;
541 /* pthread_mutex_unlock (data->gs_mutex); */
542
543 mps_secular_mnewton (s, MPS_POLYNOMIAL (s->secular_equation), s->root[i],
544 corr, mpc_get_prec (s->root[i]->mvalue));
545
546 /* Apply Aberth correction */
547 mps_maberth_s_wl (s, i, cluster, abcorr, data->aberth_mutex);
548 mpc_mul_eq (abcorr, corr);
549 mpc_ui_sub (abcorr, 1U, 0U, abcorr);
550
551 if (!mpc_eq_zero (abcorr))
552 {
553 mpc_div (abcorr, corr, abcorr);
554
555 pthread_mutex_lock (&data->aberth_mutex[i]);
556 mpc_sub_eq (mroot, abcorr);
557 pthread_mutex_unlock (&data->aberth_mutex[i]);
558 }
559 else
560 s->root[i]->again = true;
561
562
563 if (!s->root[i]->again || s->root[i]->approximated)
564 {
565 if (s->debug_level & MPS_DEBUG_APPROXIMATIONS)
566 MPS_DEBUG (s, "Root %d again was set to false on iteration %d by thread %d", i, *data->it, data->thread);
567
568 (*data->nzeros)++;
569 }
570 else
571 {
572 pthread_mutex_lock (&data->aberth_mutex[i]);
573 mpc_set (s->root[i]->mvalue, mroot);
574 pthread_mutex_unlock (&data->aberth_mutex[i]);
575
576 /* Correct the radius */
577 mpc_rmod (modcorr, abcorr);
578 rdpe_add_eq (s->root[i]->drad, modcorr);
579
580 mpc_rmod (modcorr, mroot);
581 rdpe_mul_eq (modcorr, s->mp_epsilon);
582 rdpe_add_eq (s->root[i]->drad, modcorr);
583 }
584 }
585
586 pthread_mutex_unlock (&data->roots_mutex[i]);
587 }
588
589 cleanup:
590 mpc_clear (mroot);
591 mpc_clear (abcorr);
592 mpc_clear (corr);
593
594 return NULL;
595 }
596
597 /**
598 * @brief Routine that performs a block of iteration
599 * in floating point on the secular equation using
600 * CDPE
601 *
602 * @param s the pointer to the mps_context struct.
603 * @param maxit Maximum number of iteration to perform.
604 * @param just_regenerated true if this is the first iteration after a coefficient
605 * regeneration. If just_regenerated is true and the iteration packet is completed
606 * in less than 2 * (n - computed_roots) iterations that best_approx is set to true
607 * in s->secular_equation so a raise in the precision will be triggered.
608 * @return The number of approximated roots after the iteration.
609 */
610 MPS_PRIVATE int
mps_secular_ga_miterate(mps_context * s,int maxit,mps_boolean just_regenerated)611 mps_secular_ga_miterate (mps_context * s, int maxit, mps_boolean just_regenerated)
612 {
613 int computed_roots = 0;
614 int approximated_roots = 0;
615 int root_neighborhood_roots = 0;
616 int i;
617 int nit = 0;
618 int it_threshold = 0;
619
620 s->operation = MPS_OPERATION_ABERTH_MP_ITERATIONS;
621
622 #ifndef DISABLE_DEBUG
623 clock_t *my_clock = mps_start_timer ();
624 #endif
625
626 mps_thread_worker_data *data;
627 pthread_mutex_t *aberth_mutex =
628 (pthread_mutex_t*)mps_malloc (sizeof(pthread_mutex_t) * s->n);
629 pthread_mutex_t *roots_mutex =
630 (pthread_mutex_t*)mps_malloc (sizeof(pthread_mutex_t) * s->n);
631
632 pthread_mutex_t gs_mutex = PTHREAD_MUTEX_INITIALIZER;
633
634 for (i = 0; i < s->n; i++)
635 {
636 pthread_mutex_init (roots_mutex + i, NULL);
637 pthread_mutex_init (aberth_mutex + i, NULL);
638 }
639
640 data = mps_newv (mps_thread_worker_data, s->n_threads);
641
642 MPS_DEBUG_THIS_CALL (s);
643
644 s->best_approx = false;
645
646 it_threshold = s->n - computed_roots;
647
648 /* Mark the approximated roots as ready for output */
649 for (i = 0; i < s->n; i++)
650 {
651 /* Set again to false if the root is already approximated. If a root is approximated but
652 * it has less digits than the current precision don't stop the iterations on that component. */
653 if (s->root[i]->status == MPS_ROOT_STATUS_ISOLATED ||
654 s->root[i]->status == MPS_ROOT_STATUS_APPROXIMATED)
655 {
656 if (s->debug_level & MPS_DEBUG_APPROXIMATIONS)
657 {
658 MPS_DEBUG_WITH_INFO (s, "Setting again[%d] to false since the root is ready for output (or isolated)", i);
659 }
660 s->root[i]->again = false;
661 s->root[i]->approximated = true;
662 }
663
664 if (!s->root[i]->again || s->root[i]->approximated)
665 computed_roots++;
666 }
667
668 mps_thread_job_queue *queue = mps_thread_job_queue_new (s);
669
670 for (i = 0; i < s->n_threads; i++)
671 {
672 data[i].it = &nit;
673 data[i].nzeros = &computed_roots;
674 data[i].s = s;
675 data[i].thread = i;
676 data[i].n_threads = s->n_threads;
677 data[i].aberth_mutex = aberth_mutex;
678 data[i].roots_mutex = roots_mutex;
679 data[i].queue = queue;
680 data[i].gs_mutex = &gs_mutex;
681
682 mps_thread_pool_assign (s, s->pool, __mps_secular_ga_miterate_worker, data + i);
683 }
684
685 mps_thread_pool_wait (s, s->pool);
686
687 /* Check if the roots are improvable in floating point */
688 MPS_DEBUG_WITH_INFO (s, "Performed %d iterations with MP arithmetic",
689 nit);
690
691 if (s->debug_level & MPS_DEBUG_APPROXIMATIONS)
692 mps_dump (s);
693
694 /* Check if we need to get higher precision for the roots */
695 s->best_approx = true;
696 for (i = 0; i < s->n; i++)
697 {
698 if (!s->root[i]->approximated)
699 s->best_approx = false;
700 if (s->root[i]->approximated)
701 approximated_roots++;
702 if (!s->root[i]->again)
703 root_neighborhood_roots++;
704 }
705
706 if (just_regenerated && (nit <= it_threshold))
707 s->best_approx = true;
708
709 MPS_DEBUG_WITH_INFO (s, "%d roots are approximated with the current precision", approximated_roots);
710 MPS_DEBUG_WITH_INFO (s, "%d roots are in the root neighborhood", root_neighborhood_roots);
711 MPS_DEBUG_WITH_INFO (s, "%d roots have reached a stop condition", computed_roots);
712
713 /* These lines are used to debug the again vector, but are not useful
714 * at the moment being */
715 if (s->debug_level & MPS_DEBUG_APPROXIMATIONS)
716 {
717 __MPS_DEBUG (s, "Again vector = ");
718 for (i = 0; i < s->n; i++)
719 {
720 fprintf (s->logstr, "%d ", s->root[i]->again);
721 }
722 fprintf (s->logstr, "\n");
723 }
724
725 /* Clock the routine */
726 #ifndef DISABLE_DEBUG
727 s->mp_iteration_time += mps_stop_timer (my_clock);
728 #endif
729
730 mps_thread_job_queue_free (queue);
731 free (aberth_mutex);
732 free (roots_mutex);
733
734 free (data);
735
736 /* Return the number of approximated roots */
737 return computed_roots;
738 }
739
740