/* nag_zeros_poly_complex_fpml (c02aac) Example Program.
*
* Copyright 2024 Numerical Algorithms Group.
*
* Mark 30.1, 2024.
*/
#include <nag.h>
#include <stdio.h>
/* Declarations */
int ex1_basic(void);
int ex2_polishing(void);
int main(void) {
int exit_status = 0;
int polish_example = 0;
printf("nag_zeros_poly_complex_fpml (c02aac) Example Program Results\n");
exit_status += ex1_basic();
if (polish_example)
exit_status += ex2_polishing();
return exit_status;
}
int ex1_basic() {
Complex *a = 0, *z = 0;
double *berr = 0, *cond = 0;
Integer *conv = 0, exit_status = 0, i, itmax, n;
/* Nag Types */
Nag_Root_Polish polish;
NagError fail;
printf("\n Basic Problem\n\n");
/* Skip heading in data file */
scanf(" %*[^\n]\n %*[^\n]");
/* Read polynomial degree and allocate */
scanf("%" NAG_IFMT " %*[^\n]", &n);
if (n < 0) {
printf(" Invalid n = %" NAG_IFMT "\n", n);
return -1;
}
if (!(a = NAG_ALLOC(n + 1, Complex)) || !(berr = NAG_ALLOC(n, double)) ||
!(cond = NAG_ALLOC(n, double)) || !(conv = NAG_ALLOC(n, Integer)) ||
!(z = NAG_ALLOC(n, Complex))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Read polynomial coefficients */
for (i = 0; i <= n; i++) {
scanf(" (%lf, %lf)%*[^\n]", &a[i].re, &a[i].im);
}
/* Find roots of the polynomial */
itmax = 30;
polish = Nag_Root_Polish_Simple;
INIT_FAIL(fail);
nag_zeros_poly_complex_fpml(a, n, itmax, polish, z, berr, cond, conv, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_zeros_poly_complex_fpml (c02aac).\n%s\n",
fail.message);
exit_status = -1;
goto END;
}
/* Print output */
printf(" i z conv berr cond\n");
printf(" -----------------------------------------------------\n");
for (i = 0; i < n; i++) {
if (berr[i]<nag_machine_precision) {
printf(" %2" NAG_IFMT " ( %9.2e, %9.2e) %3" NAG_IFMT " < epsilon %9.2e\n",
i, z[i].re, z[i].im, conv[i], cond[i]);
} else {
printf(" %2" NAG_IFMT " ( %9.2e, %9.2e) %3" NAG_IFMT " %9.2e %9.2e\n",
i, z[i].re, z[i].im, conv[i], berr[i], cond[i]);
}
}
END:
NAG_FREE(a);
NAG_FREE(berr);
NAG_FREE(cond);
NAG_FREE(conv);
NAG_FREE(z);
return exit_status;
}
int ex2_polishing() {
const double eps = nag_machine_precision, rmax = nag_machine_real_largest;
const Complex zero = nag_complex_create(0.0, 0.0);
char polish_names[3][7] = {"None ", "Simple", "Comp "};
Complex *a = 0, az, pz, *z = 0, *zact = 0;
double *berr = 0, *cond = 0, delta, err, fwderr, maxfwderr, maxrelerr, relerr;
Integer *conv = 0, exit_status = 0, i, j, k, n, p;
/* Nag Types */
Nag_Boolean *matched = 0;
Nag_Root_Polish polish[3] = {Nag_Root_Polish_None, Nag_Root_Polish_Simple,
Nag_Root_Polish_Compensated};
NagError fail;
printf("\n Polishing Modes\n\n");
/* Skip heading in file */
scanf(" %*[^\n]");
/* Read polynomial degree and allocate */
scanf("%" NAG_IFMT " %*[^\n]", &n);
if (n < 0) {
printf(" Invalid n = %" NAG_IFMT "\n", n);
return -1;
}
if (!(a = NAG_ALLOC(n + 1, Complex)) || !(berr = NAG_ALLOC(n, double)) ||
!(cond = NAG_ALLOC(n, double)) || !(conv = NAG_ALLOC(n, Integer)) ||
!(matched = NAG_ALLOC(n, Nag_Boolean)) || !(z = NAG_ALLOC(n, Complex)) ||
!(zact = NAG_ALLOC(n, Complex))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Set expected roots */
for (i = 0; i < n; i++)
zact[i] = nag_complex_create(i + 1, 0.0);
/* Multiply out (z-1)...(z-n) for coefficients */
INIT_FAIL(fail);
nag_blast_zload(n, zero, a, 1, &fail);
a[n] = nag_complex_create(1.0, 0.0);
for (i = 0; i < n; i++) {
for (j = n - i; j <= n; j++) {
az = nag_complex_multiply(a[j - 1], zact[i]);
a[j - 1] = nag_complex_subtract(a[j], az);
}
az = nag_complex_multiply(a[n], zact[i]);
a[n] = nag_complex_negate(az);
}
printf(" polish relerr fwderr\n");
printf(" ----------------------------\n");
for (p = 0; p <= 2; p++) {
/* Find roots */
INIT_FAIL(fail);
nag_zeros_poly_complex_fpml(a, n, 30, polish[p], z, berr, cond, conv,
&fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_zeros_poly_complex_fpml (c02aac).\n%s\n",
fail.message);
exit_status = -1;
goto END;
}
/* Calculate the maximum relative errors of the roots, and the maximum
* forward error evaluating the polynomial at those roots. Errors are
* capped at machine precision. */
maxfwderr = maxrelerr = eps;
for (i = 0; i < n; i++)
matched[i] = Nag_FALSE;
for (i = 0; i < n; i++) {
/* Evaluate polynomial at this root */
pz = a[0];
for (j = 1; j <= n; j++) {
pz = nag_complex_multiply(pz, z[i]);
pz = nag_complex_add(pz, a[j]);
}
/* Match to an expected root */
k = 0;
err = rmax;
for (j = 0; j < n; j++) {
if (!matched[j]) {
delta = nag_complex_abs(nag_complex_subtract(z[i], zact[j]));
if (delta <= err) {
err = delta;
k = j;
}
}
}
/* Mark as matched and update max errors */
matched[k] = Nag_TRUE;
relerr = err / nag_complex_abs(zact[k]);
fwderr = nag_complex_abs(pz);
if (relerr > maxrelerr)
maxrelerr = relerr;
if (fwderr > maxfwderr)
maxfwderr = fwderr;
}
/* Print output */
printf(" %s%10.2e%10.2e\n", polish_names[p], maxrelerr, maxfwderr);
}
END:
NAG_FREE(a);
NAG_FREE(berr);
NAG_FREE(cond);
NAG_FREE(conv);
NAG_FREE(matched);
NAG_FREE(z);
NAG_FREE(zact);
return exit_status;
}