ida.c

Go to the documentation of this file.

/*
 * -----------------------------------------------------------------
 * $Revision: 1.3 $
 * $Date: 2006/10/05 22:09:14 $
 * ----------------------------------------------------------------- 
 * Programmer(s): Alan Hindmarsh, Radu Serban and Aaron Collier @ LLNL
 * -----------------------------------------------------------------
 * Copyright (c) 2002, The Regents of the University of California.
 * Produced at the Lawrence Livermore National Laboratory.
 * All rights reserved.
 * For details, see the LICENSE file.
 * -----------------------------------------------------------------
 * This is the implementation file for the main IDA solver.
 * It is independent of the linear solver in use.
 * -----------------------------------------------------------------
 *
 * EXPORTED FUNCTIONS
 * ------------------
 *   Creation, allocation and re-initialization functions
 *       IDACreate
 *       IDAMalloc
 *       IDAReInit
 *       IDARootInit
 *   Main solver function
 *       IDASolve
 *   Interpolated output and extraction functions
 *       IDAGetSolution
 *   Deallocation functions
 *       IDAFree
 *
 * PRIVATE FUNCTIONS
 * -----------------
 *       IDACheckNvector
 *   Memory allocation/deallocation
 *       IDAAllocVectors
 *       IDAFreeVectors
 *   Initial setup
 *       IDAInitialSetup
 *       IDAEwtSet
 *       IDAEwtSetSS
 *       IDAEwtSetSV
 *   Stopping tests
 *       IDAStopTest1
 *       IDAStopTest2
 *   Error handler
 *       IDAHandleFailure
 *   Main IDAStep function
 *       IDAStep
 *       IDASetCoeffs
 *   Nonlinear solver functions
 *       IDANls
 *       IDAPredict
 *       IDANewtonIter
 *   Error test
 *       IDATestError
 *       IDARestore
 *   Handler for convergence and/or error test failures
 *       IDAHandleNFlag
 *       IDAReset
 *   Function called after a successful step
 *       IDACompleteStep
 *   Norm functions
 *       IDAWrmsNorm
 *   Functions for rootfinding
 *       IDARcheck1
 *       IDARcheck2
 *       IDARcheck3
 *       IDARootfind
 *   IDA Error message handling functions 
 *       IDAProcessError
 *       IDAErrHandler
 * -----------------------------------------------------------------
 */

/* 
 * =================================================================
 * IMPORTED HEADER FILES
 * =================================================================
 */

#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>

#include "ida_impl.h"
#include "sundials_math.h"

/* 
 * =================================================================
 * MACRO DEFINITIONS
 * =================================================================
 */

/* Macro: loop */
#define loop for(;;)

/* 
 * =================================================================
 * IDAS PRIVATE CONSTANTS
 * =================================================================
 */

#define ZERO      RCONST(0.0)    /* real 0.0    */
#define HALF      RCONST(0.5)    /* real 0.5    */
#define QUARTER   RCONST(0.25)   /* real 0.25   */
#define TWOTHIRDS RCONST(0.667)  /* real 2/3    */
#define ONE       RCONST(1.0)    /* real 1.0    */
#define ONEPT5    RCONST(1.5)    /* real 1.5    */
#define TWO       RCONST(2.0)    /* real 2.0    */
#define FOUR      RCONST(4.0)    /* real 4.0    */
#define FIVE      RCONST(5.0)    /* real 5.0    */
#define TEN       RCONST(10.0)   /* real 10.0   */
#define TWELVE    RCONST(12.0)   /* real 12.0   */
#define TWENTY    RCONST(20.0)   /* real 20.0   */
#define HUNDRED   RCONST(100.0)  /* real 100.0  */
#define PT9       RCONST(0.9)    /* real 0.9    */
#define PT99      RCONST(0.99)   /* real 0.99   */
#define PT1       RCONST(0.1)    /* real 0.1    */
#define PT01      RCONST(0.01)   /* real 0.01   */
#define PT001     RCONST(0.001)  /* real 0.001  */
#define PT0001    RCONST(0.0001) /* real 0.0001 */

/* 
 * =================================================================
 * IDAS ROUTINE-SPECIFIC CONSTANTS
 * =================================================================
 */

/* 
 * Control constants for lower-level functions used by IDASolve 
 * ------------------------------------------------------------
 */

/* IDAStep control constants */

#define PREDICT_AGAIN 20

