rootsl.c

Polynomial Root Finder is a reliable and fast C program (+ Matlab gateway) for finding all roots of

/******************************************************************/
/*                                                                */
/* file:          rootsl.c                                        */
/*                                                                */
/* main function: mexFunction()                                   */
/*                                                                */
/* version:       1.1                                             */
/*                                                                */
/* author:        B. Frenzel                                      */
/*                                                                */
/* date:          Jan 7 1993                                      */ 
/*                                                                */
/* input:         nrhs    number of right hand side arguments     */
/*                *prhs[] Length nrhs array for ptrs to           */
/*                        rhs arguments                           */
/* output:        nlhs    number of left hand side arguments      */
/*                *plhs[] Length nlhs array for ptrs to           */
/*                        lhs arguments                           */
/*                                                                */
/* subroutines:   mxGetM(),mxGetN(),mxGetPr(),mxGetPi(),          */
/*                mxSetPr(),mxCreateFull(),mexErrMsgTxt(),        */ 
/*                null()                                          */
/* compilation (on HP 700):                                       */
/* cmex -O -Aa rootsl.c null.c complex.c newton.c muller.c tools.c*/  
/*                                                                */
/* description:                                                   */
/* gateway routine for matlab version 4.0; invokes main rootfind. */
/* routine using the  Muller method followed by the               */
/* Newton method                                                  */
/*                                                                */
/* Copyright:                                                     */
/* Lehrstuhl fuer Nachrichtentechnik Erlangen                     */
/* Cauerstr. 7, 8520 Erlangen, FRG, 1993                          */
/* e-mail: [e-mail]                       */
/*                                                                */
/******************************************************************/
#define  GATE40
#include "header.h"







/***** gateway routine *****/
void mexFunction(int nlhs,Matrix *plhs[],int nrhs,Matrix *prhs[])
/* int    nlhs,     number of left hand side arguments          */ 
/*        nrhs;     number of right hand side arguments         */
/* Matrix *plhs[],  Length nlhs array for ptrs to lhs arguments */
/*        *prhs[];  Length nrhs array for ptrs to rhs arguments */
{
     int           i;          /* counter                             */
     int           m,          /* row dimension of prhs[0]            */
                   n;          /* column dimension of prhs[0]         */
     double        maxerr,     /* max. determined error of roots in   */
                               /* null()                              */
                   *prr,       /* ptr to real part of prhs[0]         */
                   *pri,       /* ptr to imag. part of prhs[0]        */
                   *rootr,     /* ptr to real part of plhs[0]         */
                   *rooti,     /* ptr to imag. part of plhs[0]        */
                   *prf;       /* ptr to prhs[1]                      */
     dcomplex      p[MAXCOEFF],   /* coefficient vector of original   */
                                  /* polynomial for null()            */ 
                   pred[MAXCOEFF],/* coefficient vector of deflated   */
                                  /* polynomial for null()            */
                   root[MAXCOEFF];/* vector of determined roots       */
                                  /* for null()                       */
     unsigned char flag;   /* flag = TRUE  => complex coefficients    */
                           /* flag = FALSE => real    coefficients    */
     char error;           /* indicates runtime error in null()       */

                           /* check for proper number of arguments    */
     if (nrhs > 2) {
          mexErrMsgTxt("One or two input arguments required.");
     } else if (nlhs > 2) {
          mexErrMsgTxt("One or two output arguments required.");
     }

                           /* coefficients vector in prhs[0]          */
     m = mxGetM(prhs[0]);  /* row dimension of prhs[0]                */
     n = mxGetN(prhs[0]);  /* column dimension of prhs[0]             */

     if (m>n)              /* check whether input is a row */
     n = m; m = 1;

     n--;                  /* degree of polynomial = column dimension */
                           /* of prhs[0] - 1                          */
     prr = mxGetPr(prhs[0]);
     pri = mxGetPi(prhs[0]);
     if (pri==NLL)        
          flag = FALSE;    /* real coefficients                       */
     else
          flag = TRUE;     /* complex coefficients                    */
                           /* store coefficients in p[i]              */
     for (i=0; i<=n; i++) {
          p[n-i].r = *prr;
          p[n-i].i = (flag==TRUE) ? *pri : 0.;
          prr++;
          pri++;
     }



                           /* work with complex algorithm?            */
     if (nrhs==2) {
          prf = mxGetPr(prhs[1]);
          if (*prf=='C' || *prf=='c')
               flag=TRUE;
     }


                           /* invoke main rootfinding routine         */
     error = null(p,pred,&n,root,&maxerr,flag);

                           /* create matrix for the return arguments  */
     plhs[0] = mxCreateFull(n,m,COMPLEX); /* for determined roots     */
     plhs[1] = mxCreateFull(1,1,REAL);    /* for max. error           */
  

     if (error) {          /* runtime error in null() occured         */
          switch (error) {
          case 1:
               mexErrMsgTxt("Power of polynomial lower null.");
               break;
          case 2:
               mexErrMsgTxt("Polynomial is a null vector.");
               break;
          case 3:
               mexErrMsgTxt("Polynomial is a constant unequal null.");
               break;
          }
     }                    
 
                           /* store determined roots in plhs[0]       */
     rootr = mxGetPr(plhs[0]);          
     rooti = mxGetPi(plhs[0]);          
     for (i=0;i<n;i++) {
          rootr[i] = root[i].r;
          rooti[i] = root[i].i;
     }
                           /* store maxerr in plhs[1]                 */
     mxSetPr(plhs[1],&maxerr);
}