/* nag_1d_pade (e02rac) Example Program. * * Copyright 2001 Numerical Algorithms Group. * * Mark 7, 2001. * Mark 7b revised, 2004. */ #include #include #include #include #include int main(void) { /* Scalars */ Integer exit_status, i, l, m, ia, ib, ic; NagError fail; /* Arrays */ double *aa = 0, *bb = 0, *cc = 0, *dd = 0; Complex *z = 0; INIT_FAIL(fail); exit_status = 0; printf("nag_1d_pade (e02rac) Example Program Results\n"); l = 4; m = 4; ia = l + 1; ib = m + 1; ic = ia + ib - 1; /* Allocate memory */ if (!(aa = NAG_ALLOC(ia, double)) || !(bb = NAG_ALLOC(ib, double)) || !(cc = NAG_ALLOC(ic, double)) || !(dd = NAG_ALLOC(ia + ib, double)) || !(z = NAG_ALLOC(l+m, Complex))) { printf("Allocation failure\n"); exit_status = -1; goto END; } /* Power series coefficients in cc */ cc[0] = 1.0; for (i = 1; i <= 8; ++i) cc[i] = cc[i-1] / (double) i; printf("\n"); printf("The given series coefficients are\n"); for (i = 1; i <= ic; ++i) { printf("%13.4e", cc[i-1]); printf(i%5 == 0 || i == ic?"\n":" "); } /* nag_1d_pade (e02rac). * Pade-approximants */ nag_1d_pade(ia, ib, cc, aa, bb, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_1d_pade (e02rac).\n%s\n", fail.message); exit_status = 1; goto END; } printf("\n"); printf("Numerator coefficients\n"); for (i = 1; i <= ia; ++i) { printf("%13.4e", aa[i-1]); printf(i%5 == 0 || i == ia?"\n":" "); } printf("\n"); printf("Denominator coefficients\n"); for (i = 1; i <= ib; ++i) { printf("%13.4e", bb[i-1]); printf(i%5 == 0 || i == ib?"\n":" "); } /* Calculate zeros of the approximant using nag_zeros_real_poly (c02agc) */ /* First need to reverse order of coefficients */ for (i = 1; i <= ia; ++i) dd[ia-i] = aa[i-1]; /* nag_zeros_real_poly (c02agc). * Zeros of a polynomial with real coefficients */ nag_zeros_real_poly(l, dd, Nag_TRUE, z, &fail); if (fail.code != NE_NOERROR) { printf( "Error from nag_zeros_real_poly (c02agc), 1st call.\n%s\n", fail.message); exit_status = 1; goto END; } printf("\n"); printf("Zeros of approximant are at\n"); printf(" Real part Imag part\n"); for (i = 1; i <= l; ++i) printf("%13.4e%13.4e\n", z[i-1].re, z[i-1].im); /* Calculate poles of the approximant using nag_zeros_real_poly (c02agc) */ /* Reverse order of coefficients */ for (i = 1; i <= ib; ++i) dd[ib-i] = bb[i-1]; /* nag_zeros_real_poly (c02agc), see above. */ nag_zeros_real_poly(m, dd, Nag_TRUE, z, &fail); if (fail.code != NE_NOERROR) { printf( "Error from nag_zeros_real_poly (c02agc), 2nd call.\n%s\n", fail.message); exit_status = 1; goto END; } printf("\n"); printf("Poles of approximant are at\n"); printf(" Real part Imag part\n"); for (i = 1; i <= m; ++i) printf("%13.4e%13.4e\n", z[i-1].re, z[i-1].im); END: if (aa) NAG_FREE(aa); if (bb) NAG_FREE(bb); if (cc) NAG_FREE(cc); if (dd) NAG_FREE(dd); if (z) NAG_FREE(z); return exit_status; }