/* Return values for lower level routines used by IDASolve */

#define IDA_RES_RECVR    +1
#define IDA_LSETUP_RECVR +2
#define IDA_LSOLVE_RECVR +3

#define IDA_NCONV_RECVR  +4
#define IDA_CONSTR_RECVR +5
#define CONTINUE_STEPS   +99

/* IDACompleteStep constants */

#define UNSET            -1
#define LOWER            +1 
#define RAISE            +2 
#define MAINTAIN         +3

/* IDATestError constants */

#define ERROR_TEST_FAIL  +7

/*
 * Control constants for lower-level rootfinding functions
 * -------------------------------------------------------
 */

#define RTFOUND          +1
#define INITROOT         +2
#define CLOSERT          +3
#define ZERODETACHING    +4
#define MASKED           55
/*
 * Algorithmic constants
 * ---------------------
 */

#define MXNCF           10  /* max number of convergence failures allowed */
#define MXNEF           10  /* max number of error test failures allowed  */
#define MAXNH            5  /* max. number of h tries in IC calc. */
#define MAXNJ            4  /* max. number of J tries in IC calc. */
#define MAXNI           10  /* max. Newton iterations in IC calc. */
#define EPCON RCONST(0.33)  /* Newton convergence test constant */

/* IDANewtonIter constants */

#define MAXIT   4
#define RATEMAX RCONST(0.9)
#define XRATE   RCONST(0.25)        

/* 
 * =================================================================
 * PRIVATE FUNCTION PROTOTYPES
 * =================================================================
 */

static booleantype IDACheckNvector(N_Vector tmpl);

/* Memory allocation/deallocation */

static booleantype IDAAllocVectors(IDAMem IDA_mem, N_Vector tmpl, int tol);
static void IDAFreeVectors(IDAMem IDA_mem);

/* Initial setup */

int IDAInitialSetup(IDAMem IDA_mem);
static int IDAEwtSetSS(IDAMem IDA_mem, N_Vector ycur, N_Vector weight);
static int IDAEwtSetSV(IDAMem IDA_mem, N_Vector ycur, N_Vector weight);

/* Main IDAStep function */

static int IDAStep(IDAMem IDA_mem);

/* Function called at beginning of step */

static void IDASetCoeffs(IDAMem IDA_mem, realtype *ck);

/* Nonlinear solver functions */

static void IDAPredict(IDAMem IDA_mem);
static int IDANls(IDAMem IDA_mem);
static int IDANewtonIter(IDAMem IDA_mem);

/* Error test */

static int IDATestError(IDAMem IDA_mem, realtype ck, 
                        realtype *err_k, realtype *err_km1);

/* Handling of convergence and/or error test failures */

static void IDARestore(IDAMem IDA_mem, realtype saved_t);
static int IDAHandleNFlag(IDAMem IDA_mem, int nflag, realtype err_k, realtype err_km1,
                          long int *ncfnPtr, int *ncfPtr, long int *netfPtr, int *nefPtr);
static void IDAReset(IDAMem IDA_mem);

/* Function called after a successful step */

static void IDACompleteStep(IDAMem IDA_mem, realtype err_k, realtype err_km1);

/* Stopping tests and failure handling */

static int IDAStopTest1(IDAMem IDA_mem, realtype tout,realtype *tret, 
                        N_Vector yret, N_Vector ypret, int itask);
static int IDAStopTest2(IDAMem IDA_mem, realtype tout, realtype *tret, 
                        N_Vector yret, N_Vector ypret, int itask);
static int IDAHandleFailure(IDAMem IDA_mem, int sflag);

/* Functions for rootfinding */

static int IDARcheck1(IDAMem IDA_mem);
static int IDARcheck2(IDAMem IDA_mem);
static int IDARcheck3(IDAMem IDA_mem);
static int IDARootfind(IDAMem IDA_mem);

/* Norm functions */

realtype IDAWrmsNorm(IDAMem IDA_mem, N_Vector x, N_Vector w, booleantype mask);

/* 
 * =================================================================
 * EXPORTED FUNCTIONS IMPLEMENTATION
 * =================================================================
 */

/* 
 * -----------------------------------------------------------------
 * Creation, allocation and re-initialization functions
 * -----------------------------------------------------------------
 */

