00001
00002
00003
00004 #include "solve_poly_root.h"
00005 gsl_poly_complex_workspace *
00006 gsl_poly_complex_workspace_alloc (size_t n)
00007 {
00008 size_t nc ;
00009
00010 gsl_poly_complex_workspace * w ;
00011
00012 if (n == 0)
00013 {
00014 cpl_msg_error ("gsl_poly_complex_workspace","matrix size n must be positive integer");
00015 return NULL ;
00016 }
00017
00018 w = (gsl_poly_complex_workspace *)
00019 malloc (sizeof(gsl_poly_complex_workspace));
00020
00021 if (w == 0)
00022 {
00023 cpl_msg_error ("gsl_poly_complex_workspace","failed to allocate space for struct");
00024 return NULL ;
00025 }
00026
00027 nc = n - 1;
00028
00029 w->nc = nc;
00030
00031 w->matrix = (double *) malloc (nc * nc * sizeof(double));
00032
00033 if (w->matrix == 0)
00034 {
00035 free (w) ;
00036 cpl_msg_error("gsl_poly_complex_workspace","failed to allocate for workspace matrix") ;
00037 return NULL ;
00038 }
00039
00040 return w ;
00041 }
00042
00043 void
00044 gsl_poly_complex_workspace_free (gsl_poly_complex_workspace * w)
00045 {
00046 free(w->matrix) ;
00047 free(w);
00048 }
00049
00050
00051 int
00052 gsl_poly_complex_solve (const double *a, size_t n,
00053 gsl_poly_complex_workspace * w,
00054 gsl_complex_packed_ptr z)
00055 {
00056 int status;
00057 double *m;
00058
00059 if (n == 0)
00060 {
00061 cpl_msg_error ("gsl_poly_complex_workspace","number of terms must be a positive integer");
00062 return -1 ;
00063 }
00064
00065 if (n == 1)
00066 {
00067 cpl_msg_error ("gsl_poly_complex_workspace","cannot solve for only one term");
00068 return -1 ;
00069 }
00070
00071 if (a[n - 1] == 0)
00072 {
00073 cpl_msg_error ("gsl_poly_complex_workspace","leading term of polynomial must be non-zero") ;
00074 return -1 ;
00075 }
00076
00077 if (w->nc != n - 1)
00078 {
00079 cpl_msg_error ("gsl_poly_complex_workspace","size of workspace does not match polynomial");
00080 return -1 ;
00081 }
00082
00083 m = w->matrix;
00084
00085 set_companion_matrix (a, n - 1, m);
00086
00087 balance_companion_matrix (m, n - 1);
00088
00089 status = qr_companion (m, n - 1, z);
00090
00091 if (status == -1)
00092 {
00093 cpl_msg_error("gsl_poly_complex_workspace","root solving qr method failed to converge");
00094 return -1 ;
00095 }
00096
00097 return 1;
00098 }
00099
00100
00101