/******************************************************************//*                                                                *//* 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: int@nt.e-technik.uni-erlangen.de                       *//*                                                                *//******************************************************************/#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);}