solve_poly_root.c

00001 /*#include "companion.c"*/
00002 /*#include "balance.c"*/
00003 /*#include "qr.c"*/
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) ;       /* error in constructor, avoid memory leak */
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 

Generated on Wed Oct 26 13:08:55 2005 for SINFONI Pipeline Reference Manual by doxygen1.2.13.1 written by Dimitri van Heesch, © 1997-2001