Implements Durand Kerner iterative algorithm to compute all roots of a polynomial.
111 {
112 long double tol_condition = 1;
113 uint32_t iter = 0;
114 int n;
116
117 if (write_log) {
118
119
120
121 log_file.
open(
"durand_kerner.log.csv");
123 perror(
"Unable to create a storage log file!");
125 }
126 log_file << "iter#,";
127
128 for (n = 0; n < roots->size(); n++) log_file << "root_" << n << ",";
129
130 log_file << "avg. correction";
131 log_file << "\n0,";
132 for (n = 0; n < roots->size(); n++)
134 }
135
136 bool break_loop = false;
138 !break_loop) {
139 tol_condition = 0;
140 iter++;
141 break_loop = false;
142
144 log_file << "\n" << iter << ",";
145
146#ifdef _OPENMP
147#pragma omp parallel for shared(break_loop, tol_condition)
148#endif
149 for (n = 0; n < roots->size(); n++) {
150 if (break_loop)
151 continue;
152
155 denominator = 1.0;
156 for (int i = 0; i < roots->size(); i++)
157 if (i != n)
158 denominator *= (*roots)[n] - (*roots)[i];
159
161
163 std::cerr <<
"\n\nOverflow/underrun error - got value = "
164 << std::abs(delta) << "\n";
165
166 break_loop = true;
167 }
168
169 (*roots)[n] -= delta;
170
171#ifdef _OPENMP
172#pragma omp critical
173#endif
174 tol_condition =
std::max(tol_condition, std::abs(std::abs(delta)));
175 }
176
177
178 if (break_loop)
179 break;
180
182 for (n = 0; n < roots->size(); n++)
184 }
185
186#if defined(DEBUG) || !defined(NDEBUG)
187 if (iter % 500 == 0) {
189 for (n = 0; n < roots->size(); n++)
191 std::cout <<
"\t\tabsolute average change: " << tol_condition
192 << "\n";
193 }
194#endif
195
197 log_file << tol_condition;
198 }
199
201}
bool check_termination(long double delta)
Definition durand_kerner_roots.cpp:91
const char * complex_str(const std::complex< double > &x)
Definition durand_kerner_roots.cpp:76
std::complex< double > poly_function(const std::valarray< double > &coeffs, std::complex< double > x)
Definition durand_kerner_roots.cpp:53