Porting continues.
authorilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Mon, 26 Feb 2018 23:16:37 +0000 (23:16 +0000)
committerilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Mon, 26 Feb 2018 23:16:37 +0000 (23:16 +0000)
git-svn-id: https://citdcbmipav.cit.nih.gov/repos-pub/mipav/trunk@15391 ba61647d-9d00-f842-95cd-605cb4296b96

mipav/src/gov/nih/mipav/model/algorithms/CVODES.java

index 8b08cad..e3bad28 100644 (file)
@@ -94,6 +94,18 @@ public abstract class CVODES {
        final int CV_FUNCTIONAL = 1;\r
        final int CV_NEWTON = 2;\r
        \r
        final int CV_FUNCTIONAL = 1;\r
        final int CV_NEWTON = 2;\r
        \r
+       //itask: The itask input parameter to CVode indicates the job\r
+       //        of the solver for the next user step. The CV_NORMAL\r
+       //        itask is to have the solver take internal steps until\r
+       //        it has reached or just passed the user specified tout\r
+       //        parameter. The solver then interpolates in order to\r
+       //        return an approximate value of y(tout). The CV_ONE_STEP\r
+       //        option tells the solver to just take one internal step\r
+       //        and return the solution at the point reached by that step.\r
+       final int CV_NORMAL = 1;\r
+       final int CV_ONE_STEP = 2;\r
+\r
+       \r
        /*\r
         * Control constants for type of sensitivity RHS\r
         * ---------------------------------------------\r
        /*\r
         * Control constants for type of sensitivity RHS\r
         * ---------------------------------------------\r
@@ -327,6 +339,46 @@ public abstract class CVODES {
        final String MSGCV_BACK_ERROR  = "Error occured while integrating backward problem # %d"; \r
        final String MSGCV_BAD_TINTERP = "Bad t = %g for interpolation.";\r
        final String MSGCV_WRONG_INTERP = "This function cannot be called for the specified interp type.";\r
        final String MSGCV_BACK_ERROR  = "Error occured while integrating backward problem # %d"; \r
        final String MSGCV_BAD_TINTERP = "Bad t = %g for interpolation.";\r
        final String MSGCV_WRONG_INTERP = "This function cannot be called for the specified interp type.";\r
+       \r
+       final int CVDLS_SUCCESS = 0;\r
+       final int CVDLS_MEM_NULL = -1;\r
+       final int CVDLS_LMEM_NULL = -2;\r
+       final int CVDLS_ILL_INPUT = -3;\r
+       final int  CVDLS_MEM_FAIL = -4;\r
+\r
+       /* Additional last_flag values */\r
+\r
+       final int CVDLS_JACFUNC_UNRECVR = -5;\r
+       final int CVDLS_JACFUNC_RECVR = -6;\r
+       final int CVDLS_SUNMAT_FAIL = -7;\r
+\r
+       /* Return values for the adjoint module */\r
+\r
+       final int CVDLS_NO_ADJ = -101;\r
+       final int CVDLS_LMEMB_NULL = -102;\r
+\r
+       final String MSGD_CVMEM_NULL = "Integrator memory is NULL.";\r
+       final String MSGD_BAD_NVECTOR = "A required vector operation is not implemented.";\r
+       final String MSGD_BAD_SIZES = "Illegal bandwidth parameter(s). Must have 0 <=  ml, mu <= N-1.";\r
+       final String MSGD_MEM_FAIL = "A memory request failed.";\r
+       final String MSGD_LMEM_NULL ="Linear solver memory is NULL.";\r
+       final String MSGD_JACFUNC_FAILED = "The Jacobian routine failed in an unrecoverable manner.";\r
+       final String MSGD_MATCOPY_FAILED = "The SUNMatCopy routine failed in an unrecoverable manner.";\r
+       final String MSGD_MATZERO_FAILED = "The SUNMatZero routine failed in an unrecoverable manner.";\r
+       final String MSGD_MATSCALEADDI_FAILED = "The SUNMatScaleAddI routine failed in an unrecoverable manner.";\r
+\r
+       final String MSGD_NO_ADJ = "Illegal attempt to call before calling CVodeAdjMalloc.";\r
+       final String MSGD_BAD_WHICH = "Illegal value for which.";\r
+       final String MSGD_LMEMB_NULL = "Linear solver memory is NULL for the backward integration.";\r
+       final String MSGD_BAD_TINTERP = "Bad t for interpolation.";\r
+       \r
+       public enum SUNLinearSolver_Type{\r
+                 SUNLINEARSOLVER_DIRECT,\r
+                 SUNLINEARSOLVER_ITERATIVE,\r
+                 SUNLINEARSOLVER_CUSTOM\r
+               };\r
+\r
+    final int cvDlsDQJac = -1;\r
 \r
 \r
     final int cvsRoberts_dns = 1;\r
 \r
 \r
     final int cvsRoberts_dns = 1;\r
@@ -403,6 +455,7 @@ public abstract class CVODES {
                final int NOUT = 12; // number of output times\r
                int f = cvsRoberts_dns;\r
                int g = cvsRoberts_dns;\r
                final int NOUT = 12; // number of output times\r
                int f = cvsRoberts_dns;\r
                int g = cvsRoberts_dns;\r
+               int Jac = cvsRoberts_dns;\r
                \r
                // Set initial conditions\r
                NVector y = new NVector();\r
                \r
                // Set initial conditions\r
                NVector y = new NVector();\r
@@ -419,6 +472,9 @@ public abstract class CVODES {
                int flag;\r
                double A[][];\r
                SUNLinearSolver LS;\r
                int flag;\r
                double A[][];\r
                SUNLinearSolver LS;\r
+               double tout;\r
+               int iout;\r
+               double t[] = new double[1];\r
                \r
                // Call CVodeCreate to create the solver memory and specify the\r
                // Backward Differentiation Formula and the use of a Newton\r
                \r
                // Call CVodeCreate to create the solver memory and specify the\r
                // Backward Differentiation Formula and the use of a Newton\r
@@ -466,7 +522,27 @@ public abstract class CVODES {
                }\r
                \r
                // Call CVDlsSetLinearSolver to attach the matrix and linear solver to CVode\r
                }\r
                \r
                // Call CVDlsSetLinearSolver to attach the matrix and linear solver to CVode\r
-               //flag = CVDlsSetLinearSolver(cvode_mem, LS, A);\r
+               flag = CVDlsSetLinearSolver(cvode_mem, LS, A);\r
+               if (flag != CVDLS_SUCCESS) {\r
+                       return;\r
+               }\r
+               \r
+               // Set the user-supplied Jacobian routine Jac\r
+               flag = CVDlsSetJacFn(cvode_mem, Jac);\r
+               if (flag != CVDLS_SUCCESS) {\r
+                       return;\r
+               }\r
+               \r
+               // In loop, call CVode, print results, and test for error.\r
+               // Break out of loop when NOUT preset output times have been reached.\r
+               System.out.println(" \n3-species kinetics problem\n");\r
+               \r
+               iout = 0;\r
+               tout = T1;\r
+               \r
+               while (true) {\r
+                       flag = CVode(cvode_mem, tout, y, t, CV_NORMAL);\r
+               } // while (true)\r
        }\r
        \r
        private void N_VNew_Serial(NVector y, int length) {\r
        }\r
        \r
        private void N_VNew_Serial(NVector y, int length) {\r
@@ -510,6 +586,33 @@ public abstract class CVODES {
                return 0;\r
        }\r
        \r
                return 0;\r
        }\r
        \r
+       /**\r
+        * Jacobian routine.  Compute J(t,y) = df/dy.\r
+        * @param t\r
+        * @param y\r
+        * @param fy\r
+        * @param J\r
+        * @param tmp1\r
+        * @param tmp2\r
+        * @return\r
+        */\r
+       private int Jac(double t, NVector yv, NVector fy, double J[][], NVector tmp1, NVector tmp2) {\r
+               double y[] = yv.getData();\r
+           J[0][0] = -0.04;\r
+           J[0][1] = 1.0E4 * y[2];\r
+           J[0][2] = 1.0E4 * y[1];\r
+           \r
+           J[1][0] = 0.04;\r
+           J[1][1] = -1.0E4 * y[2] - 6.0E7*y[1];\r
+           J[1][2] = -1.0E4 * y[1];\r
+           \r
+           J[2][0] = ZERO;\r
+           J[2][1] = 6.0E7 * y[1];\r
+           J[2][2] = ZERO;\r
+           \r
+           return 0;\r
+       }\r
+       \r
        // Types: struct CVodeMemRec, CVodeMem\r
        // -----------------------------------------------------------------\r
        // The type CVodeMem is type pointer to struct CVodeMemRec.\r
        // Types: struct CVodeMemRec, CVodeMem\r
        // -----------------------------------------------------------------\r
        // The type CVodeMem is type pointer to struct CVodeMemRec.\r
@@ -539,6 +642,7 @@ public abstract class CVODES {
          //CVRhsFn cv_f;               /* y' = f(t,y(t))                                */\r
          int cv_f;\r
          //void *cv_user_data;         /* user pointer passed to f                      */\r
          //CVRhsFn cv_f;               /* y' = f(t,y(t))                                */\r
          int cv_f;\r
          //void *cv_user_data;         /* user pointer passed to f                      */\r
+         CVodeMemRec cv_user_data; \r
 \r
          int cv_lmm;                 /* lmm = ADAMS or BDF                            */\r
          int cv_iter;                /* iter = FUNCTIONAL or NEWTON                   */\r
 \r
          int cv_lmm;                 /* lmm = ADAMS or BDF                            */\r
          int cv_iter;                /* iter = FUNCTIONAL or NEWTON                   */\r
@@ -810,7 +914,8 @@ public abstract class CVODES {
 \r
          /* Linear Solver specific memory */\r
 \r
 \r
          /* Linear Solver specific memory */\r
 \r
-         //void *cv_lmem;           \r
+         //void *cv_lmem; \r
+         CVDlsMemRec cv_lmem;\r
 \r
          /* Flag to request a call to the setup routine */\r
 \r
 \r
          /* Flag to request a call to the setup routine */\r
 \r
@@ -1624,6 +1729,7 @@ public abstract class CVODES {
          long N; // Size of the linear system, the number of matrix rows\r
          long pivots[]; // Array of size N, index array for partial pivoting in LU factorization\r
          long last_flag; // last error return flag from internal setup/solve\r
          long N; // Size of the linear system, the number of matrix rows\r
          long pivots[]; // Array of size N, index array for partial pivoting in LU factorization\r
          long last_flag; // last error return flag from internal setup/solve\r
+         SUNLinearSolver_Type type;\r
      }\r
      \r
      private SUNLinearSolver SUNDenseLinearSolver(NVector y, double A[][]) {\r
      }\r
      \r
      private SUNLinearSolver SUNDenseLinearSolver(NVector y, double A[][]) {\r
@@ -1642,8 +1748,771 @@ public abstract class CVODES {
         S.N = matrixRows;\r
         S.last_flag = 0;\r
         S.pivots = new long[matrixRows];\r
         S.N = matrixRows;\r
         S.last_flag = 0;\r
         S.pivots = new long[matrixRows];\r
+        S.type = SUNLinearSolver_Type.SUNLINEARSOLVER_DIRECT;\r
         return S;\r
      }\r
         return S;\r
      }\r
+     \r
+     /*---------------------------------------------------------------\r
+     CVDlsSetLinearSolver specifies the direct linear solver.\r
+    ---------------------------------------------------------------*/\r
+    int CVDlsSetLinearSolver(CVodeMemRec cv_mem, SUNLinearSolver LS,\r
+                           double A[][])\r
+    {\r
+      CVDlsMemRec cvdls_mem;\r
+\r
+      /* Return immediately if any input is NULL */\r
+      if (cv_mem == null) {\r
+        cvProcessError(null, CVDLS_MEM_NULL, "CVSDLS", \r
+                       "CVDlsSetLinearSolver", MSGD_CVMEM_NULL);\r
+        return(CVDLS_MEM_NULL);\r
+      }\r
+      if ( (LS == null)  || (A == null) ) {\r
+        cvProcessError(null, CVDLS_ILL_INPUT, "CVSDLS", \r
+                       "CVDlsSetLinearSolver",\r
+                        "Both LS and A must be non-NULL");\r
+        return(CVDLS_ILL_INPUT);\r
+      }\r
+\r
+      /* Test if solver and vector are compatible with DLS */\r
+      if (LS.type != SUNLinearSolver_Type.SUNLINEARSOLVER_DIRECT) {\r
+        cvProcessError(cv_mem, CVDLS_ILL_INPUT, "CVSDLS", \r
+                       "CVDlsSetLinearSolver", \r
+                       "Non-direct LS supplied to CVDls interface");\r
+        return(CVDLS_ILL_INPUT);\r
+      }\r
+      //if (cv_mem.cv_tempv.ops.nvgetarraypointer == null ||\r
+          //cv_mem.cv_tempv.ops.nvsetarraypointer ==  null) {\r
+        //cvProcessError(cv_mem, CVDLS_ILL_INPUT, "CVSDLS", \r
+                       //"CVDlsSetLinearSolver", MSGD_BAD_NVECTOR);\r
+        //return(CVDLS_ILL_INPUT);\r
+      //}\r
+\r
+      /* free any existing system solver attached to CVode */\r
+      //if (cv_mem.cv_lfree)  cv_mem.cv_lfree(cv_mem);\r
+\r
+      /* Set four main system linear solver function fields in cv_mem */\r
+      //cv_mem.cv_linit  = cvDlsInitialize;\r
+      //cv_mem.cv_lsetup = cvDlsSetup;\r
+      //cv_mem.cv_lsolve = cvDlsSolve;\r
+      //cv_mem.cv_lfree  = cvDlsFree;\r
+      \r
+      /* Get memory for CVDlsMemRec */\r
+      cvdls_mem = null;\r
+      cvdls_mem = new CVDlsMemRec();\r
+      if (cvdls_mem == null) {\r
+        cvProcessError(cv_mem, CVDLS_MEM_FAIL, "CVSDLS", \r
+                        "CVDlsSetLinearSolver", MSGD_MEM_FAIL);\r
+        return(CVDLS_MEM_FAIL);\r
+      }\r
+\r
+      /* set SUNLinearSolver pointer */\r
+      cvdls_mem.LS = LS;\r
+      \r
+      /* Initialize Jacobian-related data */\r
+      cvdls_mem.jacDQ = true;\r
+      //cvdls_mem.jac = cvDlsDQJac;\r
+      cvdls_mem.J_data = cv_mem;\r
+      cvdls_mem.last_flag = CVDLS_SUCCESS;\r
+\r
+      /* Initialize counters */\r
+      cvDlsInitializeCounters(cvdls_mem);\r
+\r
+      /* Store pointer to A and create saved_J */\r
+      cvdls_mem.A = A;\r
+      cvdls_mem.savedJ = A.clone();\r
+      if (cvdls_mem.savedJ == null) {\r
+        cvProcessError(cv_mem, CVDLS_MEM_FAIL, "CVSDLS", \r
+                        "CVDlsSetLinearSolver", MSGD_MEM_FAIL);\r
+        cvdls_mem = null;\r
+        return(CVDLS_MEM_FAIL);\r
+      }\r
+\r
+      /* Allocate memory for x */\r
+      cvdls_mem.x = N_VClone(cv_mem.cv_tempv);\r
+      if (cvdls_mem.x == null) {\r
+        cvProcessError(cv_mem, CVDLS_MEM_FAIL, "CVSDLS", \r
+                        "CVDlsSetLinearSolver", MSGD_MEM_FAIL);\r
+        cvdls_mem.savedJ = null;\r
+        cvdls_mem = null;\r
+        return(CVDLS_MEM_FAIL);\r
+      }\r
+      /* Attach linear solver memory to integrator memory */\r
+      cv_mem.cv_lmem = cvdls_mem;\r
+\r
+      return(CVDLS_SUCCESS);\r
+    }\r
+    \r
+    //-----------------------------------------------------------------\r
+    //cvDlsInitializeCounters\r
+    //-----------------------------------------------------------------\r
+    //This routine resets the counters inside the CVDlsMem object.\r
+    //-----------------------------------------------------------------*/\r
+  int cvDlsInitializeCounters(CVDlsMemRec cvdls_mem)\r
+  {\r
+    cvdls_mem.nje   = 0;\r
+    cvdls_mem.nfeDQ = 0;\r
+    cvdls_mem.nstlj = 0;\r
+    return(0);\r
+  }\r
+\r
+    \r
+    \r
+    \r
+    //-----------------------------------------------------------------\r
+    //CVDlsMem is pointer to a CVDlsMemRec structure.\r
+    //-----------------------------------------------------------------*/\r
+\r
+  class CVDlsMemRec {\r
+\r
+    boolean jacDQ;    /* true if using internal DQ Jacobian approx. */\r
+    //CVDlsJacFn jac;       /* Jacobian routine to be called                 */\r
+    int jac;\r
+    //void *J_data;         /* data pointer passed to jac                    */\r
+    CVodeMemRec J_data;\r
+\r
+    double A[][];          /* A = I - gamma * df/dy                         */\r
+    double savedJ[][];     /* savedJ = old Jacobian                         */\r
+\r
+    SUNLinearSolver LS;   /* generic direct linear solver object           */\r
+\r
+    NVector x;           /* solution vector used by SUNLinearSolver       */\r
+    \r
+    long nstlj;       /* nstlj = nst at last Jacobian eval.            */\r
+\r
+    long nje;         /* nje = no. of calls to jac                     */\r
+\r
+    long nfeDQ;       /* no. of calls to f due to DQ Jacobian approx.  */\r
+\r
+    long last_flag;   /* last error return flag                        */\r
+\r
+  };\r
+\r
+\r
+  /* CVDlsSetJacFn specifies the Jacobian function. */\r
+  int CVDlsSetJacFn(CVodeMemRec cv_mem, int jac)\r
+  {\r
+    CVDlsMemRec cvdls_mem;\r
+\r
+    /* Return immediately if cvode_mem or cv_mem->cv_lmem are NULL */\r
+    if (cv_mem == null) {\r
+      cvProcessError(null, CVDLS_MEM_NULL, "CVSDLS",\r
+                     "CVDlsSetJacFn", MSGD_CVMEM_NULL);\r
+      return(CVDLS_MEM_NULL);\r
+    }\r
+  \r
+    if (cv_mem.cv_lmem == null) {\r
+      cvProcessError(cv_mem, CVDLS_LMEM_NULL, "CVSDLS",\r
+                     "CVDlsSetJacFn", MSGD_LMEM_NULL);\r
+      return(CVDLS_LMEM_NULL);\r
+    }\r
+    cvdls_mem = cv_mem.cv_lmem;\r
+\r
+    if (jac >= 0) {\r
+      cvdls_mem.jacDQ  = false;\r
+      cvdls_mem.jac    = jac;\r
+      cvdls_mem.J_data = cv_mem.cv_user_data;\r
+    } else {\r
+      cvdls_mem.jacDQ  = true;\r
+      cvdls_mem.jac    = cvDlsDQJac;\r
+      cvdls_mem.J_data = cv_mem;\r
+    }\r
+\r
+    return(CVDLS_SUCCESS);\r
+  }\r
+  \r
+  /*\r
+   * CVode\r
+   *\r
+   * This routine is the main driver of the CVODES package. \r
+   *\r
+   * It integrates over a time interval defined by the user, by calling\r
+   * cvStep to do internal time steps.\r
+   *\r
+   * The first time that CVode is called for a successfully initialized\r
+   * problem, it computes a tentative initial step size h.\r
+   *\r
+   * CVode supports two modes, specified by itask: CV_NORMAL, CV_ONE_STEP.\r
+   * In the CV_NORMAL mode, the solver steps until it reaches or passes tout\r
+   * and then interpolates to obtain y(tout).\r
+   * In the CV_ONE_STEP mode, it takes one internal step and returns.\r
+   */\r
+\r
+  int CVode(CVodeMemRec cv_mem, double tout, NVector yout, \r
+            double tret[], int itask)\r
+  {\r
+    long nstloc; \r
+    int retval, hflag, kflag, istate, is, ir, ier, irfndp;\r
+    double troundoff, tout_hin, rh, nrm;\r
+    boolean inactive_roots;\r
+\r
+    /*\r
+     * -------------------------------------\r
+     * 1. Check and process inputs\r
+     * -------------------------------------\r
+     */\r
+\r
+    /* Check if cvode_mem exists */\r
+    if (cv_mem == null) {\r
+      cvProcessError(null, CV_MEM_NULL, "CVODES", "CVode",\r
+                     MSGCV_NO_MEM);\r
+      return(CV_MEM_NULL);\r
+    }\r
+\r
+    /* Check if cvode_mem was allocated */\r
+    if (cv_mem.cv_MallocDone == false) {\r
+      cvProcessError(cv_mem, CV_NO_MALLOC, "CVODES", "CVode",\r
+                     MSGCV_NO_MALLOC);\r
+      return(CV_NO_MALLOC);\r
+    }\r
+    \r
+    /* Check for yout != NULL */\r
+    if ((cv_mem.cv_y = yout) == null) {\r
+      cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode",\r
+                     MSGCV_YOUT_NULL);\r
+      return(CV_ILL_INPUT);\r
+    }\r
+    \r
+    /* Check for tret != NULL */\r
+    if (tret == null) {\r
+      cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode",\r
+                     MSGCV_TRET_NULL);\r
+      return(CV_ILL_INPUT);\r
+    }\r
+\r
+    /* Check for valid itask */\r
+    if ( (itask != CV_NORMAL) && (itask != CV_ONE_STEP) ) {\r
+      cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode",\r
+                     MSGCV_BAD_ITASK);\r
+      return(CV_ILL_INPUT);\r
+    }\r
+\r
+    if (itask == CV_NORMAL) cv_mem.cv_toutc = tout;\r
+    cv_mem.cv_taskc = itask;\r
+\r
+    /*\r
+     * ----------------------------------------\r
+     * 2. Initializations performed only at\r
+     *    the first step (nst=0):\r
+     *    - initial setup\r
+     *    - initialize Nordsieck history array\r
+     *    - compute initial step size\r
+     *    - check for approach to tstop\r
+     *    - check for approach to a root\r
+     * ----------------------------------------\r
+     */\r
+\r
+   /* if (cv_mem.cv_nst == 0) {\r
+\r
+      cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tn;*/\r
+\r
+      /* Check inputs for corectness */\r
+\r
+      /*ier = cvInitialSetup(cv_mem);\r
+      if (ier!= CV_SUCCESS) return(ier);*/\r
+\r
+      /* \r
+       * Call f at (t0,y0), set zn[1] = y'(t0). \r
+       * If computing any quadratures, call fQ at (t0,y0), set znQ[1] = yQ'(t0)\r
+       * If computing sensitivities, call fS at (t0,y0,yS0), set znS[1][is] = yS'(t0), is=1,...,Ns.\r
+       * If computing quadr. sensi., call fQS at (t0,y0,yS0), set znQS[1][is] = yQS'(t0), is=1,...,Ns.\r
+       */\r
+\r
+      /*retval = cv_mem->cv_f(cv_mem->cv_tn, cv_mem->cv_zn[0],\r
+                            cv_mem->cv_zn[1], cv_mem->cv_user_data); \r
+      cv_mem->cv_nfe++;\r
+      if (retval < 0) {\r
+        cvProcessError(cv_mem, CV_RHSFUNC_FAIL, "CVODES", "CVode",\r
+                       MSGCV_RHSFUNC_FAILED, cv_mem->cv_tn);\r
+        return(CV_RHSFUNC_FAIL);\r
+      }\r
+      if (retval > 0) {\r
+        cvProcessError(cv_mem, CV_FIRST_RHSFUNC_ERR, "CVODES", "CVode",\r
+                       MSGCV_RHSFUNC_FIRST);\r
+        return(CV_FIRST_RHSFUNC_ERR);\r
+      }\r
+\r
+      if (cv_mem->cv_quadr) {\r
+        retval = cv_mem->cv_fQ(cv_mem->cv_tn, cv_mem->cv_zn[0],\r
+                               cv_mem->cv_znQ[1], cv_mem->cv_user_data);\r
+        cv_mem->cv_nfQe++;\r
+        if (retval < 0) {\r
+          cvProcessError(cv_mem, CV_QRHSFUNC_FAIL, "CVODES", "CVode",\r
+                         MSGCV_QRHSFUNC_FAILED, cv_mem->cv_tn);\r
+          return(CV_QRHSFUNC_FAIL);\r
+        }\r
+        if (retval > 0) {\r
+          cvProcessError(cv_mem, CV_FIRST_QRHSFUNC_ERR, "CVODES",\r
+                         "CVode", MSGCV_QRHSFUNC_FIRST);\r
+          return(CV_FIRST_QRHSFUNC_ERR);\r
+        }\r
+      }\r
+\r
+      if (cv_mem->cv_sensi) {\r
+        retval = cvSensRhsWrapper(cv_mem, cv_mem->cv_tn, cv_mem->cv_zn[0],\r
+                                  cv_mem->cv_zn[1], cv_mem->cv_znS[0],\r
+                                  cv_mem->cv_znS[1], cv_mem->cv_tempv,\r
+                                  cv_mem->cv_ftemp);\r
+        if (retval < 0) {\r
+          cvProcessError(cv_mem, CV_SRHSFUNC_FAIL, "CVODES", "CVode",\r
+                         MSGCV_SRHSFUNC_FAILED, cv_mem->cv_tn);\r
+          return(CV_SRHSFUNC_FAIL);\r
+        } \r
+        if (retval > 0) {\r
+          cvProcessError(cv_mem, CV_FIRST_SRHSFUNC_ERR, "CVODES",\r
+                         "CVode", MSGCV_SRHSFUNC_FIRST);\r
+          return(CV_FIRST_SRHSFUNC_ERR);\r
+        }\r
+      }\r
+\r
+      if (cv_mem->cv_quadr_sensi) {\r
+        retval = cv_mem->cv_fQS(cv_mem->cv_Ns, cv_mem->cv_tn, cv_mem->cv_zn[0],\r
+                                cv_mem->cv_znS[0], cv_mem->cv_znQ[1],\r
+                                cv_mem->cv_znQS[1], cv_mem->cv_fQS_data,\r
+                                cv_mem->cv_tempv, cv_mem->cv_tempvQ); \r
+        cv_mem->cv_nfQSe++;\r
+        if (retval < 0) {\r
+          cvProcessError(cv_mem, CV_QSRHSFUNC_FAIL, "CVODES", "CVode",\r
+                         MSGCV_QSRHSFUNC_FAILED, cv_mem->cv_tn);\r
+          return(CV_QSRHSFUNC_FAIL);\r
+        } \r
+        if (retval > 0) {\r
+          cvProcessError(cv_mem, CV_FIRST_QSRHSFUNC_ERR, "CVODES",\r
+                         "CVode", MSGCV_QSRHSFUNC_FIRST);\r
+          return(CV_FIRST_QSRHSFUNC_ERR);\r
+        }\r
+      }*/\r
+\r
+      /* Test input tstop for legality. */\r
+\r
+     /* if (cv_mem->cv_tstopset) {\r
+        if ( (cv_mem->cv_tstop - cv_mem->cv_tn)*(tout - cv_mem->cv_tn) <= ZERO ) {\r
+          cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode",\r
+                         MSGCV_BAD_TSTOP, cv_mem->cv_tstop, cv_mem->cv_tn);\r
+          return(CV_ILL_INPUT);\r
+        }\r
+      }*/\r
+\r
+      /* Set initial h (from H0 or cvHin). */\r
+      \r
+      /*cv_mem->cv_h = cv_mem->cv_hin;\r
+      if ( (cv_mem->cv_h != ZERO) && ((tout-cv_mem->cv_tn)*cv_mem->cv_h < ZERO) ) {\r
+        cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode", MSGCV_BAD_H0);\r
+        return(CV_ILL_INPUT);\r
+      }\r
+      if (cv_mem->cv_h == ZERO) {\r
+        tout_hin = tout;\r
+        if ( cv_mem->cv_tstopset &&\r
+             (tout-cv_mem->cv_tn)*(tout-cv_mem->cv_tstop) > ZERO )\r
+          tout_hin = cv_mem->cv_tstop; \r
+        hflag = cvHin(cv_mem, tout_hin);\r
+        if (hflag != CV_SUCCESS) {\r
+          istate = cvHandleFailure(cv_mem, hflag);\r
+          return(istate);\r
+        }\r
+      }\r
+      rh = SUNRabs(cv_mem->cv_h)*cv_mem->cv_hmax_inv;\r
+      if (rh > ONE) cv_mem->cv_h /= rh;\r
+      if (SUNRabs(cv_mem->cv_h) < cv_mem->cv_hmin)\r
+        cv_mem->cv_h *= cv_mem->cv_hmin/SUNRabs(cv_mem->cv_h);*/\r
+\r
+      /* Check for approach to tstop */\r
+\r
+      /*if (cv_mem->cv_tstopset) {\r
+        if ( (cv_mem->cv_tn + cv_mem->cv_h - cv_mem->cv_tstop)*cv_mem->cv_h > ZERO ) \r
+          cv_mem->cv_h = (cv_mem->cv_tstop - cv_mem->cv_tn)*(ONE-FOUR*cv_mem->cv_uround);\r
+      }*/\r
+\r
+      /* \r
+       * Scale zn[1] by h.\r
+       * If computing any quadratures, scale znQ[1] by h.\r
+       * If computing sensitivities,  scale znS[1][is] by h. \r
+       * If computing quadrature sensitivities,  scale znQS[1][is] by h. \r
+       */\r
+\r
+      /*cv_mem->cv_hscale = cv_mem->cv_h;\r
+      cv_mem->cv_h0u    = cv_mem->cv_h;\r
+      cv_mem->cv_hprime = cv_mem->cv_h;\r
+\r
+      N_VScale(cv_mem->cv_h, cv_mem->cv_zn[1], cv_mem->cv_zn[1]);\r
+      \r
+      if (cv_mem->cv_quadr)\r
+        N_VScale(cv_mem->cv_h, cv_mem->cv_znQ[1], cv_mem->cv_znQ[1]);\r
+\r
+      if (cv_mem->cv_sensi)\r
+        for (is=0; is<cv_mem->cv_Ns; is++) \r
+          N_VScale(cv_mem->cv_h, cv_mem->cv_znS[1][is], cv_mem->cv_znS[1][is]);\r
+\r
+      if (cv_mem->cv_quadr_sensi)\r
+        for (is=0; is<cv_mem->cv_Ns; is++) \r
+          N_VScale(cv_mem->cv_h, cv_mem->cv_znQS[1][is], cv_mem->cv_znQS[1][is]);*/\r
+      \r
+      /* Check for zeros of root function g at and near t0. */\r
+\r
+      /*if (cv_mem->cv_nrtfn > 0) {\r
+\r
+        retval = cvRcheck1(cv_mem);\r
+\r
+        if (retval == CV_RTFUNC_FAIL) {\r
+          cvProcessError(cv_mem, CV_RTFUNC_FAIL, "CVODES", "cvRcheck1",\r
+                         MSGCV_RTFUNC_FAILED, cv_mem->cv_tn);\r
+          return(CV_RTFUNC_FAIL);\r
+        }\r
+\r
+      }\r
+\r
+    } *//* end first call block */\r
+\r
+    /*\r
+     * ------------------------------------------------------\r
+     * 3. At following steps, perform stop tests:\r
+     *    - check for root in last step\r
+     *    - check if we passed tstop\r
+     *    - check if we passed tout (NORMAL mode)\r
+     *    - check if current tn was returned (ONE_STEP mode)\r
+     *    - check if we are close to tstop\r
+     *      (adjust step size if needed)\r
+     * -------------------------------------------------------\r
+     */\r
+\r
+   /* if (cv_mem->cv_nst > 0) {*/\r
+\r
+      /* Estimate an infinitesimal time interval to be used as\r
+         a roundoff for time quantities (based on current time \r
+         and step size) */\r
+      /*troundoff = FUZZ_FACTOR * cv_mem->cv_uround *\r
+        (SUNRabs(cv_mem->cv_tn) + SUNRabs(cv_mem->cv_h));*/\r
+\r
+      /* First check for a root in the last step taken, other than the\r
+         last root found, if any.  If itask = CV_ONE_STEP and y(tn) was not\r
+         returned because of an intervening root, return y(tn) now.     */\r
+      /*if (cv_mem->cv_nrtfn > 0) {\r
+        \r
+        irfndp = cv_mem->cv_irfnd;\r
+        \r
+        retval = cvRcheck2(cv_mem);\r
+\r
+        if (retval == CLOSERT) {\r
+          cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "cvRcheck2",\r
+                         MSGCV_CLOSE_ROOTS, cv_mem->cv_tlo);\r
+          return(CV_ILL_INPUT);\r
+        } else if (retval == CV_RTFUNC_FAIL) {\r
+          cvProcessError(cv_mem, CV_RTFUNC_FAIL, "CVODES", "cvRcheck2",\r
+                         MSGCV_RTFUNC_FAILED, cv_mem->cv_tlo);\r
+          return(CV_RTFUNC_FAIL);\r
+        } else if (retval == RTFOUND) {\r
+          cv_mem->cv_tretlast = *tret = cv_mem->cv_tlo;\r
+          return(CV_ROOT_RETURN);\r
+        }*/\r
+        \r
+        /* If tn is distinct from tretlast (within roundoff),\r
+           check remaining interval for roots */\r
+        /*if ( SUNRabs(cv_mem->cv_tn - cv_mem->cv_tretlast) > troundoff ) {\r
+\r
+          retval = cvRcheck3(cv_mem);\r
+\r
+          if (retval == CV_SUCCESS) {*/     /* no root found */\r
+            /*cv_mem->cv_irfnd = 0;\r
+            if ((irfndp == 1) && (itask == CV_ONE_STEP)) {\r
+              cv_mem->cv_tretlast = *tret = cv_mem->cv_tn;\r
+              N_VScale(ONE, cv_mem->cv_zn[0], yout);\r
+              return(CV_SUCCESS);\r
+            }\r
+          } else if (retval == RTFOUND) {*/  /* a new root was found */\r
+            /*cv_mem->cv_irfnd = 1;\r
+            cv_mem->cv_tretlast = *tret = cv_mem->cv_tlo;\r
+            return(CV_ROOT_RETURN);\r
+          } else if (retval == CV_RTFUNC_FAIL) { */ /* g failed */\r
+            /*cvProcessError(cv_mem, CV_RTFUNC_FAIL, "CVODES", "cvRcheck3", \r
+                           MSGCV_RTFUNC_FAILED, cv_mem->cv_tlo);\r
+            return(CV_RTFUNC_FAIL);\r
+          }\r
+\r
+        }\r
+        \r
+      } *//* end of root stop check */\r
+      \r
+      /* In CV_NORMAL mode, test if tout was reached */\r
+     /* if ( (itask == CV_NORMAL) && ((cv_mem->cv_tn-tout)*cv_mem->cv_h >= ZERO) ) {\r
+        cv_mem->cv_tretlast = *tret = tout;\r
+        ier =  CVodeGetDky(cv_mem, tout, 0, yout);\r
+        if (ier != CV_SUCCESS) {\r
+          cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode",\r
+                         MSGCV_BAD_TOUT, tout);\r
+          return(CV_ILL_INPUT);\r
+        }\r
+        return(CV_SUCCESS);\r
+      }*/\r
+      \r
+      /* In CV_ONE_STEP mode, test if tn was returned */\r
+     /* if ( itask == CV_ONE_STEP &&\r
+           SUNRabs(cv_mem->cv_tn - cv_mem->cv_tretlast) > troundoff ) {\r
+        cv_mem->cv_tretlast = *tret = cv_mem->cv_tn;\r
+        N_VScale(ONE, cv_mem->cv_zn[0], yout);\r
+        return(CV_SUCCESS);\r
+      }\r
+      \r
+      /* Test for tn at tstop or near tstop */\r
+     /* if ( cv_mem->cv_tstopset ) {\r
+        \r
+        if ( SUNRabs(cv_mem->cv_tn - cv_mem->cv_tstop) <= troundoff ) {\r
+          ier =  CVodeGetDky(cv_mem, cv_mem->cv_tstop, 0, yout);\r
+          if (ier != CV_SUCCESS) {\r
+            cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode",\r
+                           MSGCV_BAD_TSTOP, cv_mem->cv_tstop, cv_mem->cv_tn);\r
+            return(CV_ILL_INPUT);\r
+          }\r
+          cv_mem->cv_tretlast = *tret = cv_mem->cv_tstop;\r
+          cv_mem->cv_tstopset = SUNFALSE;\r
+          return(CV_TSTOP_RETURN);\r
+        }*/\r
+        \r
+        /* If next step would overtake tstop, adjust stepsize */\r
+        /*if ( (cv_mem->cv_tn + cv_mem->cv_hprime - cv_mem->cv_tstop)*cv_mem->cv_h > ZERO ) {\r
+          cv_mem->cv_hprime = (cv_mem->cv_tstop - cv_mem->cv_tn)*(ONE-FOUR*cv_mem->cv_uround);\r
+          cv_mem->cv_eta = cv_mem->cv_hprime / cv_mem->cv_h;\r
+        }\r
+        \r
+      }\r
+      \r
+    }*/ /* end stopping tests block at nst>0 */  \r
+    \r
+    /*\r
+     * --------------------------------------------------\r
+     * 4. Looping point for internal steps\r
+     *\r
+     *    4.1. check for errors (too many steps, too much\r
+     *         accuracy requested, step size too small)\r
+     *    4.2. take a new step (call cvStep)\r
+     *    4.3. stop on error \r
+     *    4.4. perform stop tests:\r
+     *         - check for root in last step\r
+     *         - check if tout was passed\r
+     *         - check if close to tstop\r
+     *         - check if in ONE_STEP mode (must return)\r
+     * --------------------------------------------------\r
+     */  \r
+    \r
+    /*nstloc = 0;\r
+    for(;;) {\r
+     \r
+      cv_mem->cv_next_h = cv_mem->cv_h;\r
+      cv_mem->cv_next_q = cv_mem->cv_q;*/\r
+      \r
+      /* Reset and check ewt, ewtQ, ewtS */   \r
+      /*if (cv_mem->cv_nst > 0) {\r
+\r
+        ier = cv_mem->cv_efun(cv_mem->cv_zn[0], cv_mem->cv_ewt, cv_mem->cv_e_data);\r
+        if(ier != 0) {\r
+          if (cv_mem->cv_itol == CV_WF)\r
+            cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode",\r
+                           MSGCV_EWT_NOW_FAIL, cv_mem->cv_tn);\r
+          else\r
+            cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode",\r
+                           MSGCV_EWT_NOW_BAD, cv_mem->cv_tn);\r
+          istate = CV_ILL_INPUT;\r
+          cv_mem->cv_tretlast = *tret = cv_mem->cv_tn;\r
+          N_VScale(ONE, cv_mem->cv_zn[0], yout);\r
+          break;\r
+        }\r
+\r
+        if (cv_mem->cv_quadr && cv_mem->cv_errconQ) {\r
+          ier = cvQuadEwtSet(cv_mem, cv_mem->cv_znQ[0], cv_mem->cv_ewtQ);\r
+          if(ier != 0) {\r
+            cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode",\r
+                           MSGCV_EWTQ_NOW_BAD, cv_mem->cv_tn);\r
+            istate = CV_ILL_INPUT;\r
+            cv_mem->cv_tretlast = *tret = cv_mem->cv_tn;\r
+            N_VScale(ONE, cv_mem->cv_zn[0], yout);\r
+            break;\r
+          }\r
+        }\r
+\r
+        if (cv_mem->cv_sensi) {\r
+          ier = cvSensEwtSet(cv_mem, cv_mem->cv_znS[0], cv_mem->cv_ewtS);\r
+          if (ier != 0) {\r
+            cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode",\r
+                           MSGCV_EWTS_NOW_BAD, cv_mem->cv_tn);\r
+            istate = CV_ILL_INPUT;\r
+            cv_mem->cv_tretlast = *tret = cv_mem->cv_tn;\r
+            N_VScale(ONE, cv_mem->cv_zn[0], yout);\r
+            break;\r
+          }\r
+        }\r
+\r
+        if (cv_mem->cv_quadr_sensi && cv_mem->cv_errconQS) {\r
+          ier = cvQuadSensEwtSet(cv_mem, cv_mem->cv_znQS[0], cv_mem->cv_ewtQS);\r
+          if (ier != 0) {\r
+            cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode",\r
+                           MSGCV_EWTQS_NOW_BAD, cv_mem->cv_tn);\r
+            istate = CV_ILL_INPUT;\r
+            cv_mem->cv_tretlast = *tret = cv_mem->cv_tn;\r
+            N_VScale(ONE, cv_mem->cv_zn[0], yout);\r
+            break;\r
+          }\r
+        }\r
+\r
+      }*/\r
+\r
+      /* Check for too many steps */\r
+      /*if ( (cv_mem->cv_mxstep>0) && (nstloc >= cv_mem->cv_mxstep) ) {\r
+        cvProcessError(cv_mem, CV_TOO_MUCH_WORK, "CVODES", "CVode",\r
+                       MSGCV_MAX_STEPS, cv_mem->cv_tn);\r
+        istate = CV_TOO_MUCH_WORK;\r
+        cv_mem->cv_tretlast = *tret = cv_mem->cv_tn;\r
+        N_VScale(ONE, cv_mem->cv_zn[0], yout);\r
+        break;\r
+      }*/\r
+\r
+      /* Check for too much accuracy requested */\r
+      /*nrm = N_VWrmsNorm(cv_mem->cv_zn[0], cv_mem->cv_ewt);\r
+      if (cv_mem->cv_quadr && cv_mem->cv_errconQ) {\r
+        nrm = cvQuadUpdateNorm(cv_mem, nrm, cv_mem->cv_znQ[0], cv_mem->cv_ewtQ); \r
+      }\r
+      if (cv_mem->cv_sensi && cv_mem->cv_errconS) {\r
+        nrm = cvSensUpdateNorm(cv_mem, nrm, cv_mem->cv_znS[0], cv_mem->cv_ewtS);\r
+      }\r
+      if (cv_mem->cv_quadr_sensi && cv_mem->cv_errconQS) {\r
+        nrm = cvQuadSensUpdateNorm(cv_mem, nrm, cv_mem->cv_znQS[0], cv_mem->cv_ewtQS);\r
+      }\r
+      cv_mem->cv_tolsf = cv_mem->cv_uround * nrm;\r
+      if (cv_mem->cv_tolsf > ONE) {\r
+        cvProcessError(cv_mem, CV_TOO_MUCH_ACC, "CVODES", "CVode",\r
+                       MSGCV_TOO_MUCH_ACC, cv_mem->cv_tn);\r
+        istate = CV_TOO_MUCH_ACC;\r
+        cv_mem->cv_tretlast = *tret = cv_mem->cv_tn;\r
+        N_VScale(ONE, cv_mem->cv_zn[0], yout);\r
+        cv_mem->cv_tolsf *= TWO;\r
+        break;\r
+      } else {\r
+        cv_mem->cv_tolsf = ONE;\r
+      }*/\r
+      \r
+      /* Check for h below roundoff level in tn */\r
+     /* if (cv_mem->cv_tn + cv_mem->cv_h == cv_mem->cv_tn) {\r
+        cv_mem->cv_nhnil++;\r
+        if (cv_mem->cv_nhnil <= cv_mem->cv_mxhnil) \r
+          cvProcessError(cv_mem, CV_WARNING, "CVODES", "CVode", MSGCV_HNIL,\r
+                         cv_mem->cv_tn, cv_mem->cv_h);\r
+        if (cv_mem->cv_nhnil == cv_mem->cv_mxhnil) \r
+          cvProcessError(cv_mem, CV_WARNING, "CVODES", "CVode", MSGCV_HNIL_DONE);\r
+      }*/\r
+\r
+      /* Call cvStep to take a step */\r
+      /*kflag = cvStep(cv_mem); */\r
+\r
+      /* Process failed step cases, and exit loop */\r
+     /* if (kflag != CV_SUCCESS) {\r
+        istate = cvHandleFailure(cv_mem, kflag);\r
+        cv_mem->cv_tretlast = *tret = cv_mem->cv_tn;\r
+        N_VScale(ONE, cv_mem->cv_zn[0], yout);\r
+        break;\r
+      }\r
+      \r
+      nstloc++;*/\r
+\r
+      /* If tstop is set and was reached, reset tn = tstop */\r
+      /*if ( cv_mem->cv_tstopset ) {\r
+        troundoff = FUZZ_FACTOR * cv_mem->cv_uround *\r
+          (SUNRabs(cv_mem->cv_tn) + SUNRabs(cv_mem->cv_h));\r
+        if ( SUNRabs(cv_mem->cv_tn - cv_mem->cv_tstop) <= troundoff)\r
+          cv_mem->cv_tn = cv_mem->cv_tstop;\r
+      }*/\r
+\r
+      /* Check for root in last step taken. */    \r
+      /*if (cv_mem->cv_nrtfn > 0) {\r
+        \r
+        retval = cvRcheck3(cv_mem);\r
+        \r
+        if (retval == RTFOUND) { */ /* A new root was found */\r
+         /* cv_mem->cv_irfnd = 1;\r
+          istate = CV_ROOT_RETURN;\r
+          cv_mem->cv_tretlast = *tret = cv_mem->cv_tlo;\r
+          break;\r
+        } else if (retval == CV_RTFUNC_FAIL) {*/ /* g failed */\r
+         /* cvProcessError(cv_mem, CV_RTFUNC_FAIL, "CVODES", "cvRcheck3",\r
+                         MSGCV_RTFUNC_FAILED, cv_mem->cv_tlo);\r
+          istate = CV_RTFUNC_FAIL;\r
+          break;\r
+        }*/\r
+\r
+        /* If we are at the end of the first step and we still have\r
+         * some event functions that are inactive, issue a warning\r
+         * as this may indicate a user error in the implementation\r
+         * of the root function. */\r
+\r
+       /* if (cv_mem->cv_nst==1) {\r
+          inactive_roots = SUNFALSE;\r
+          for (ir=0; ir<cv_mem->cv_nrtfn; ir++) { \r
+            if (!cv_mem->cv_gactive[ir]) {\r
+              inactive_roots = SUNTRUE;\r
+              break;\r
+            }\r
+          }\r
+          if ((cv_mem->cv_mxgnull > 0) && inactive_roots) {\r
+            cvProcessError(cv_mem, CV_WARNING, "CVODES", "CVode",\r
+                           MSGCV_INACTIVE_ROOTS);\r
+          }\r
+        }\r
+\r
+      }*/\r
+\r
+      /* In NORMAL mode, check if tout reached */\r
+     /* if ( (itask == CV_NORMAL) &&  (cv_mem->cv_tn-tout)*cv_mem->cv_h >= ZERO ) {\r
+        istate = CV_SUCCESS;\r
+        cv_mem->cv_tretlast = *tret = tout;\r
+        (void) CVodeGetDky(cv_mem, tout, 0, yout);\r
+        cv_mem->cv_next_q = cv_mem->cv_qprime;\r
+        cv_mem->cv_next_h = cv_mem->cv_hprime;\r
+        break;\r
+      }*/\r
+\r
+      /* Check if tn is at tstop, or about to pass tstop */\r
+     /* if ( cv_mem->cv_tstopset ) {\r
+\r
+        troundoff = FUZZ_FACTOR * cv_mem->cv_uround *\r
+          (SUNRabs(cv_mem->cv_tn) + SUNRabs(cv_mem->cv_h));\r
+        if ( SUNRabs(cv_mem->cv_tn - cv_mem->cv_tstop) <= troundoff) {\r
+          (void) CVodeGetDky(cv_mem, cv_mem->cv_tstop, 0, yout);\r
+          cv_mem->cv_tretlast = *tret = cv_mem->cv_tstop;\r
+          cv_mem->cv_tstopset = SUNFALSE;\r
+          istate = CV_TSTOP_RETURN;\r
+          break;\r
+        }\r
+\r
+        if ( (cv_mem->cv_tn + cv_mem->cv_hprime - cv_mem->cv_tstop)*cv_mem->cv_h > ZERO ) {\r
+          cv_mem->cv_hprime = (cv_mem->cv_tstop - cv_mem->cv_tn)*(ONE-FOUR*cv_mem->cv_uround);\r
+          cv_mem->cv_eta = cv_mem->cv_hprime / cv_mem->cv_h;\r
+        }\r
+\r
+      }*/\r
+\r
+      /* In ONE_STEP mode, copy y and exit loop */\r
+      /*if (itask == CV_ONE_STEP) {\r
+        istate = CV_SUCCESS;\r
+        cv_mem->cv_tretlast = *tret = cv_mem->cv_tn;\r
+        N_VScale(ONE, cv_mem->cv_zn[0], yout);\r
+        cv_mem->cv_next_q = cv_mem->cv_qprime;\r
+        cv_mem->cv_next_h = cv_mem->cv_hprime;\r
+        break;\r
+      }\r
+\r
+    } *//* end looping for internal steps */\r
+    \r
+    /* Load optional output */\r
+   /* if (cv_mem->cv_sensi && (cv_mem->cv_ism==CV_STAGGERED1)) { \r
+      cv_mem->cv_nniS  = 0;\r
+      cv_mem->cv_ncfnS = 0;\r
+      for (is=0; is<cv_mem->cv_Ns; is++) {\r
+        cv_mem->cv_nniS  += cv_mem->cv_nniS1[is];\r
+        cv_mem->cv_ncfnS += cv_mem->cv_ncfnS1[is];\r
+      }\r
+    }*/\r
+    \r
+    /*return(istate);*/\r
+    return -100;\r
+\r
+  }\r
+\r
 \r
 \r
 }
\ No newline at end of file
 \r
 \r
 }
\ No newline at end of file