/* 
 * IDACreate
 *
 * IDACreate creates an internal memory block for a problem to 
 * be solved by IDA.
 * If successful, IDACreate returns a pointer to the problem memory. 
 * This pointer should be passed to IDAMalloc.  
 * If an initialization error occurs, IDACreate prints an error 
 * message to standard err and returns NULL. 
 */

void *IDACreate(void)
{
  IDAMem IDA_mem;

  IDA_mem = NULL;
  IDA_mem = (IDAMem) malloc(sizeof(struct IDAMemRec));
  if (IDA_mem == NULL) {
    IDAProcessError(NULL, 0, "IDA", "IDACreate", MSG_MEM_FAIL);
    return (NULL);
  }

  /* Set unit roundoff in IDA_mem */
  IDA_mem->ida_uround = UNIT_ROUNDOFF;

  /* Set default values for integrator optional inputs */
  IDA_mem->ida_res         = NULL;
  IDA_mem->ida_rdata       = NULL;
  IDA_mem->ida_efun        = NULL;
  IDA_mem->ida_edata       = NULL;
  IDA_mem->ida_ehfun       = IDAErrHandler;
  IDA_mem->ida_eh_data     = (void *) IDA_mem;
  IDA_mem->ida_errfp       = stderr;
  IDA_mem->ida_maxord      = MAXORD_DEFAULT;
  IDA_mem->ida_mxstep      = MXSTEP_DEFAULT;
  IDA_mem->ida_hmax_inv    = HMAX_INV_DEFAULT;
  IDA_mem->ida_hin         = ZERO;
  IDA_mem->ida_epcon       = EPCON;
  IDA_mem->ida_maxnef      = MXNEF;
  IDA_mem->ida_maxncf      = MXNCF;
  IDA_mem->ida_maxcor      = MAXIT;
  IDA_mem->ida_suppressalg = FALSE;
  IDA_mem->ida_id          = NULL;
  IDA_mem->ida_constraints = NULL;
  IDA_mem->ida_constraintsSet = FALSE;
  IDA_mem->ida_tstopset    = FALSE;

  /* set the saved value maxord_alloc */
  IDA_mem->ida_maxord_alloc = MAXORD_DEFAULT;

  /* Set default values for IC optional inputs */
  IDA_mem->ida_epiccon = PT01 * EPCON;
  IDA_mem->ida_maxnh   = MAXNH;
  IDA_mem->ida_maxnj   = MAXNJ;
  IDA_mem->ida_maxnit  = MAXNI;
  IDA_mem->ida_lsoff   = FALSE;
  IDA_mem->ida_steptol = RPowerR(IDA_mem->ida_uround, TWOTHIRDS);

  /* Initialize lrw and liw */
  IDA_mem->ida_lrw = 25 + 5*MXORDP1;
  IDA_mem->ida_liw = 38;

  /* No mallocs have been done yet */
  IDA_mem->ida_VatolMallocDone = FALSE;
  IDA_mem->ida_constraintsMallocDone = FALSE;
  IDA_mem->ida_idMallocDone = FALSE;
  IDA_mem->ida_MallocDone = FALSE;

  /* Return pointer to IDA memory block */
  return((void *)IDA_mem);
}

/*-----------------------------------------------------------------*/

#define lrw   (IDA_mem->ida_lrw)
#define liw   (IDA_mem->ida_liw)

/*-----------------------------------------------------------------*/

/*
 * IDAMalloc
 *
 * IDAMalloc allocates and initializes memory for a problem. All
 * problem specification inputs are checked for errors. If any
 * error occurs during initialization, it is reported to the 
 * error handler function.
 */

