Porting continues.
authorilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Tue, 3 Apr 2018 22:15:13 +0000 (22:15 +0000)
committerilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Tue, 3 Apr 2018 22:15:13 +0000 (22:15 +0000)
git-svn-id: https://citdcbmipav.cit.nih.gov/repos-pub/mipav/trunk@15438 ba61647d-9d00-f842-95cd-605cb4296b96

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

index 0ef94ce..babd29a 100644 (file)
@@ -4,6 +4,7 @@ package gov.nih.mipav.model.algorithms;
 import gov.nih.mipav.view.MipavUtil;\r
 \r
 public abstract class CVODES_ASA extends CVODES {\r
+       final int CVArhs_select = 100;\r
        public CVODES_ASA() {\r
                super();\r
                runcvsRoberts_ASAi_dns();\r
@@ -147,6 +148,14 @@ public abstract class CVODES_ASA extends CVODES {
                int steps;\r
                double time[] = new double[1];\r
                int ncheck[] = new int[1];\r
+               long nst;\r
+               int i;\r
+               CVadjCheckPointRec ckpnt[];\r
+               double reltolB;\r
+               double abstolB;\r
+               double abstolQB;\r
+               int indexB[] = new int[1];\r
+               int fB = cvsRoberts_ASAi_dns;\r
                \r
                // Print problem description \r
                System.out.printf("\nAdjoint Sensitivity Example for Chemical Kinetics\n");\r
@@ -289,6 +298,66 @@ public abstract class CVODES_ASA extends CVODES {
         if (flag < 0) {\r
                return;\r
         }\r
+        // Get the current number of integration steps\r
+        nst = cvode_mem.cv_nst;\r
+\r
+        System.out.printf("done ( nst = " + nst+")");\r
+        System.out.printf("\nncheck = " + ncheck[0] + "\n\n");\r
+        \r
+        flag = CVodeGetQuad(cvode_mem, time, q);\r
+        \r
+        System.out.printf("------------------------------------------------\n");\r
+        System.out.printf("G: " + q.data[0] + "\n");\r
+        System.out.printf("------------------------------------------------\n\n");\r
+        \r
+        // Test check point linked list\r
+        // (Uncomment next block to print check point information)\r
+        /*System.out.printf("\nList of Check Points (ncheck = " + ncheck[0] + ")\n\n");\r
+        ckpnt = new CVadjCheckPointRec[ncheck[0]+1];\r
+        for (i = 0; i < ncheck[0]+1; i++) {\r
+               ckpnt[i] = new CVadjCheckPointRec();\r
+        }\r
+        CVodeGetAdjCheckPointsInfo(cvode_mem, ckpnt);\r
+        for (i = 0; i <= ncheck[0]; i++) {\r
+               System.out.printf("Time interval: " + ckpnt[i].t0 + ", " +  ckpnt[i].t1 + "\n");\r
+               System.out.printf("Order: " + ckpnt[i].nstep + "\n");\r
+               System.out.printf("Step size: " + ckpnt[i].step + "\n");\r
+               System.out.printf("\n");\r
+        }*/\r
+        \r
+        // Initialize yB\r
+        NVector yB = new NVector();\r
+               double yBr[] = new double[]{0.0,0.0,0.0};\r
+               N_VNew_Serial(yB, NEQ);\r
+               yB.setData(yBr);\r
+               \r
+               // Initialize qB\r
+               NVector qB = new NVector();\r
+               double qBr[] = new double[]{0.0,0.0,0.0};\r
+               N_VNew_Serial(qB,NP);\r
+               qB.setData(qBr);\r
+               \r
+               // Set the scalar relative tolerance reltolB\r
+               reltolB = RTOL;\r
+               \r
+               // Set the scalar absolute tolerance abstolB\r
+               abstolB = ATOLl;\r
+               \r
+               // Set the scalar absolute tolerance abstolQB\r
+               abstolQB = ATOLq;\r
+               \r
+               // Create and allocate CVODES memory for backward run\r
+               System.out.printf("Create and allocate CVODES memory for backward run\n");\r
+               \r
+               // Call CVodeCreateB to specify the solution method for the backward problem.\r
+               flag = CVodeCreateB(cvode_mem, CV_BDF, CV_NEWTON, indexB);\r
+               if (flag != CV_SUCCESS) {\r
+                       return;\r
+               }\r
+               \r
+               // Call CVodeInitB to allocate internal memory and initialize the\r
+               // backward problem.\r
+               flag = CVodeInitB(cvode_mem, indexB[0], fB, TB1, yB);\r
        }\r
        \r
        /*\r
@@ -553,7 +622,7 @@ public abstract class CVODES_ASA extends CVODES {
            \r
            ca_mem.ca_IMmalloc = CVAhermiteMalloc_select;\r
            //ca_mem.ca_IMfree   = CVAhermiteFree;\r
-           //ca_mem.ca_IMget    = CVAhermiteGetY;\r
+           ca_mem.ca_IMget    = CVAhermiteGetY_select;\r
            ca_mem.ca_IMstore  = CVAhermiteStorePnt_select;\r
 \r
            break;\r
@@ -562,7 +631,7 @@ public abstract class CVODES_ASA extends CVODES {
          \r
            ca_mem.ca_IMmalloc = CVApolynomialMalloc_select;\r
            //ca_mem.ca_IMfree   = CVApolynomialFree;\r
-           //ca_mem.ca_IMget    = CVApolynomialGetY;\r
+           ca_mem.ca_IMget    = CVApolynomialGetY_select;\r
            ca_mem.ca_IMstore  = CVApolynomialStorePnt_select;\r
 \r
            break;\r
@@ -1250,7 +1319,6 @@ public abstract class CVODES_ASA extends CVODES {
                  /* Allocate space for ckdata */\r
                  ck_mem = null;\r
                  ck_mem = new CkpntMemRec();\r
-                 if (ck_mem == null) return(null);\r
 \r
                  /* Set cv_next to NULL */\r
                  ck_mem.ck_next = null;\r
@@ -1432,6 +1500,572 @@ public abstract class CVODES_ASA extends CVODES {
 \r
                  return(ck_mem);\r
                }\r
+               \r
+               /* \r
+                * CVodeGetQuad\r
+                *\r
+                * This routine extracts quadrature solution into yQout at the\r
+                * time which CVode returned the solution.\r
+                * This is just a wrapper that calls CVodeGetQuadDky with k=0.\r
+                */\r
+                \r
+               private int CVodeGetQuad(CVodeMemRec cv_mem, double tret[], NVector yQout)\r
+               {\r
+                 int flag;\r
+\r
+                 if (cv_mem == null) {\r
+                   cvProcessError(null, CV_MEM_NULL, "CVODES", "CVodeGetQuad", MSGCV_NO_MEM);\r
+                   return(CV_MEM_NULL);\r
+                 } \r
+\r
+                 tret[0] = cv_mem.cv_tretlast;\r
+                 \r
+                 flag = CVodeGetQuadDky(cv_mem,cv_mem.cv_tretlast,0,yQout);\r
+\r
+                 return(flag);\r
+               }\r
+               \r
+               /*\r
+                * CVodeGetQuadDky\r
+                *\r
+                * CVodeQuadDky computes the kth derivative of the yQ function at\r
+                * time t, where tn-hu <= t <= tn, tn denotes the current         \r
+                * internal time reached, and hu is the last internal step size   \r
+                * successfully used by the solver. The user may request \r
+                * k=0, 1, ..., qu, where qu is the current order. \r
+                * The derivative vector is returned in dky. This vector \r
+                * must be allocated by the caller. It is only legal to call this         \r
+                * function after a successful return from CVode with quadrature\r
+                * computation enabled.\r
+                */\r
+\r
+               private int CVodeGetQuadDky(CVodeMemRec cv_mem, double t, int k, NVector dkyQ)\r
+               { \r
+                 double s, c, r;\r
+                 double tfuzz, tp, tn1;\r
+                 int i, j;\r
+                 \r
+                 /* Check all inputs for legality */\r
+                 \r
+                 if (cv_mem == null) {\r
+                   cvProcessError(null, CV_MEM_NULL, "CVODES", "CVodeGetQuadDky", MSGCV_NO_MEM);\r
+                   return(CV_MEM_NULL);\r
+                 } \r
+\r
+                 if(cv_mem.cv_quadr != true) {\r
+                   cvProcessError(cv_mem, CV_NO_QUAD, "CVODES", "CVodeGetQuadDky", MSGCV_NO_QUAD);\r
+                   return(CV_NO_QUAD);\r
+                 }\r
+\r
+                 if (dkyQ == null) {\r
+                   cvProcessError(cv_mem, CV_BAD_DKY, "CVODES", "CVodeGetQuadDky", MSGCV_NULL_DKY);\r
+                   return(CV_BAD_DKY);\r
+                 }\r
+                 \r
+                 if ((k < 0) || (k > cv_mem.cv_q)) {\r
+                   cvProcessError(cv_mem, CV_BAD_K, "CVODES", "CVodeGetQuadDky", MSGCV_BAD_K);\r
+                   return(CV_BAD_K);\r
+                 }\r
+                 \r
+                 /* Allow for some slack */\r
+                 tfuzz = FUZZ_FACTOR * cv_mem.cv_uround *\r
+                   (Math.abs(cv_mem.cv_tn) + Math.abs(cv_mem.cv_hu));\r
+                 if (cv_mem.cv_hu < ZERO) tfuzz = -tfuzz;\r
+                 tp = cv_mem.cv_tn - cv_mem.cv_hu - tfuzz;\r
+                 tn1 = cv_mem.cv_tn + tfuzz;\r
+                 if ((t-tp)*(t-tn1) > ZERO) {\r
+                   cvProcessError(cv_mem, CV_BAD_T, "CVODES", "CVodeGetQuadDky", MSGCV_BAD_T);\r
+                   return(CV_BAD_T);\r
+                 }\r
+                 \r
+                 /* Sum the differentiated interpolating polynomial */\r
+                 \r
+                 s = (t - cv_mem.cv_tn) / cv_mem.cv_h;\r
+                 for (j=cv_mem.cv_q; j >= k; j--) {\r
+                   c = ONE;\r
+                   for (i=j; i >= j-k+1; i--) c *= i;\r
+                   if (j == cv_mem.cv_q) {\r
+                     N_VScale_Serial(c, cv_mem.cv_znQ[cv_mem.cv_q], dkyQ);\r
+                   } else {\r
+                     N_VLinearSum_Serial(c, cv_mem.cv_znQ[j], s, dkyQ, dkyQ);\r
+                   }\r
+                 }\r
+                 if (k == 0) return(CV_SUCCESS);\r
+                 r = Math.pow(cv_mem.cv_h, -k);\r
+                 N_VScale_Serial(r, dkyQ, dkyQ);\r
+                 return(CV_SUCCESS);\r
+                 \r
+               }\r
+\r
+               class CVadjCheckPointRec {\r
+                         CkpntMemRec my_addr;\r
+                         CkpntMemRec next_addr;\r
+                         double t0;\r
+                         double t1;\r
+                         long nstep;\r
+                         int order;\r
+                         double step;\r
+                       };\r
+\r
+                       /*\r
+                        * CVodeGetAdjCheckPointsInfo\r
+                        *\r
+                        * This routine loads an array of nckpnts structures of type CVadjCheckPointRec.\r
+                        * The user must allocate space for ckpnt.\r
+                        */\r
+\r
+                       private int CVodeGetAdjCheckPointsInfo(CVodeMemRec cv_mem, CVadjCheckPointRec ckpnt[])\r
+                       {\r
+                         CVadjMemRec ca_mem;\r
+                         CkpntMemRec ck_mem;\r
+                         int i;\r
+\r
+                         /* Check if cvode_mem exists */\r
+                         if (cv_mem == null) {\r
+                           cvProcessError(null, CV_MEM_NULL, "CVODEA", "CVodeGetAdjCheckPointsInfo", MSGCV_NO_MEM);\r
+                           return(CV_MEM_NULL);\r
+                         }\r
+\r
+                         /* Was ASA initialized? */\r
+                         if (cv_mem.cv_adjMallocDone == false) {\r
+                           cvProcessError(cv_mem, CV_NO_ADJ, "CVODEA", "CVodeGetAdjCheckPointsInfo", MSGCV_NO_ADJ);\r
+                           return(CV_NO_ADJ);\r
+                         } \r
+                         ca_mem = cv_mem.cv_adj_mem;\r
+\r
+                         ck_mem = ca_mem.ck_mem;\r
+\r
+                         i = 0;\r
+\r
+                         while (ck_mem != null) {\r
+\r
+                           ckpnt[i].my_addr = ck_mem;\r
+                           ckpnt[i].next_addr = ck_mem.ck_next;\r
+                           ckpnt[i].t0 = ck_mem.ck_t0;\r
+                           ckpnt[i].t1 = ck_mem.ck_t1;\r
+                           ckpnt[i].nstep = ck_mem.ck_nst;\r
+                           ckpnt[i].order = ck_mem.ck_q;\r
+                           ckpnt[i].step = ck_mem.ck_h;\r
+\r
+                           ck_mem = ck_mem.ck_next;\r
+                           i++;\r
+\r
+                         }\r
+\r
+                         return(CV_SUCCESS);\r
+\r
+                       }\r
+                       \r
+                       private int CVodeCreateB(CVodeMemRec cv_mem, int lmmB, int iterB, int which[])\r
+                       {\r
+                         CVadjMemRec ca_mem;\r
+                         CVodeBMemRec new_cvB_mem;\r
+                         CVodeMemRec cvodeB_mem;\r
+\r
+                         /* Check if cvode_mem exists */\r
+                         if (cv_mem == null) {\r
+                           cvProcessError(null, CV_MEM_NULL, "CVODEA", "CVodeCreateB", MSGCV_NO_MEM);\r
+                           return(CV_MEM_NULL);\r
+                         }\r
+\r
+                         /* Was ASA initialized? */\r
+                         if (cv_mem.cv_adjMallocDone == false) {\r
+                           cvProcessError(cv_mem, CV_NO_ADJ, "CVODEA", "CVodeCreateB", MSGCV_NO_ADJ);\r
+                           return(CV_NO_ADJ);\r
+                         }\r
+                         ca_mem = cv_mem.cv_adj_mem;\r
+\r
+                         /* Allocate space for new CVodeBMem object */\r
+\r
+                         new_cvB_mem = null;\r
+                         new_cvB_mem = new CVodeBMemRec();\r
+\r
+                         /* Create and set a new CVODES object for the backward problem */\r
+\r
+                         cvodeB_mem = CVodeCreate(lmmB, iterB);\r
+                         if (cvodeB_mem == null) {\r
+                           cvProcessError(cv_mem, CV_MEM_FAIL, "CVODEA", "CVodeCreateB", MSGCV_MEM_FAIL);\r
+                           return(CV_MEM_FAIL);\r
+                         }\r
+\r
+                         UserData userData = new UserData();\r
+                         userData.memRec = cv_mem;\r
+                         cvodeB_mem.cv_user_data = userData;\r
+\r
+                         cvodeB_mem.cv_mxhnil = -1;\r
+\r
+                         //CVodeSetErrHandlerFn(cvodeB_mem, cv_mem->cv_ehfun, cv_mem->cv_eh_data);\r
+                         //cvodeB_mem.cv_ehfun = cv_mem.cv_ehfun;\r
+                         cvodeB_mem.cv_eh_data = cv_mem.cv_eh_data;\r
+\r
+                         //CVodeSetErrFile(cvodeB_mem, cv_mem->cv_errfp);\r
+                         cvodeB_mem.cv_errfp = cv_mem.cv_errfp;\r
+\r
+                         /* Set/initialize fields in the new CVodeBMem object, new_cvB_mem */\r
+\r
+                         new_cvB_mem.cv_index   = ca_mem.ca_nbckpbs;\r
+\r
+                         new_cvB_mem.cv_mem     = cvodeB_mem;\r
+\r
+                         //new_cvB_mem.cv_f       = null;\r
+                         //new_cvB_mem.cv_fs      = null;\r
+\r
+                         //new_cvB_mem.cv_fQ      = null;\r
+                         //new_cvB_mem.cv_fQs     = null;\r
+\r
+                         new_cvB_mem.cv_user_data  = null;\r
+\r
+                         //new_cvB_mem.cv_lmem    = null;\r
+                         //new_cvB_mem.cv_lfree   = null;\r
+                         //new_cvB_mem.cv_pmem    = null;\r
+                         //new_cvB_mem.cv_pfree   = null;\r
+\r
+                         new_cvB_mem.cv_y       = null;\r
+\r
+                         new_cvB_mem.cv_f_withSensi = false;\r
+                         new_cvB_mem.cv_fQ_withSensi = false;\r
+\r
+                         /* Attach the new object to the linked list cvB_mem */\r
+\r
+                         new_cvB_mem.cv_next = ca_mem.cvB_mem;\r
+                         ca_mem.cvB_mem = new_cvB_mem;\r
+                         \r
+                         /* Return the index of the newly created CVodeBMem object.\r
+                          * This must be passed to CVodeInitB and to other ***B \r
+                          * functions to set optional inputs for this backward problem */\r
+\r
+                         which[0] = ca_mem.ca_nbckpbs;\r
+\r
+                         ca_mem.ca_nbckpbs++;\r
+\r
+                         return(CV_SUCCESS);\r
+                       }\r
+\r
+                       private int CVodeInitB(CVodeMemRec cv_mem, int which, \r
+                              int fB,\r
+                              double tB0, NVector yB0)\r
+               {\r
+                 CVadjMemRec ca_mem;\r
+                 CVodeBMemRec cvB_mem;\r
+                 CVodeMemRec cvodeB_mem;\r
+                 int flag;\r
+\r
+                 /* Check if cvode_mem exists */\r
+\r
+                 if (cv_mem == null) {\r
+                   cvProcessError(null, CV_MEM_NULL, "CVODEA", "CVodeInitB", MSGCV_NO_MEM);\r
+                   return(CV_MEM_NULL);\r
+                 }\r
+\r
+                 /* Was ASA initialized? */\r
+\r
+                 if (cv_mem.cv_adjMallocDone == false) {\r
+                   cvProcessError(cv_mem, CV_NO_ADJ, "CVODEA", "CVodeInitB", MSGCV_NO_ADJ);\r
+                   return(CV_NO_ADJ);\r
+                 } \r
+                 ca_mem = cv_mem.cv_adj_mem;\r
+\r
+                 /* Check the value of which */\r
+\r
+                 if ( which >= ca_mem.ca_nbckpbs ) {\r
+                   cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeInitB", MSGCV_BAD_WHICH);\r
+                   return(CV_ILL_INPUT);\r
+                 }\r
+\r
+                 /* Find the CVodeBMem entry in the linked list corresponding to which */\r
+\r
+                 cvB_mem = ca_mem.cvB_mem;\r
+                 while (cvB_mem != null) {\r
+                   if ( which == cvB_mem.cv_index ) break;\r
+                   cvB_mem = cvB_mem.cv_next;\r
+                 }\r
+\r
+                 cvodeB_mem = cvB_mem.cv_mem;\r
+                 \r
+                 /* Allocate and set the CVODES object */\r
+\r
+                 flag = CVodeInit(cvodeB_mem, CVArhs_select, tB0, yB0);\r
+\r
+                 if (flag != CV_SUCCESS) return(flag);\r
+\r
+                 /* Copy fB function in cvB_mem */\r
+\r
+                 cvB_mem.cv_f_withSensi = false;\r
+                 //cvB_mem.cv_f = fB;\r
+\r
+                 /* Allocate space and initialize the y Nvector in cvB_mem */\r
+\r
+                 cvB_mem.cv_t0 = tB0;\r
+                 cvB_mem.cv_y = N_VClone(yB0);\r
+                 N_VScale_Serial(ONE, yB0, cvB_mem.cv_y);\r
+\r
+                 return(CV_SUCCESS);\r
+               }\r
+\r
+                       /*\r
+                        * CVArhs\r
+                        *\r
+                        * This routine interfaces to the CVRhsFnB (or CVRhsFnBS) routine \r
+                        * provided by the user.\r
+                        */\r
+\r
+\r
+                       private int CVArhs(double t, NVector yB, \r
+                                         NVector yBdot, CVodeMemRec cv_mem)\r
+                       {\r
+                         CVadjMemRec ca_mem;\r
+                         CVodeBMemRec cvB_mem;\r
+                         int flag = 0;\r
+                         int retval = 0;\r
+\r
+                         ca_mem = cv_mem.cv_adj_mem;\r
+\r
+                         cvB_mem = ca_mem.ca_bckpbCrt;\r
+\r
+                         /* Get forward solution from interpolation */\r
+\r
+                         if (ca_mem.ca_IMinterpSensi)\r
+                               if (ca_mem.ca_IMget == CVAhermiteGetY_select) {\r
+                               flag = CVAhermiteGetY(cv_mem, t, ca_mem.ca_ytmp, ca_mem.ca_yStmp);\r
+                               }\r
+                         else \r
+                                 if (ca_mem.ca_IMget == CVAhermiteGetY_select) {\r
+                                 flag = CVAhermiteGetY(cv_mem, t, ca_mem.ca_ytmp, null);\r
+                                 }\r
+\r
+                         if (flag != CV_SUCCESS) {\r
+                           cvProcessError(cv_mem, -1, "CVODEA", "CVArhs", MSGCV_BAD_TINTERP, t);\r
+                           return(-1);\r
+                         }\r
+\r
+                         /* Call the user's RHS function */\r
+\r
+                         if (cvB_mem.cv_f_withSensi)\r
+                               if (testMode) {\r
+                                       //retval = fBS1TestMode(t, ca_mem.ca_ytmp, ca_mem.ca_yStmp, yB, yBdot, cvB_mem.cv_user_data);   \r
+                               }\r
+                               else {\r
+                               //retval = fBS1(t, ca_mem.ca_ytmp, ca_mem.ca_yStmp, yB, yBdot, cvB_mem.cv_user_data);\r
+                               }\r
+                         else\r
+                                 if (testMode) {\r
+                                         retval = fBTestMode(t, ca_mem.ca_ytmp, yB, yBdot, cvB_mem.cv_user_data);\r
+                                 }\r
+                                 else {\r
+                                 retval = fB(t, ca_mem.ca_ytmp, yB, yBdot, cvB_mem.cv_user_data);\r
+                                 }\r
+\r
+                         return(retval);\r
+                       }\r
+                       \r
+                       /*\r
+                        * CVAhermiteGetY ( -> IMget )\r
+                        *\r
+                        * This routine uses cubic piece-wise Hermite interpolation for \r
+                        * the forward solution vector. \r
+                        * It is typically called by the wrapper routines before calling\r
+                        * user provided routines (fB, djacB, bjacB, jtimesB, psolB) but\r
+                        * can be directly called by the user through CVodeGetAdjY\r
+                        */\r
+\r
+                       private int CVAhermiteGetY(CVodeMemRec cv_mem, double t,\r
+                                                 NVector y, NVector yS[])\r
+                       {\r
+                         CVadjMemRec ca_mem;\r
+                         DtpntMemRec dt_mem[];\r
+                         HermiteDataMemRec content0, content1;\r
+\r
+                         double t0, t1, delta;\r
+                         double factor1, factor2, factor3;\r
+\r
+                         NVector y0, yd0, y1, yd1;\r
+                         NVector yS0[]=null;\r
+                         NVector ySd0[]= null;\r
+                         NVector yS1[], ySd1[];\r
+\r
+                         int flag, is, NS;\r
+                         long indx[] = new long[1];\r
+                         boolean newpoint[] = new boolean[1];\r
+\r
+                        \r
+                         ca_mem = cv_mem.cv_adj_mem;\r
+                         dt_mem = ca_mem.dt_mem;\r
+                        \r
+                         /* Local value of Ns */\r
+                        \r
+                         NS = (ca_mem.ca_IMinterpSensi && (yS != null)) ? cv_mem.cv_Ns : 0;\r
+\r
+                         /* Get the index in dt_mem */\r
+\r
+                         flag = CVAfindIndex(cv_mem, t, indx, newpoint);\r
+                         if (flag != CV_SUCCESS) return(flag);\r
+\r
+                         /* If we are beyond the left limit but close enough,\r
+                            then return y at the left limit. */\r
+\r
+                         if (indx[0] == 0) {\r
+                           content0 = dt_mem[0].hermiteContent;\r
+                           N_VScale_Serial(ONE, content0.y, y);\r
+                           for (is=0; is<NS; is++) N_VScale_Serial(ONE, content0.yS[is], yS[is]);\r
+                           return(CV_SUCCESS);\r
+                         }\r
+\r
+                         /* Extract stuff from the appropriate data points */\r
+\r
+                         t0 = dt_mem[(int)(indx[0]-1)].t;\r
+                         t1 = dt_mem[(int)indx[0]].t;\r
+                         delta = t1 - t0;\r
+\r
+                         content0 = dt_mem[(int)(indx[0]-1)].hermiteContent;\r
+                         y0  = content0.y;\r
+                         yd0 = content0.yd;\r
+                         if (ca_mem.ca_IMinterpSensi) {\r
+                           yS0  = content0.yS;\r
+                           ySd0 = content0.ySd;\r
+                         }\r
+\r
+                         if (newpoint[0]) {\r
+                           \r
+                           /* Recompute Y0 and Y1 */\r
+\r
+                           content1 = dt_mem[(int)indx[0]].hermiteContent;\r
+\r
+                           y1  = content1.y;\r
+                           yd1 = content1.yd;\r
+\r
+                           N_VLinearSum_Serial(ONE, y1, -ONE, y0, ca_mem.ca_Y[0]);\r
+                           N_VLinearSum_Serial(ONE, yd1,  ONE, yd0, ca_mem.ca_Y[1]);\r
+                           N_VLinearSum_Serial(delta, ca_mem.ca_Y[1], -TWO, ca_mem.ca_Y[0], ca_mem.ca_Y[1]);\r
+                           N_VLinearSum_Serial(ONE, ca_mem.ca_Y[0], -delta, yd0, ca_mem.ca_Y[0]);\r
+\r
+\r
+                           yS1  = content1.yS;\r
+                           ySd1 = content1.ySd;\r
+                             \r
+                           for (is=0; is<NS; is++) {\r
+                             N_VLinearSum_Serial(ONE, yS1[is], -ONE, yS0[is], ca_mem.ca_YS[0][is]);\r
+                             N_VLinearSum_Serial(ONE, ySd1[is],  ONE, ySd0[is], ca_mem.ca_YS[1][is]);\r
+                             N_VLinearSum_Serial(delta, ca_mem.ca_YS[1][is], -TWO, ca_mem.ca_YS[0][is], ca_mem.ca_YS[1][is]);\r
+                             N_VLinearSum_Serial(ONE, ca_mem.ca_YS[0][is], -delta, ySd0[is], ca_mem.ca_YS[0][is]);\r
+                           }\r
+\r
+                         }\r
+\r
+                         /* Perform the actual interpolation. */\r
+\r
+                         factor1 = t - t0;\r
+\r
+                         factor2 = factor1/delta;\r
+                         factor2 = factor2*factor2;\r
+\r
+                         factor3 = factor2*(t-t1)/delta;\r
+\r
+                         N_VLinearSum_Serial(ONE, y0, factor1, yd0, y);\r
+                         N_VLinearSum_Serial(ONE, y, factor2, ca_mem.ca_Y[0], y);\r
+                         N_VLinearSum_Serial(ONE, y, factor3, ca_mem.ca_Y[1], y);\r
+\r
+                         for (is=0; is<NS; is++) {\r
+                           N_VLinearSum_Serial(ONE, yS0[is], factor1, ySd0[is], yS[is]);\r
+                           N_VLinearSum_Serial(ONE, yS[is], factor2, ca_mem.ca_YS[0][is], yS[is]);\r
+                           N_VLinearSum_Serial(ONE, yS[is], factor3, ca_mem.ca_YS[1][is], yS[is]);\r
+                         }\r
+\r
+\r
+                         return(CV_SUCCESS);\r
+                       }\r
+                       \r
+                       /*\r
+                        * CVAfindIndex\r
+                        *\r
+                        * Finds the index in the array of data point strctures such that\r
+                        *     dt_mem[indx-1].t <= t < dt_mem[indx].t\r
+                        * If indx is changed from the previous invocation, then newpoint = SUNTRUE\r
+                        *\r
+                        * If t is beyond the leftmost limit, but close enough, indx=0.\r
+                        *\r
+                        * Returns CV_SUCCESS if successful and CV_GETY_BADT if unable to\r
+                        * find indx (t is too far beyond limits).\r
+                        */\r
+\r
+                       static int CVAfindIndex(CVodeMemRec cv_mem, double t, \r
+                                               long indx[], boolean newpoint[])\r
+                       {\r
+                         CVadjMemRec ca_mem;\r
+                         long ilast = 0;\r
+                         DtpntMemRec dt_mem[];\r
+                         int sign;\r
+                         boolean to_left, to_right;\r
+\r
+                         ca_mem = cv_mem.cv_adj_mem;\r
+                         dt_mem = ca_mem.dt_mem;\r
+\r
+                         newpoint[0] = false;\r
+\r
+                         /* Find the direction of integration */\r
+                         sign = (ca_mem.ca_tfinal - ca_mem.ca_tinitial > ZERO) ? 1 : -1;\r
+\r
+                         /* If this is the first time we use new data */\r
+                         if (ca_mem.ca_IMnewData) {\r
+                           ilast     = ca_mem.ca_np-1;\r
+                           newpoint[0] = true;\r
+                           ca_mem.ca_IMnewData   = false;\r
+                         }\r
+\r
+                         /* Search for indx starting from ilast */\r
+                         to_left  = ( sign*(t - dt_mem[(int)(ilast-1)].t) < ZERO);\r
+                         to_right = ( sign*(t - dt_mem[(int)(ilast)].t)   > ZERO);\r
+\r
+                         if ( to_left ) {\r
+                           /* look for a new indx to the left */\r
+\r
+                           newpoint[0] = true;\r
+                           \r
+                           indx[0] = ilast;\r
+                           for(;;) {\r
+                             if ( indx[0] == 0 ) break;\r
+                             if ( sign*(t - dt_mem[(int)(indx[0]-1)].t) <= ZERO ) (indx[0])--;\r
+                             else                                         break;\r
+                           }\r
+\r
+                           if ( indx[0] == 0 )\r
+                             ilast = 1;\r
+                           else\r
+                             ilast = indx[0];\r
+\r
+                           if ( indx[0] == 0 ) {\r
+                             /* t is beyond leftmost limit. Is it too far? */  \r
+                             if ( Math.abs(t - dt_mem[0].t) > FUZZ_FACTOR * cv_mem.cv_uround ) {\r
+                               return(CV_GETY_BADT);\r
+                             }\r
+                           }\r
+\r
+                         } else if ( to_right ) {\r
+                           /* look for a new indx to the right */\r
+\r
+                           newpoint[0] = true;\r
+\r
+                           indx [0]= ilast;\r
+                           for(;;) {\r
+                             if ( sign*(t - dt_mem[(int)indx[0]].t) > ZERO) (indx[0])++;\r
+                             else                                     break;\r
+                           }\r
+\r
+                           ilast = indx[0];\r
+\r
+\r
+                         } else {\r
+                           /* ilast is still OK */\r
+\r
+                           indx[0] = ilast;\r
+\r
+                         }\r
+\r
+                         return(CV_SUCCESS);\r
+\r
+\r
+                       }\r
+\r
+\r
 \r
 \r
 }
\ No newline at end of file