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

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

index a2f8050..0aba980 100644 (file)
@@ -160,6 +160,7 @@ public abstract class CVODES_ASA extends CVODES {
                int fB = cvsRoberts_ASAi_dns;\r
                double AB[][];\r
                SUNLinearSolver LSB;\r
+               long nstB[] = new long[1];\r
                \r
                // Print problem description \r
                System.out.printf("\nAdjoint Sensitivity Example for Chemical Kinetics\n");\r
@@ -429,6 +430,81 @@ public abstract class CVODES_ASA extends CVODES {
                \r
                // Call CVodeBto integrate the backward ODE problem.\r
                flag = CVodeB(cvode_mem, TBout1, CV_NORMAL);\r
+               if (flag < 0) {\r
+                       return;\r
+               }\r
+               \r
+               // Call CVodeGetB to get yB of the backward ODE problem\r
+               flag = CVodeGetB(cvode_mem, indexB[0], time, yB);\r
+               if (flag != CV_SUCCESS) {\r
+                       return;\r
+               }\r
+               \r
+               // Call CvodeGetAdjY to get the interpolated value of the forward solution\r
+               // y during a backward integration.\r
+               flag = CVodeGetAdjY(cvode_mem, TBout1, y);\r
+               if (flag < 0) {\r
+                       return;\r
+               }\r
+               \r
+               System.out.printf("returned t: " + time[0] + "\n");\r
+               System.out.printf("tout: " + TBout1 + "\n");\r
+               System.out.printf("lambda(t): " + yB.data[0] + "  " + yB.data[1] + "  " + yB.data[2] + "\n");\r
+               System.out.printf("y(t): " + y.data[0] + "  " + y.data[1] + "  " + y.data[2] + "\n");\r
+               \r
+               // Then at t = T0\r
+               \r
+               flag = CVodeB(cvode_mem, T0, CV_NORMAL);\r
+               if (flag < 0) {\r
+                       return;\r
+               }\r
+               \r
+               CVodeGetNumSteps(CVodeGetAdjCVodeBmem(cvode_mem, indexB[0]), nstB);\r
+               System.out.printf("Done (nst = " + nstB[0] + ")\n");\r
+               \r
+               flag = CVodeGetB(cvode_mem, indexB[0], time, yB);\r
+               if (flag != CV_SUCCESS) {\r
+                       return;\r
+               }\r
+               \r
+               // Call CVodeGetQuadB to get the quadrature solution vector after a\r
+               // successful return from CVodeB.\r
+               flag = CVodeGetQuadB(cvode_mem, indexB[0], time, qB);\r
+               if (flag < 0) {\r
+                       return;\r
+               }\r
+               \r
+               flag = CVodeGetAdjY(cvode_mem, T0, y);\r
+               if (flag < 0) {\r
+                       return;\r
+               }\r
+               \r
+               System.out.printf("returned t: " + time[0] + "\n");\r
+               System.out.printf("lambda(t0): " + yB.data[0] + "  " + yB.data[1] + "  " + yB.data[2] + "\n");\r
+               System.out.printf("y(t0): " + y.data[0] + "  " + y.data[1] + "  " + y.data[2] + "\n");\r
+               System.out.printf("dG/dp: " + (-qB.data[0]) + "  " + (-qB.data[1]) + "  " + (-qB.data[2]) + "\n");\r
+               \r
+               // Reinitialize backward phase (new tB0)\r
+               \r
+               yB.data[0] = ZERO;\r
+               yB.data[1] = ZERO;\r
+               yB.data[2] = ZERO;\r
+               \r
+               qB.data[0] = ZERO;\r
+               qB.data[1] = ZERO;\r
+               qB.data[2] = ZERO;\r
+               \r
+               System.out.printf("Re-initialize CVODES memory for backward run\n");\r
+               \r
+               flag = CVodeReInitB(cvode_mem, indexB[0], TB2, yB);\r
+               if (flag < 0) {\r
+                       return;\r
+               }\r
+               \r
+               flag = CVodeQuadReInitB(cvode_mem, indexB[0], qB);\r
+               if (flag != CV_SUCCESS) {\r
+                       return;\r
+               }\r
        } \r
        \r
        /*\r
@@ -2166,7 +2242,6 @@ public abstract class CVODES_ASA extends CVODES {
          CVadjMemRec ca_mem;\r
          CVodeBMemRec cvB_mem;\r
          CVodeMemRec cvodeB_mem;\r
-         int flag;\r
 \r
          /* Check if cvode_mem exists */\r
          if (cv_mem == null) {\r
@@ -2511,13 +2586,18 @@ public abstract class CVODES_ASA extends CVODES {
          dt_mem = ca_mem.dt_mem;\r
 \r
          /* Initialize cv_mem with data from ck_mem */\r
-         //flag = CVAckpntGet(cv_mem, ck_mem);\r
-         //if (flag != CV_SUCCESS)\r
-           //return(CV_REIFWD_FAIL);\r
+         flag = CVAckpntGet(cv_mem, ck_mem);\r
+         if (flag != CV_SUCCESS)\r
+           return(CV_REIFWD_FAIL);\r
 \r
          /* Set first structure in dt_mem[0] */\r
          dt_mem[0].t = ck_mem.ck_t0;\r
-         //ca_mem.ca_IMstore(cv_mem, dt_mem[0]);\r
+         if (ca_mem.ca_IMstore == CVAhermiteStorePnt_select) {\r
+               CVAhermiteStorePnt(cv_mem, dt_mem[0]);\r
+         }\r
+         else {\r
+                   CVApolynomialStorePnt(cv_mem, dt_mem[0]);\r
+         }\r
 \r
          /* Decide whether TSTOP must be activated */\r
          if (ca_mem.ca_tstopCVodeFcall) {\r
@@ -2535,7 +2615,12 @@ public abstract class CVODES_ASA extends CVODES {
            if (flag < 0) return(CV_FWD_FAIL);\r
 \r
            dt_mem[(int)i].t = t[0];\r
-           //ca_mem.ca_IMstore(cv_mem, dt_mem[(int)i]);\r
+           if (ca_mem.ca_IMstore == CVAhermiteStorePnt_select) {\r
+               CVAhermiteStorePnt(cv_mem, dt_mem[(int)i]);\r
+         }\r
+         else {\r
+                   CVApolynomialStorePnt(cv_mem, dt_mem[(int)i]);\r
+         }\r
            i++;\r
 \r
          } while ( sign*(ck_mem.ck_t1 - t[0]) > ZERO );\r
@@ -2581,6 +2666,678 @@ public abstract class CVODES_ASA extends CVODES {
          return(CV_SUCCESS);\r
        }\r
 \r
+       /*\r
+        * CVAckpntGet\r
+        *\r
+        * This routine prepares CVODES for a hot restart from\r
+        * the check point ck_mem\r
+        */\r
+\r
+       private int CVAckpntGet(CVodeMemRec cv_mem, CkpntMemRec ck_mem) \r
+       {\r
+         int flag, j, is, qmax;\r
+\r
+         if (ck_mem.ck_next == null) {\r
+\r
+           /* In this case, we just call the reinitialization routine,\r
+            * but make sure we use the same initial stepsize as on \r
+            * the first run. */\r
+\r
+           // Specify the initial step size\r
+                 cv_mem.cv_hin = cv_mem.cv_h0u;\r
+\r
+           flag = CVodeReInit(cv_mem, ck_mem.ck_t0, ck_mem.ck_zn[0]);\r
+           if (flag != CV_SUCCESS) return(flag);\r
+\r
+           if (ck_mem.ck_quadr) {\r
+             flag = CVodeQuadReInit(cv_mem, ck_mem.ck_znQ[0]);\r
+             if (flag != CV_SUCCESS) return(flag);\r
+           }\r
+\r
+           if (ck_mem.ck_sensi) {\r
+             flag = CVodeSensReInit(cv_mem, cv_mem.cv_ism, ck_mem.ck_znS[0]);\r
+             if (flag != CV_SUCCESS) return(flag);\r
+           }\r
+\r
+           if (ck_mem.ck_quadr_sensi) {\r
+             flag = CVodeQuadSensReInit(cv_mem, ck_mem.ck_znQS[0]);\r
+             if (flag != CV_SUCCESS) return(flag);\r
+           }\r
+\r
+         } else {\r
+           \r
+           qmax = cv_mem.cv_qmax;\r
+\r
+           /* Copy parameters from check point data structure */\r
+\r
+           cv_mem.cv_nst       = ck_mem.ck_nst;\r
+           cv_mem.cv_tretlast  = ck_mem.ck_tretlast;\r
+           cv_mem.cv_q         = ck_mem.ck_q;\r
+           cv_mem.cv_qprime    = ck_mem.ck_qprime;\r
+           cv_mem.cv_qwait     = ck_mem.ck_qwait;\r
+           cv_mem.cv_L         = ck_mem.ck_L;\r
+           cv_mem.cv_gammap    = ck_mem.ck_gammap;\r
+           cv_mem.cv_h         = ck_mem.ck_h;\r
+           cv_mem.cv_hprime    = ck_mem.ck_hprime;\r
+           cv_mem.cv_hscale    = ck_mem.ck_hscale;\r
+           cv_mem.cv_eta       = ck_mem.ck_eta;\r
+           cv_mem.cv_etamax    = ck_mem.ck_etamax;\r
+           cv_mem.cv_tn        = ck_mem.ck_t0;\r
+           cv_mem.cv_saved_tq5 = ck_mem.ck_saved_tq5;\r
+           \r
+           /* Copy the arrays from check point data structure */\r
+\r
+           for (j=0; j<=cv_mem.cv_q; j++) N_VScale_Serial(ONE, ck_mem.ck_zn[j], cv_mem.cv_zn[j]);\r
+           if ( cv_mem.cv_q < qmax ) N_VScale_Serial(ONE, ck_mem.ck_zn[qmax], cv_mem.cv_zn[qmax]);\r
+\r
+           if (ck_mem.ck_quadr) {\r
+             for (j=0; j<=cv_mem.cv_q; j++) N_VScale_Serial(ONE, ck_mem.ck_znQ[j], cv_mem.cv_znQ[j]);\r
+             if ( cv_mem.cv_q < qmax ) N_VScale_Serial(ONE, ck_mem.ck_znQ[qmax], cv_mem.cv_znQ[qmax]);\r
+           }\r
+\r
+           if (ck_mem.ck_sensi) {\r
+             for (is=0; is<cv_mem.cv_Ns; is++) {\r
+               for (j=0; j<=cv_mem.cv_q; j++) N_VScale_Serial(ONE, ck_mem.ck_znS[j][is], cv_mem.cv_znS[j][is]);\r
+               if ( cv_mem.cv_q < qmax ) N_VScale_Serial(ONE, ck_mem.ck_znS[qmax][is], cv_mem.cv_znS[qmax][is]);\r
+             }\r
+           }\r
+\r
+           if (ck_mem.ck_quadr_sensi) {\r
+             for (is=0; is<cv_mem.cv_Ns; is++) {\r
+               for (j=0; j<=cv_mem.cv_q; j++) N_VScale_Serial(ONE, ck_mem.ck_znQS[j][is], cv_mem.cv_znQS[j][is]);\r
+               if ( cv_mem.cv_q < qmax ) N_VScale_Serial(ONE, ck_mem.ck_znQS[qmax][is], cv_mem.cv_znQS[qmax][is]);\r
+             }\r
+           }\r
+\r
+           for (j=0; j<=L_MAX; j++)        cv_mem.cv_tau[j] = ck_mem.ck_tau[j];\r
+           for (j=0; j<=NUM_TESTS; j++)    cv_mem.cv_tq[j] = ck_mem.ck_tq[j];\r
+           for (j=0; j<=cv_mem.cv_q; j++) cv_mem.cv_l[j] = ck_mem.ck_l[j];\r
+           \r
+           /* Force a call to setup */\r
+\r
+           cv_mem.cv_forceSetup = true;\r
+\r
+         }\r
+\r
+         return(CV_SUCCESS);\r
+       }\r
+       \r
+       /*\r
+        * CVodeReInit\r
+        *\r
+        * CVodeReInit re-initializes CVODES's memory for a problem, assuming\r
+        * it has already been allocated in a prior CVodeInit call.\r
+        * All problem specification inputs are checked for errors.\r
+        * If any error occurs during initialization, it is reported to the\r
+        * file whose file pointer is errfp.\r
+        * The return value is CV_SUCCESS = 0 if no errors occurred, or\r
+        * a negative value otherwise.\r
+        */\r
+\r
+       private int CVodeReInit(CVodeMemRec cv_mem, double t0, NVector y0)\r
+       {\r
+         int i,k;\r
+        \r
+         /* Check cvode_mem */\r
+\r
+         if (cv_mem==null) {\r
+           cvProcessError(null, CV_MEM_NULL, "CVODES", "CVodeReInit",\r
+                          MSGCV_NO_MEM);\r
+           return(CV_MEM_NULL);\r
+         }\r
+\r
+         /* Check if cvode_mem was allocated */\r
+\r
+         if (cv_mem.cv_MallocDone == false) {\r
+           cvProcessError(cv_mem, CV_NO_MALLOC, "CVODES", "CVodeReInit",\r
+                          MSGCV_NO_MALLOC);\r
+           return(CV_NO_MALLOC);\r
+         }\r
+\r
+         /* Check for legal input parameters */\r
+\r
+         if (y0 == null) {\r
+           cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVodeReInit",\r
+                          MSGCV_NULL_Y0);\r
+           return(CV_ILL_INPUT);\r
+         }\r
+         \r
+         /* Copy the input parameters into CVODES state */\r
+\r
+         cv_mem.cv_tn = t0;\r
+\r
+         /* Set step parameters */\r
+\r
+         cv_mem.cv_q      = 1;\r
+         cv_mem.cv_L      = 2;\r
+         cv_mem.cv_qwait  = cv_mem.cv_L;\r
+         cv_mem.cv_etamax = ETAMX1;\r
 \r
+         cv_mem.cv_qu     = 0;\r
+         cv_mem.cv_hu     = ZERO;\r
+         cv_mem.cv_tolsf  = ONE;\r
+\r
+         /* Set forceSetup to SUNFALSE */\r
+\r
+         cv_mem.cv_forceSetup = false;\r
+\r
+         /* Initialize zn[0] in the history array */\r
+\r
+         N_VScale_Serial(ONE, y0, cv_mem.cv_zn[0]);\r
+        \r
+         /* Initialize all the counters */\r
+\r
+         cv_mem.cv_nst     = 0;\r
+         cv_mem.cv_nfe     = 0;\r
+         cv_mem.cv_ncfn[0]    = 0;\r
+         cv_mem.cv_netf[0]    = 0;\r
+         cv_mem.cv_nni     = 0;\r
+         cv_mem.cv_nsetups = 0;\r
+         cv_mem.cv_nhnil   = 0;\r
+         cv_mem.cv_nstlp   = 0;\r
+         cv_mem.cv_nscon   = 0;\r
+         cv_mem.cv_nge     = 0;\r
+\r
+         cv_mem.cv_irfnd   = 0;\r
+\r
+         /* Initialize other integrator optional outputs */\r
+\r
+         cv_mem.cv_h0u      = ZERO;\r
+         cv_mem.cv_next_h   = ZERO;\r
+         cv_mem.cv_next_q   = 0;\r
+\r
+         /* Initialize Stablilty Limit Detection data */\r
+\r
+         cv_mem.cv_nor = 0;\r
+         for (i = 1; i <= 5; i++)\r
+           for (k = 1; k <= 3; k++) \r
+             cv_mem.cv_ssdat[i-1][k-1] = ZERO;\r
+         \r
+         /* Problem has been successfully re-initialized */\r
+\r
+         return(CV_SUCCESS);\r
+       }\r
+\r
+       /*\r
+        * CVodeQuadReInit\r
+        *\r
+        * CVodeQuadReInit re-initializes CVODES's quadrature related memory \r
+        * for a problem, assuming it has already been allocated in prior \r
+        * calls to CVodeInit and CVodeQuadInit. \r
+        * All problem specification inputs are checked for errors.\r
+        * If any error occurs during initialization, it is reported to the\r
+        * file whose file pointer is errfp.\r
+        * The return value is CV_SUCCESS = 0 if no errors occurred, or\r
+        * a negative value otherwise.\r
+        */\r
+\r
+       private int CVodeQuadReInit(CVodeMemRec cv_mem, NVector yQ0)\r
+       {\r
+\r
+         /* Check cvode_mem */\r
+         if (cv_mem==null) {\r
+           cvProcessError(null, CV_MEM_NULL, "CVODES",\r
+                          "CVodeQuadReInit", MSGCV_NO_MEM);\r
+           return(CV_MEM_NULL);\r
+         }\r
+\r
+         /* Ckeck if quadrature was initialized? */\r
+         if (cv_mem.cv_QuadMallocDone == false) {\r
+           cvProcessError(cv_mem, CV_NO_QUAD, "CVODES",\r
+                          "CVodeQuadReInit", MSGCV_NO_QUAD);\r
+           return(CV_NO_QUAD);\r
+         }\r
+\r
+         /* Initialize znQ[0] in the history array */\r
+         N_VScale_Serial(ONE, yQ0, cv_mem.cv_znQ[0]);\r
+\r
+         /* Initialize counters */\r
+         cv_mem.cv_nfQe  = 0;\r
+         cv_mem.cv_netfQ[0] = 0;\r
+\r
+         /* Quadrature integration turned ON */\r
+         cv_mem.cv_quadr = true;\r
+\r
+         /* Quadrature re-initialization was successfull */\r
+         return(CV_SUCCESS);\r
+       }\r
+\r
+       /*\r
+        * CVodeSensReInit\r
+        *\r
+        * CVodeSensReInit re-initializes CVODES's sensitivity related memory \r
+        * for a problem, assuming it has already been allocated in prior \r
+        * calls to CVodeInit and CVodeSensInit/CVodeSensInit1. \r
+        * All problem specification inputs are checked for errors.\r
+        * The number of sensitivities Ns is assumed to be unchanged since\r
+        * the previous call to CVodeSensInit.\r
+        * If any error occurs during initialization, it is reported to the\r
+        * file whose file pointer is errfp.\r
+        * The return value is CV_SUCCESS = 0 if no errors occurred, or\r
+        * a negative value otherwise.\r
+        */ \r
+\r
+       private int CVodeSensReInit(CVodeMemRec cv_mem, int ism, NVector yS0[])\r
+       {\r
+         int is;  \r
+\r
+         /* Check cvode_mem */\r
+\r
+         if (cv_mem==null) {\r
+           cvProcessError(null, CV_MEM_NULL, "CVODES", "CVodeSensReInit",\r
+                          MSGCV_NO_MEM);\r
+           return(CV_MEM_NULL);\r
+         }\r
+\r
+         /* Was sensitivity initialized? */\r
+\r
+         if (cv_mem.cv_SensMallocDone == false) {\r
+           cvProcessError(cv_mem, CV_NO_SENS, "CVODES", "CVodeSensReInit",\r
+                          MSGCV_NO_SENSI);\r
+           return(CV_NO_SENS);\r
+         } \r
+\r
+         /* Check if ism is compatible */\r
+\r
+         if ((cv_mem.cv_ifS==CV_ALLSENS) && (ism==CV_STAGGERED1)) {\r
+           cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES",\r
+                          "CVodeSensReInit", MSGCV_BAD_ISM_IFS);\r
+           return(CV_ILL_INPUT);\r
+         }\r
+         \r
+         /* Check if ism is legal */\r
+\r
+         if ((ism!=CV_SIMULTANEOUS) && (ism!=CV_STAGGERED) && (ism!=CV_STAGGERED1)) {\r
+           cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES",\r
+                          "CVodeSensReInit", MSGCV_BAD_ISM);\r
+           return(CV_ILL_INPUT);\r
+         }\r
+         cv_mem.cv_ism = ism;\r
+\r
+         /* Check if yS0 is non-null */\r
+\r
+         if (yS0 == null) {\r
+           cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES",\r
+                          "CVodeSensReInit", MSGCV_NULL_YS0);\r
+           return(CV_ILL_INPUT);\r
+         }  \r
+\r
+         /* Allocate ncfS1, ncfnS1, and nniS1 if needed */\r
+\r
+         if ( (ism==CV_STAGGERED1) && (cv_mem.cv_stgr1alloc==false) ) {\r
+           cv_mem.cv_stgr1alloc = true;\r
+           cv_mem.cv_ncfS1 = null;\r
+           cv_mem.cv_ncfS1 = new int[cv_mem.cv_Ns][1];\r
+           cv_mem.cv_ncfnS1 = null;\r
+           cv_mem.cv_ncfnS1 = new long[cv_mem.cv_Ns][1];\r
+           cv_mem.cv_nniS1 = null;\r
+           cv_mem.cv_nniS1 = new long[cv_mem.cv_Ns];\r
+           if ( (cv_mem.cv_ncfS1==null) ||\r
+                (cv_mem.cv_ncfnS1==null) ||\r
+                (cv_mem.cv_nniS1==null) ) {\r
+             cvProcessError(cv_mem, CV_MEM_FAIL, "CVODES",\r
+                            "CVodeSensReInit", MSGCV_MEM_FAIL);\r
+             return(CV_MEM_FAIL);\r
+           }\r
+         }\r
+\r
+         /*---------------------------------------------- \r
+           All error checking is complete at this point \r
+           -----------------------------------------------*/\r
+\r
+         /* Initialize znS[0] in the history array */\r
+\r
+         for (is=0; is<cv_mem.cv_Ns; is++) \r
+           N_VScale_Serial(ONE, yS0[is], cv_mem.cv_znS[0][is]);\r
+\r
+         /* Initialize all sensitivity related counters */\r
+\r
+         cv_mem.cv_nfSe     = 0;\r
+         cv_mem.cv_nfeS     = 0;\r
+         cv_mem.cv_ncfnS[0]    = 0;\r
+         cv_mem.cv_netfS[0]    = 0;\r
+         cv_mem.cv_nniS     = 0;\r
+         cv_mem.cv_nsetupsS = 0;\r
+         if (ism==CV_STAGGERED1)\r
+           for (is=0; is<cv_mem.cv_Ns; is++) {\r
+             cv_mem.cv_ncfnS1[is][0] = 0;\r
+             cv_mem.cv_nniS1[is] = 0;\r
+           }\r
+\r
+         /* Problem has been successfully re-initialized */\r
+\r
+         cv_mem.cv_sensi = true;\r
+\r
+         return(CV_SUCCESS);\r
+       }\r
+       \r
+       /*\r
+        * CVodeQuadSensReInit\r
+        *\r
+        */\r
+\r
+       private int CVodeQuadSensReInit(CVodeMemRec cv_mem, NVector yQS0[])\r
+       {\r
+         int is;\r
+\r
+         /* Check cvode_mem */\r
+         if (cv_mem==null) {\r
+           cvProcessError(null, CV_MEM_NULL, "CVODES", "CVodeQuadSensReInit",\r
+                          MSGCV_NO_MEM);\r
+           return(CV_MEM_NULL);\r
+         }\r
+\r
+         /* Check if sensitivity analysis is active */\r
+         if (!cv_mem.cv_sensi) {\r
+           cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES",\r
+                          "CVodeQuadSensReInit", MSGCV_NO_SENSI);\r
+           return(CV_NO_SENS);\r
+         }\r
+\r
+         /* Was quadrature sensitivity initialized? */\r
+         if (cv_mem.cv_QuadSensMallocDone == false) {\r
+           cvProcessError(cv_mem, CV_NO_QUADSENS, "CVODES",\r
+                          "CVodeQuadSensReInit", MSGCV_NO_QUADSENSI);\r
+           return(CV_NO_QUADSENS);\r
+         } \r
+\r
+         /* Check if yQS0 is non-null */\r
+         if (yQS0 == null) {\r
+           cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES",\r
+                          "CVodeQuadSensReInit", MSGCV_NULL_YQS0);\r
+           return(CV_ILL_INPUT);\r
+         }\r
+\r
+         /*---------------------------------------------- \r
+           All error checking is complete at this point \r
+           -----------------------------------------------*/\r
+\r
+         /* Initialize znQS[0] in the history array */\r
+         for (is=0; is<cv_mem.cv_Ns; is++) \r
+           N_VScale_Serial(ONE, yQS0[is], cv_mem.cv_znQS[0][is]);\r
+\r
+         /* Initialize all sensitivity related counters */\r
+         cv_mem.cv_nfQSe  = 0;\r
+         cv_mem.cv_nfQeS  = 0;\r
+         cv_mem.cv_netfQS[0] = 0;\r
+\r
+         /* Quadrature sensitivities will be computed */\r
+         cv_mem.cv_quadr_sensi = true;\r
+\r
+         /* Problem has been successfully re-initialized */\r
+         return(CV_SUCCESS);\r
+       }\r
+\r
+       private int CVodeGetB(CVodeMemRec cv_mem, int which, double tret[], NVector yB)\r
+       {\r
+         CVadjMemRec ca_mem;\r
+         CVodeBMemRec cvB_mem;\r
+\r
+         /* Check if cvode_mem exists */\r
+         if (cv_mem == null) {\r
+           cvProcessError(null, CV_MEM_NULL, "CVODEA", "CVodeGetB", 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", "CVodeGetB", MSGCV_NO_ADJ);\r
+           return(CV_NO_ADJ);\r
+         } \r
+\r
+         ca_mem = cv_mem.cv_adj_mem;\r
+\r
+         /* Check the value of which */\r
+         if ( which >= ca_mem.ca_nbckpbs ) {\r
+           cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeGetB", MSGCV_BAD_WHICH);\r
+           return(CV_ILL_INPUT);\r
+         }\r
+\r
+         /* Find the CVodeBMem entry in the linked list corresponding to which */\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
+         N_VScale_Serial(ONE, cvB_mem.cv_y, yB);\r
+         tret[0] = cvB_mem.cv_tout;\r
+\r
+         return(CV_SUCCESS);\r
+       }\r
+       \r
+       /*\r
+        * CVodeGetAdjY\r
+        *\r
+        * This routine returns the interpolated forward solution at time t.\r
+        * The user must allocate space for y.\r
+        */\r
+\r
+       private int CVodeGetAdjY(CVodeMemRec cv_mem, double t, NVector y)\r
+       {\r
+         CVadjMemRec ca_mem;\r
+         int flag;\r
+\r
+         if (cv_mem == null) {\r
+           cvProcessError(null, CV_MEM_NULL, "CVODEA", "CVodeGetAdjY", MSGCV_NO_MEM);\r
+           return(CV_MEM_NULL);\r
+         }\r
+\r
+         ca_mem = cv_mem.cv_adj_mem;\r
+\r
+         if (ca_mem.ca_IMget == CVAhermiteGetY_select) {\r
+               flag = CVAhermiteGetY(cv_mem, t, y, null);\r
+         }\r
+         else {\r
+                   flag = CVApolynomialGetY(cv_mem, t, y, null);\r
+         }\r
+         \r
+         return(flag);\r
+       }\r
+       \r
+       /*\r
+        * CVodeGetAdjCVodeBmem\r
+        *\r
+        * This function returns a (void *) pointer to the CVODES     \r
+        * memory allocated for the backward problem. This pointer can    \r
+        * then be used to call any of the CVodeGet* CVODES routines to  \r
+        * extract optional output for the backward integration phase.\r
+        */\r
+\r
+       private CVodeMemRec CVodeGetAdjCVodeBmem(CVodeMemRec cv_mem, int which)\r
+       {\r
+         CVadjMemRec ca_mem;\r
+         CVodeBMemRec cvB_mem;\r
+         CVodeMemRec cvodeB_mem;\r
+\r
+         /* Check if cvode_mem exists */\r
+         if (cv_mem == null) {\r
+           cvProcessError(null, 0, "CVODEA", "CVodeGetAdjCVodeBmem", MSGCV_NO_MEM);\r
+           return(null);\r
+         }\r
+\r
+         /* Was ASA initialized? */\r
+         if (cv_mem.cv_adjMallocDone == false) {\r
+           cvProcessError(cv_mem, 0, "CVODEA", "CVodeGetAdjCVodeBmem", MSGCV_NO_ADJ);\r
+           return(null);\r
+         } \r
+         ca_mem = cv_mem.cv_adj_mem;\r
+\r
+         /* Check which */\r
+         if ( which >= ca_mem.ca_nbckpbs ) {\r
+           cvProcessError(cv_mem, 0, "CVODEA", "CVodeGetAdjCVodeBmem", MSGCV_BAD_WHICH);\r
+           return(null);\r
+         }\r
+\r
+         /* Find the CVodeBMem entry in the linked list corresponding to which */\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
+         return(cvodeB_mem);\r
+       }\r
+\r
+       /*\r
+        * CVodeGetNumSteps\r
+        *\r
+        * Returns the current number of integration steps\r
+        */\r
+\r
+       private int CVodeGetNumSteps(CVodeMemRec cv_mem, long nsteps[])\r
+       {\r
+\r
+         if (cv_mem== null) {\r
+           cvProcessError(null, CV_MEM_NULL, "CVODES", "CVodeGetNumSteps", MSGCV_NO_MEM);\r
+           return(CV_MEM_NULL);\r
+         }\r
+\r
+         nsteps[0] = cv_mem.cv_nst;\r
+\r
+         return(CV_SUCCESS);\r
+       }\r
+\r
+       /*\r
+        * CVodeGetQuadB\r
+        */\r
+\r
+       private int CVodeGetQuadB(CVodeMemRec cv_mem, int which, double tret[], NVector qB)\r
+       {\r
+         CVadjMemRec ca_mem;\r
+         CVodeBMemRec cvB_mem;\r
+         CVodeMemRec cvodeB_mem;\r
+         long nstB[] = new long[1];\r
+         int flag;\r
+\r
+         /* Check if cvode_mem exists */\r
+         if (cv_mem == null) {\r
+           cvProcessError(null, CV_MEM_NULL, "CVODEA", "CVodeGetQuadB", 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", "CVodeGetQuadB", MSGCV_NO_ADJ);\r
+           return(CV_NO_ADJ);\r
+         } \r
+\r
+         ca_mem = cv_mem.cv_adj_mem;\r
+\r
+         /* Check the value of which */\r
+         if ( which >= ca_mem.ca_nbckpbs ) {\r
+           cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeGetQuadB", MSGCV_BAD_WHICH);\r
+           return(CV_ILL_INPUT);\r
+         }\r
+\r
+         /* Find the CVodeBMem entry in the linked list corresponding to which */\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
+         /* If the integration for this backward problem has not started yet,\r
+          * simply return the current value of qB (i.e. the final conditions) */\r
+\r
+         flag = CVodeGetNumSteps(cvodeB_mem, nstB);\r
+         \r
+         if (nstB[0] == 0) {\r
+           N_VScale_Serial(ONE, cvB_mem.cv_mem.cv_znQ[0], qB);\r
+           tret[0] = cvB_mem.cv_tout;\r
+         } else {\r
+           flag = CVodeGetQuad(cvodeB_mem, tret, qB);\r
+         }\r
+\r
+         return(flag);\r
+       }\r
+       \r
+               private int CVodeReInitB(CVodeMemRec cv_mem, int which,\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
+       if (cv_mem == null) {\r
+       cvProcessError(null, CV_MEM_NULL, "CVODEA", "CVodeReInitB", 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", "CVodeReInitB", 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
+       if ( which >= ca_mem.ca_nbckpbs ) {\r
+       cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeReInitB", MSGCV_BAD_WHICH);\r
+       return(CV_ILL_INPUT);\r
+       }\r
+       \r
+       /* Find the CVodeBMem entry in the linked list corresponding to which */\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
+       /* Reinitialize CVODES object */\r
+       \r
+       flag = CVodeReInit(cvodeB_mem, tB0, yB0);\r
+       \r
+       return(flag);\r
+       }\r
+\r
+       private int CVodeQuadReInitB(CVodeMemRec cv_mem, int which, NVector yQB0)\r
+       {\r
+         CVadjMemRec ca_mem;\r
+         CVodeBMemRec cvB_mem;\r
+         CVodeMemRec cvodeB_mem;\r
+         int flag;\r
+\r
+         /* Check if cvode_mem exists */\r
+         if (cv_mem == null) {\r
+           cvProcessError(null, CV_MEM_NULL, "CVODEA", "CVodeQuadReInitB", 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", "CVodeQuadReInitB", 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
+         if ( which >= ca_mem.ca_nbckpbs ) {\r
+           cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeQuadReInitB", MSGCV_BAD_WHICH);\r
+           return(CV_ILL_INPUT);\r
+         }\r
+\r
+         /* Find the CVodeBMem entry in the linked list corresponding to which */\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
+         flag = CVodeQuadReInit(cvodeB_mem, yQB0);\r
+         if (flag != CV_SUCCESS) return(flag);\r
+\r
+         return(CV_SUCCESS);\r
+       }\r
 \r
 }
\ No newline at end of file