int IDAMalloc(void *ida_mem, IDAResFn res,
              realtype t0, N_Vector yy0, N_Vector yp0, 
              int itol, realtype rtol, void *atol)
{
  IDAMem IDA_mem;
  booleantype nvectorOK, allocOK, neg_atol;
  long int lrw1, liw1;

  /* Check ida_mem */
  if (ida_mem == NULL) {
    IDAProcessError(NULL, IDA_MEM_NULL, "IDA", "IDAMalloc", MSG_NO_MEM);
    return(IDA_MEM_NULL);
  }
  IDA_mem = (IDAMem) ida_mem;
  
  /* Check for legal input parameters */
  
  if (yy0 == NULL) { 
    IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDAMalloc", MSG_Y0_NULL);
    /*  return(IDA_ILL_INPUT); Masoud */
    return(IDA_NULL_Y0); 

  }
  
  if (yp0 == NULL) { 
    IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDAMalloc", MSG_YP0_NULL);
    /*  return(IDA_ILL_INPUT); Masoud */
    return(IDA_YP0_NULL); 
  }

  if ((itol != IDA_SS) && (itol != IDA_SV) && (itol != IDA_WF)) {
    IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDAMalloc", MSG_BAD_ITOL);
    /*  return(IDA_ILL_INPUT); Masoud */
    return(IDA_BAD_ITOL);
  }

  if (res == NULL) { 
    IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDAMalloc", MSG_RES_NULL);
    /*  return(IDA_ILL_INPUT); Masoud */
    return(IDA_RES_NULL); 
  }

  /* Test if all required vector operations are implemented */
  nvectorOK = IDACheckNvector(yy0);
  if (!nvectorOK) {
    IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDAMalloc", MSG_BAD_NVECTOR);
    /*  return(IDA_ILL_INPUT); Masoud */
    return(IDA_BAD_NVECTOR);
  }

  /* Test tolerances */

  if (itol != IDA_WF) {

    if (atol == NULL) { 
      IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDAMalloc", MSG_ATOL_NULL);
    /*  return(IDA_ILL_INPUT); Masoud */
      return(IDA_NULL_ABSTOL); 
    }

    if (rtol < ZERO) { 
      IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDAMalloc", MSG_BAD_RTOL);
    /*  return(IDA_ILL_INPUT); Masoud */
      return(IDA_BAD_RELTOL); 
    }
   
    if (itol == IDA_SS) { 
      neg_atol = (*((realtype *)atol) < ZERO); 
    } else { 
      neg_atol = (N_VMin((N_Vector)atol) < ZERO); 
    }

    if (neg_atol) { 
      IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDAMalloc", MSG_BAD_ATOL);
      /*  return(IDA_ILL_INPUT); Masoud */
      return(IDA_BAD_ABSTOL); 
    }

  }

  /* Set space requirements for one N_Vector */
  if (yy0->ops->nvspace != NULL) {
    N_VSpace(yy0, &lrw1, &liw1);
  } else {
    lrw1 = 0;
    liw1 = 0;
  }
  IDA_mem->ida_lrw1 = lrw1;
  IDA_mem->ida_liw1 = liw1;

  /* Allocate the vectors (using yy0 as a template) */
  allocOK = IDAAllocVectors(IDA_mem, yy0, itol);
  if (!allocOK) {
    IDAProcessError(IDA_mem, IDA_MEM_FAIL, "IDA", "IDAMalloc", MSG_MEM_FAIL);
    return(IDA_MEM_FAIL);
  }
 
  /* All error checking is complete at this point */

  /* Copy the input parameters into IDA memory block */

  IDA_mem->ida_res = res;
  IDA_mem->ida_tn  = t0;

  /* Copy tolerances into memory */

  IDA_mem->ida_itol = itol;
  IDA_mem->ida_rtol = rtol;      

  if (itol == IDA_SS)
    IDA_mem->ida_Satol = *((realtype *)atol);
  else if (itol == IDA_SV) 
    N_VScale(ONE, (N_Vector)atol, IDA_mem->ida_Vatol);

  /* Set the linear solver addresses to NULL */
  IDA_mem->ida_linit  = NULL;
  IDA_mem->ida_lsetup = NULL;
  IDA_mem->ida_lsolve = NULL;
  IDA_mem->ida_lperf  = NULL;
  IDA_mem->ida_lfree  = NULL;
  IDA_mem->ida_lmem   = NULL;

  /* Initialize the phi array */
  N_VScale(ONE, yy0, IDA_mem->ida_phi[0]);  
  N_VScale(ONE, yp0, IDA_mem->ida_phi[1]);  
 
  /* Initialize all the counters and other optional output values */
  IDA_mem->ida_nst     = 0;
  IDA_mem->ida_nre     = 0;
  IDA_mem->ida_ncfn    = 0;
  IDA_mem->ida_netf    = 0;
  IDA_mem->ida_nni     = 0;
  IDA_mem->ida_nsetups = 0;
  
  IDA_mem->ida_kused = 0;
  IDA_mem->ida_hused = ZERO;
  IDA_mem->ida_tolsf = ONE;

  IDA_mem->ida_nge = 0;

  /* Initialize root-finding variables */

  IDA_mem->ida_glo    = NULL;
  IDA_mem->ida_ghi    = NULL;
  IDA_mem->ida_grout  = NULL;
  IDA_mem->ida_iroots = NULL;
  IDA_mem->ida_gfun   = NULL;
  IDA_mem->ida_g_data = NULL;
  IDA_mem->ida_nrtfn  = 0;

  /* Initial setup not done yet */
  IDA_mem->ida_SetupDone = FALSE;

  /* Problem memory has been successfully allocated */
  IDA_mem->ida_MallocDone = TRUE;
  return(IDA_SUCCESS);
}

/*-----------------------------------------------------------------*/

#define lrw1 (IDA_mem->ida_lrw1)
#define liw1 (IDA_mem->ida_liw1)

/*-----------------------------------------------------------------*/

/*
 * IDAReInit
 *
 * IDAReInit re-initializes IDA's memory for a problem, assuming
 * it has already beeen allocated in a prior IDAMalloc call.
 * All problem specification inputs are checked for errors.
 * The problem size Neq is assumed to be unchaged since the call
 * to IDAMalloc, and the maximum order maxord must not be larger.
 * If any error occurs during reinitialization, it is reported to
 * the error handler function.
 * The return value is IDA_SUCCESS = 0 if no errors occurred, or
 * a negative value otherwise.
 */

int IDAReInit(void *ida_mem, IDAResFn res,
              realtype t0, N_Vector yy0, N_Vector yp0,
              int itol, realtype rtol, void *atol)
{
  IDAMem IDA_mem;
  booleantype neg_atol;

  /* Check for legal input parameters */
  
  if (ida_mem == NULL) {
    IDAProcessError(NULL, IDA_MEM_NULL, "IDA", "IDAReInit", MSG_NO_MEM);
    return(IDA_MEM_NULL);
  }
  IDA_mem = (IDAMem) ida_mem;

  /* Check if problem was malloc'ed */
  
  if (IDA_mem->ida_MallocDone == FALSE) {
    IDAProcessError(IDA_mem, IDA_NO_MALLOC, "IDA", "IDAReInit", MSG_NO_MALLOC);
    return(IDA_NO_MALLOC);
  }

  /* Check for legal input parameters */
  
  if (yy0 == NULL) { 
    IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDAReInit", MSG_Y0_NULL);
    /*  return(IDA_ILL_INPUT); Masoud */
    return(IDA_NULL_Y0); 
  }
  
  if (yp0 == NULL) { 
    IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDAReInit", MSG_YP0_NULL);
    /*  return(IDA_ILL_INPUT); Masoud */
    return(IDA_YP0_NULL); 
  }

  if ((itol != IDA_SS) && (itol != IDA_SV) && (itol != IDA_WF)) {
    IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDAReInit", MSG_BAD_ITOL);
    /*  return(IDA_ILL_INPUT); Masoud */
    return(IDA_BAD_ITOL);
  }

  if (res == NULL) { 
    IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDAReInit", MSG_RES_NULL);
    /*  return(IDA_ILL_INPUT); Masoud */
    return(IDA_RES_NULL); 
  }

  /* Test tolerances */

  if (itol != IDA_WF) {

    if (atol == NULL) { 
      IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDAReInit", MSG_ATOL_NULL);
    /*  return(IDA_ILL_INPUT); Masoud */
      return(IDA_NULL_ABSTOL); 
    }
    
    if (rtol < ZERO) {
      IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDAReInit", MSG_BAD_RTOL);
    /*  return(IDA_ILL_INPUT); Masoud */
      return(IDA_BAD_RELTOL); 
    }
   
    if (itol == IDA_SS) { 
      neg_atol = (*((realtype *)atol) < ZERO); 
    } else { 
      neg_atol = (N_VMin((N_Vector)atol) < ZERO); 
    }
    if (neg_atol) { 
      IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDAReInit", MSG_BAD_ATOL);
    /*  return(IDA_ILL_INPUT); Masoud */
      return(IDA_BAD_ABSTOL); 
    }

  }

  /* Copy the input parameters into IDA memory block */
  IDA_mem->ida_res = res;
  IDA_mem->ida_tn  = t0;

  if ( (itol != IDA_SV) && (IDA_mem->ida_VatolMallocDone) ) {
    N_VDestroy(IDA_mem->ida_Vatol);
    lrw -= lrw1;
    liw -= liw1;
    IDA_mem->ida_VatolMallocDone = FALSE;
  }

  if ( (itol == IDA_SV) && !(IDA_mem->ida_VatolMallocDone) ) {
    IDA_mem->ida_Vatol = NULL;
    IDA_mem->ida_Vatol = N_VClone(yy0);
    lrw += lrw1;
    liw += liw1;
    IDA_mem->ida_VatolMallocDone = TRUE;
  }

  IDA_mem->ida_itol = itol;
  IDA_mem->ida_rtol = rtol;      
  if (itol == IDA_SS)
    IDA_mem->ida_Satol = *((realtype *)atol);
  else if (itol == IDA_SV)
    N_VScale(ONE, (N_Vector)atol, IDA_mem->ida_Vatol);

  /* Initialize the phi array */
  N_VScale(ONE, yy0, IDA_mem->ida_phi[0]);  
  N_VScale(ONE, yp0, IDA_mem->ida_phi[1]);  
 
  /* Initialize all the counters and other optional output values */
 
  IDA_mem->ida_nst     = 0;
  IDA_mem->ida_nre     = 0;
  IDA_mem->ida_ncfn    = 0;
  IDA_mem->ida_netf    = 0;
  IDA_mem->ida_nni     = 0;
  IDA_mem->ida_nsetups = 0;
  
  IDA_mem->ida_kused = 0;
  IDA_mem->ida_hused = ZERO;
  IDA_mem->ida_tolsf = ONE;

  IDA_mem->ida_nge = 0;

  /* Initial setup not done yet */
  IDA_mem->ida_SetupDone = FALSE;
      
  /* Problem has been successfully re-initialized */

  return(IDA_SUCCESS);
}

/*-----------------------------------------------------------------*/

#define gfun   (IDA_mem->ida_gfun)
#define g_data (IDA_mem->ida_g_data) 
#define glo    (IDA_mem->ida_glo)
#define ghi    (IDA_mem->ida_ghi)
#define grout  (IDA_mem->ida_grout)
#define iroots (IDA_mem->ida_iroots)

/*-----------------------------------------------------------------*/

/*
 * IDARootInit
 *
 * IDARootInit initializes a rootfinding problem to be solved
 * during the integration of the DAE system.  It loads the root
 * function pointer and the number of root functions, and allocates
 * workspace memory.  The return value is IDA_SUCCESS = 0 if no
 * errors occurred, or a negative value otherwise.
 */

int IDARootInit(void *ida_mem, int nrtfn, IDARootFn g, void *gdata)
{
  IDAMem IDA_mem;
  int nrt;

  /* Check ida_mem pointer */
  if (ida_mem == NULL) {
    IDAProcessError(NULL, IDA_MEM_NULL, "IDA", "IDARootInit", MSG_NO_MEM);
    return(IDA_MEM_NULL);
  }
  IDA_mem = (IDAMem) ida_mem;

  nrt = (nrtfn < 0) ? 0 : nrtfn;

  /* If rerunning IDARootInit() with a different number of root
     functions (changing number of gfun components), then free
     currently held memory resources */
  if ((nrt != IDA_mem->ida_nrtfn) && (IDA_mem->ida_nrtfn > 0)) {

    free(glo); glo = NULL;
    free(ghi); ghi = NULL;
    free(grout); grout = NULL;
    free(iroots); iroots = NULL;

    lrw -= 3* (IDA_mem->ida_nrtfn);
    liw -= IDA_mem->ida_nrtfn;

  }

  /* If IDARootInit() was called with nrtfn == 0, then set ida_nrtfn to
     zero and ida_gfun to NULL before returning */
  if (nrt == 0) {
    IDA_mem->ida_nrtfn = nrt;
    gfun = NULL;
    g_data = NULL;
    return(IDA_SUCCESS);
  }

  /* Store user's data pointer */
  g_data = gdata;

  /* If rerunning IDARootInit() with the same number of root functions
     (not changing number of gfun components), then check if the root
     function argument has changed */
  /* If g != NULL then return as currently reserved memory resources
     will suffice */
  if (nrt == IDA_mem->ida_nrtfn) {
    if (g != gfun) {
      if (g == NULL) {
        free(glo); glo = NULL;
        free(ghi); ghi = NULL;
        free(grout); grout = NULL;
        free(iroots); iroots = NULL;

        lrw -= 3*nrt;
        liw -= nrt;

        IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDARootInit", MSG_ROOT_FUNC_NULL);
        /*  return(IDA_ILL_INPUT); Masoud */
        return(IDA_NULL_G);
      }
      else {
        gfun = g;
        return(IDA_SUCCESS);
      }
    }
    else return(IDA_SUCCESS);
  }

  /* Set variable values in IDA memory block */
  IDA_mem->ida_nrtfn = nrt;
  if (g == NULL) {
    IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDARootInit", MSG_ROOT_FUNC_NULL);
    /*  return(IDA_ILL_INPUT); Masoud */
    return(IDA_NULL_G);
  }
  else gfun = g;

  /* Allocate necessary memory and return */
  glo = NULL;
  glo = (realtype *) malloc(nrt*sizeof(realtype));
  if (glo == NULL) {
    IDAProcessError(IDA_mem, IDA_MEM_FAIL, "IDA", "IDARootInit", MSG_MEM_FAIL);
    return(IDA_MEM_FAIL);
  }

  ghi = NULL;
  ghi = (realtype *) malloc(nrt*sizeof(realtype));
  if (ghi == NULL) {
    free(glo); glo = NULL;
    IDAProcessError(IDA_mem, IDA_MEM_FAIL, "IDA", "IDARootInit", MSG_MEM_FAIL);
    return(IDA_MEM_FAIL);
  }

  grout = NULL;
  grout = (realtype *) malloc(nrt*sizeof(realtype));
  if (grout == NULL) {
    free(glo); glo = NULL;
    free(ghi); ghi = NULL;
    IDAProcessError(IDA_mem, IDA_MEM_FAIL, "IDA", "IDARootInit", MSG_MEM_FAIL);
    return(IDA_MEM_FAIL);
  }

  iroots = NULL;
  iroots = (int *) malloc(nrt*sizeof(int));
  if (iroots == NULL) {
    free(glo); glo = NULL;
    free(ghi); ghi = NULL;
    free(grout); grout = NULL;
    IDAProcessError(IDA_mem, IDA_MEM_FAIL, "IDA", "IDARootInit", MSG_MEM_FAIL);
    return(IDA_MEM_FAIL);
  }

  lrw += 3*nrt;
  liw += nrt;

  return(IDA_SUCCESS);
}

/*
 * -----------------------------------------------------------------
 * readability constants
 * -----------------------------------------------------------------
 */

#define res (IDA_mem->ida_res)
#define y0  (IDA_mem->ida_y0)
#define yp0 (IDA_mem->ida_yp0)

#define itol  (IDA_mem->ida_itol)
#define rtol  (IDA_mem->ida_rtol)
#define Satol (IDA_mem->ida_Satol)
#define Vatol (IDA_mem->ida_Vatol)
#define efun  (IDA_mem->ida_efun)
#define edata (IDA_mem->ida_edata)

#define rdata       (IDA_mem->ida_rdata)
#define maxord      (IDA_mem->ida_maxord)
#define mxstep      (IDA_mem->ida_mxstep)
#define hin         (IDA_mem->ida_hin)
#define hmax_inv    (IDA_mem->ida_hmax_inv)
#define istop       (IDA_mem->ida_istop)
#define tstop       (IDA_mem->ida_tstop)
#define tstopset    (IDA_mem->ida_tstopset)
#define epcon       (IDA_mem->ida_epcon)
#define maxnef      (IDA_mem->ida_maxnef)
#define maxncf      (IDA_mem->ida_maxncf)
#define maxcor      (IDA_mem->ida_maxcor)
#define suppressalg (IDA_mem->ida_suppressalg)
#define id          (IDA_mem->ida_id)
#define constraints (IDA_mem->ida_constraints)

#define epiccon (IDA_mem->ida_epiccon)
#define maxnh   (IDA_mem->ida_maxnh)
#define maxnj   (IDA_mem->ida_maxnj)
#define maxnit  (IDA_mem->ida_maxnit)
#define lsoff   (IDA_mem->ida_lsoff)
#define steptol (IDA_mem->ida_steptol)

#define uround         (IDA_mem->ida_uround)  
#define phi            (IDA_mem->ida_phi) 
#define ewt            (IDA_mem->ida_ewt)  
#define yy             (IDA_mem->ida_yy)
#define yp             (IDA_mem->ida_yp)
#define delta          (IDA_mem->ida_delta)