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

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

index e3bad28..a899409 100644 (file)
@@ -372,6 +372,33 @@ public abstract class CVODES {
        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
+       /*-----------------------------------------------------------------\r
+        * IV. SUNLinearSolver error codes\r
+        * ---------------------------------------------------------------\r
+        */\r
+\r
+       final int SUNLS_SUCCESS = 0;  /* successful/converged          */\r
+\r
+       final int SUNLS_MEM_NULL = -1;  /* mem argument is NULL          */\r
+       final int SUNLS_ILL_INPUT = -2;  /* illegal function input        */\r
+       final int SUNLS_MEM_FAIL = -3;  /* failed memory access          */\r
+       final int SUNLS_ATIMES_FAIL_UNREC = -4;  /* atimes unrecoverable failure  */\r
+       final int SUNLS_PSET_FAIL_UNREC = -5;  /* pset unrecoverable failure    */\r
+       final int SUNLS_PSOLVE_FAIL_UNREC = -6;  /* psolve unrecoverable failure  */\r
+       final int SUNLS_PACKAGE_FAIL_UNREC = -7;  /* external package unrec. fail  */\r
+       final int SUNLS_GS_FAIL = -8;  /* Gram-Schmidt failure          */        \r
+       final int SUNLS_QRSOL_FAIL = -9;  /* QRsol found singular R        */\r
+\r
+       final int SUNLS_RES_REDUCED = 1;  /* nonconv. solve, resid reduced */\r
+       final int SUNLS_CONV_FAIL = 2;  /* nonconvergent solve           */\r
+       final int SUNLS_ATIMES_FAIL_REC = 3;  /* atimes failed recoverably     */\r
+       final int SUNLS_PSET_FAIL_REC = 4;  /* pset failed recoverably       */\r
+       final int SUNLS_PSOLVE_FAIL_REC = 5;  /* psolve failed recoverably     */\r
+       final int SUNLS_PACKAGE_FAIL_REC = 6;  /* external package recov. fail  */\r
+       final int SUNLS_QRFACT_FAIL = 7;  /* QRfact found singular matrix  */\r
+       final int SUNLS_LUFACT_FAIL = 8;  /* LUfact found singular matrix  */\r
+\r
+       \r
        public enum SUNLinearSolver_Type{\r
                  SUNLINEARSOLVER_DIRECT,\r
                  SUNLINEARSOLVER_ITERATIVE,\r
@@ -379,7 +406,17 @@ public abstract class CVODES {
                };\r
 \r
     final int cvDlsDQJac = -1;\r
-\r
+    \r
+    // For cv_mem.cv_efun_select\r
+    final int cvEwtSet_select = 1;\r
+    // For cv_mem.cv_linit_select\r
+    final int cvDlsInitialize_select = 1;\r
+    // For cv_mem.cv_lsolve_select\r
+    final int cvDlsSolve_select = 1;\r
+    // For cv_mem.cv_lfree_select\r
+    final int cvDlsFree_select = 1;\r
+    // For cv_mem.cv_setup_select\r
+    final int cvDlsSetup_select = 1;\r
 \r
     final int cvsRoberts_dns = 1;\r
     int problem = cvsRoberts_dns;\r
@@ -559,7 +596,7 @@ public abstract class CVODES {
         * @param ydotv\r
         * @return\r
         */\r
-       private int f(double t, NVector yv, NVector ydotv) {\r
+       private int f(double t, NVector yv, NVector ydotv, CVodeMemRec user_data) {\r
                double y[] = yv.getData();\r
                double ydot[] = ydotv.getData();\r
                if (problem == 1) {\r
@@ -570,6 +607,10 @@ public abstract class CVODES {
                return 0;\r
        }\r
        \r
+       private int fQ(double t, NVector x, NVector y, CVodeMemRec user_data) {\r
+               return 0;\r
+       }\r
+       \r
        /**\r
         * g routine.  Compute functions g_i(t,y) for i = 0,1.\r
         * @param t\r
@@ -653,7 +694,9 @@ public abstract class CVODES {
          NVector cv_Vabstol = new NVector();        /* vector absolute tolerance                     */\r
          boolean cv_user_efun;   /* SUNTRUE if user sets efun                     */\r
          //CVEwtFn cv_efun;            /* function to set ewt                           */\r
+         int cv_efun_select;\r
          //void *cv_e_data;            /* user pointer passed to efun                   */\r
+         CVodeMemRec cv_e_data;\r
 \r
          /*-----------------------\r
            Quadrature Related Data \r
@@ -697,7 +740,7 @@ public abstract class CVODES {
          int cv_itolS;\r
          double cv_reltolS;        /* relative tolerance for sensitivities         */\r
          double cv_SabstolS[];      /* scalar absolute tolerances for sensi.        */\r
-         NVector cv_VabstolS = new NVector();      /* vector absolute tolerances for sensi.        */\r
+         NVector cv_VabstolS[];      /* vector absolute tolerances for sensi.        */\r
 \r
          /*-----------------------------------\r
            Quadrature Sensitivity Related Data \r
@@ -714,7 +757,7 @@ public abstract class CVODES {
          int cv_itolQS;\r
          double cv_reltolQS;       /* relative tolerance for yQS                   */\r
          double cv_SabstolQS[];     /* scalar absolute tolerances for yQS           */\r
-         NVector cv_VabstolQS = new NVector();     /* vector absolute tolerances for yQS           */\r
+         NVector cv_VabstolQS[];     /* vector absolute tolerances for yQS           */\r
 \r
          /*-----------------------\r
            Nordsieck History Array \r
@@ -754,7 +797,7 @@ public abstract class CVODES {
            Sensitivity Related Vectors \r
            ---------------------------*/\r
 \r
-         NVector cv_znS[] = new NVector[L_MAX];    /* Nordsieck arrays for sensitivities           */\r
+         NVector cv_znS[][] = new NVector[L_MAX][];    /* Nordsieck arrays for sensitivities           */\r
          NVector cv_ewtS[];          /* error weight vectors for sensitivities       */\r
          NVector cv_yS[];            /* yS=yS0 (allocated by the user)               */\r
          NVector cv_acorS[];         /* acorS = yS_n(m) - yS_n(0)                    */\r
@@ -767,7 +810,7 @@ public abstract class CVODES {
            Quadrature Sensitivity Related Vectors \r
            --------------------------------------*/\r
 \r
-         NVector cv_znQS[] = new NVector[L_MAX];   /* Nordsieck arrays for quadr. sensitivities    */\r
+         NVector cv_znQS[][] = new NVector[L_MAX][];   /* Nordsieck arrays for quadr. sensitivities    */\r
          NVector cv_ewtQS[];         /* error weight vectors for sensitivities       */\r
          NVector cv_yQS[];           /* Unlike yS, yQS is not allocated by the user  */\r
          NVector cv_acorQS[];        /* acorQS = yQS_n(m) - yQS_n(0)                 */\r
@@ -902,15 +945,19 @@ public abstract class CVODES {
          /* Linear Solver functions to be called */\r
 \r
          //int (*cv_linit)(struct CVodeMemRec *cv_mem);\r
+         int cv_linit_select;\r
 \r
          //int (*cv_lsetup)(struct CVodeMemRec *cv_mem, int convfail, \r
                           //N_Vector ypred, N_Vector fpred, booleantype *jcurPtr, \r
                           //N_Vector vtemp1, N_Vector vtemp2, N_Vector vtemp3); \r
+         int cv_lsetup_select;\r
 \r
          //int (*cv_lsolve)(struct CVodeMemRec *cv_mem, N_Vector b, N_Vector weight,\r
                           //N_Vector ycur, N_Vector fcur);\r
+         int cv_lsolve_select;\r
 \r
          //int (*cv_lfree)(struct CVodeMemRec *cv_mem);\r
+         int cv_lfree_select;\r
 \r
          /* Linear Solver specific memory */\r
 \r
@@ -1073,7 +1120,8 @@ public abstract class CVODES {
          cv_mem.cv_itol       = CV_NN;\r
          cv_mem.cv_user_efun  = false;\r
          //cv_mem.cv_efun       = null;\r
-         //cv_mem.cv_e_data     = null;\r
+         cv_mem.cv_efun_select = -1;\r
+         cv_mem.cv_e_data     = null;\r
          //cv_mem.cv_ehfun      = cvErrHandler;\r
          cv_mem.cv_eh_data    = cv_mem;\r
          //cv_mem.cv_errfp      = stderr;\r
@@ -1310,9 +1358,13 @@ public abstract class CVODES {
           (We check != NULL later, in CVode, if using CV_NEWTON.) */\r
 \r
        //cv_mem.cv_linit  = null;\r
+       cv_mem.cv_linit_select = -1;\r
        //cv_mem.cv_lsetup = null;\r
+       cv_mem.cv_lsetup_select = -1;\r
        //cv_mem.cv_lsolve = null;\r
+       cv_mem.cv_lsolve_select = -1;\r
        //cv_mem.cv_lfree  = null;\r
+       cv_mem.cv_lfree_select = -1;\r
        //cv_mem.cv_lmem   = null;\r
 \r
        /* Set forceSetup to SUNFALSE */\r
@@ -1321,7 +1373,7 @@ public abstract class CVODES {
 \r
        /* Initialize zn[0] in the history array */\r
 \r
-       N_VScale(ONE, y0, cv_mem.cv_zn[0]);\r
+       N_VScale_Serial(ONE, y0, cv_mem.cv_zn[0]);\r
 \r
        /* Initialize all the counters */\r
 \r
@@ -1486,16 +1538,6 @@ public abstract class CVODES {
         x = null;\r
      }\r
      \r
-     private void N_VScale(double c, NVector x, NVector z) {\r
-        int i;\r
-        double xa[] = x.getData();\r
-        int n = xa.length;\r
-        double za[] = z.getData();\r
-        for (i = 0; i < n; i++) {\r
-                za[i] = c * xa[i];\r
-        }\r
-     }\r
-     \r
      int CVodeSVtolerances(CVodeMemRec cv_mem, double reltol, NVector abstol)\r
      {\r
 \r
@@ -1535,13 +1577,14 @@ public abstract class CVODES {
        }\r
 \r
        cv_mem.cv_reltol = reltol;\r
-       N_VScale(ONE, abstol, cv_mem.cv_Vabstol);\r
+       N_VScale_Serial(ONE, abstol, cv_mem.cv_Vabstol);\r
 \r
        cv_mem.cv_itol = CV_SV;\r
 \r
        cv_mem.cv_user_efun = false;\r
        //cv_mem.cv_efun = cvEwtSet;\r
-       //cv_mem.cv_e_data = null; /* will be set to cvode_mem in InitialSetup */\r
+       cv_mem.cv_efun_select = cvEwtSet_select;\r
+       cv_mem.cv_e_data = null; /* will be set to cvode_mem in InitialSetup */\r
 \r
        return(CV_SUCCESS);\r
      }\r
@@ -1727,8 +1770,8 @@ public abstract class CVODES {
      \r
      class SUNLinearSolver {\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
+         int pivots[]; // Array of size N, index array for partial pivoting in LU factorization\r
+         int last_flag; // last error return flag from internal setup/solve\r
          SUNLinearSolver_Type type;\r
      }\r
      \r
@@ -1747,7 +1790,7 @@ public abstract class CVODES {
         SUNLinearSolver S = new SUNLinearSolver();\r
         S.N = matrixRows;\r
         S.last_flag = 0;\r
-        S.pivots = new long[matrixRows];\r
+        S.pivots = new int[matrixRows];\r
         S.type = SUNLinearSolver_Type.SUNLINEARSOLVER_DIRECT;\r
         return S;\r
      }\r
@@ -1788,29 +1831,28 @@ public abstract class CVODES {
       //}\r
 \r
       /* free any existing system solver attached to CVode */\r
-      //if (cv_mem.cv_lfree)  cv_mem.cv_lfree(cv_mem);\r
+      if (cv_mem.cv_lfree_select > 0)  cv_lfree(cv_mem, cv_mem.cv_lfree_select);\r
 \r
       /* Set four main system linear solver function fields in cv_mem */\r
       //cv_mem.cv_linit  = cvDlsInitialize;\r
+      cv_mem.cv_linit_select = cvDlsInitialize_select;\r
       //cv_mem.cv_lsetup = cvDlsSetup;\r
+      cv_mem.cv_lsetup_select = cvDlsSetup_select;\r
       //cv_mem.cv_lsolve = cvDlsSolve;\r
+      cv_mem.cv_lsolve_select = cvDlsSolve_select;\r
       //cv_mem.cv_lfree  = cvDlsFree;\r
+      cv_mem.cv_lfree_select = cvDlsFree_select;\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.jac = cvDlsDQJac;\r
       cvdls_mem.J_data = cv_mem;\r
       cvdls_mem.last_flag = CVDLS_SUCCESS;\r
 \r
@@ -1883,7 +1925,7 @@ public abstract class CVODES {
 \r
     long nfeDQ;       /* no. of calls to f due to DQ Jacobian approx.  */\r
 \r
-    long last_flag;   /* last error return flag                        */\r
+    int last_flag;   /* last error return flag                        */\r
 \r
   };\r
 \r
@@ -2001,7 +2043,7 @@ public abstract class CVODES {
      * ----------------------------------------\r
      */\r
 \r
-   /* if (cv_mem.cv_nst == 0) {\r
+    /*if (cv_mem.cv_nst == 0) {\r
 \r
       cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tn;*/\r
 \r
@@ -2017,12 +2059,13 @@ public abstract class CVODES {
        * 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
+      /*retval = 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
+                       //MSGCV_RHSFUNC_FAILED, cv_mem.cv_tn);\r
+                         MSGCV_RHSFUNC_FAILED);\r
         return(CV_RHSFUNC_FAIL);\r
       }\r
       if (retval > 0) {\r
@@ -2031,10 +2074,10 @@ public abstract class CVODES {
         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 (cv_mem.cv_quadr) {\r
+        retval = 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
@@ -2047,7 +2090,7 @@ public abstract class CVODES {
         }\r
       }\r
 \r
-      if (cv_mem->cv_sensi) {\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
@@ -2080,11 +2123,11 @@ public abstract class CVODES {
                          "CVode", MSGCV_QSRHSFUNC_FIRST);\r
           return(CV_FIRST_QSRHSFUNC_ERR);\r
         }\r
-      }*/\r
+      } */\r
 \r
       /* Test input tstop for legality. */\r
 \r
-     /* if (cv_mem->cv_tstopset) {\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
@@ -2160,7 +2203,7 @@ public abstract class CVODES {
 \r
       }\r
 \r
-    } *//* end first call block */\r
+    }*/ /* end first call block */ // if (cv_mem.cv_nst == 0)\r
 \r
     /*\r
      * ------------------------------------------------------\r
@@ -2512,7 +2555,936 @@ public abstract class CVODES {
     return -100;\r
 \r
   }\r
+  \r
+  /*  \r
+   * cvInitialSetup\r
+   *\r
+   * This routine performs input consistency checks at the first step.\r
+   * If needed, it also checks the linear solver module and calls the\r
+   * linear solver initialization routine.\r
+   */\r
+\r
+  private int cvInitialSetup(CVodeMemRec cv_mem)\r
+  {\r
+    int ier;\r
+\r
+    /* Did the user specify tolerances? */\r
+    if (cv_mem.cv_itol == CV_NN) {\r
+      cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "cvInitialSetup",\r
+                     MSGCV_NO_TOL);\r
+      return(CV_ILL_INPUT);\r
+    }\r
+\r
+    /* Set data for efun */\r
+    if (cv_mem.cv_user_efun) cv_mem.cv_e_data = cv_mem.cv_user_data;\r
+    else                      cv_mem.cv_e_data = cv_mem;\r
+\r
+    /* Load initial error weights */\r
+    ier = cv_efun(cv_mem.cv_zn[0], cv_mem.cv_ewt,\r
+                          cv_mem.cv_e_data, cv_mem.cv_efun_select);\r
+    if (ier != 0) {\r
+      if (cv_mem.cv_itol == CV_WF) \r
+        cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "cvInitialSetup",\r
+                       MSGCV_EWT_FAIL);\r
+      else\r
+        cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "cvInitialSetup",\r
+                       MSGCV_BAD_EWT);\r
+      return(CV_ILL_INPUT);\r
+    }\r
+    \r
+    /* Quadrature initial setup */\r
+\r
+    if (cv_mem.cv_quadr && cv_mem.cv_errconQ) {\r
+\r
+      /* Did the user specify tolerances? */\r
+      if (cv_mem.cv_itolQ == CV_NN) {\r
+        cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "cvInitialSetup",\r
+                       MSGCV_NO_TOLQ);\r
+        return(CV_ILL_INPUT);\r
+      }\r
+\r
+      /* Load ewtQ */\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", "cvInitialSetup",\r
+                       MSGCV_BAD_EWTQ);\r
+        return(CV_ILL_INPUT);\r
+      }\r
+\r
+    }\r
+\r
+    if (!cv_mem.cv_quadr) cv_mem.cv_errconQ = false;\r
+\r
+    /* Forward sensitivity initial setup */\r
+\r
+    if (cv_mem.cv_sensi) {\r
+\r
+      /* Did the user specify tolerances? */\r
+      if (cv_mem.cv_itolS == CV_NN) {\r
+        cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "cvInitialSetup",\r
+                       MSGCV_NO_TOLS);\r
+        return(CV_ILL_INPUT);\r
+      }\r
+\r
+      /* If using the internal DQ functions, we must have access to the problem parameters */\r
+      if(cv_mem.cv_fSDQ && (cv_mem.cv_p == null)) {\r
+        cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "cvInitialSetup",\r
+                       MSGCV_NULL_P);\r
+        return(CV_ILL_INPUT);\r
+      }\r
+\r
+      /* Load ewtS */\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", "cvInitialSetup",\r
+                       MSGCV_BAD_EWTS);\r
+        return(CV_ILL_INPUT);\r
+      }\r
+\r
+    }\r
+\r
+    /* FSA of quadrature variables */\r
+\r
+    if (cv_mem.cv_quadr_sensi) {\r
+\r
+      /* If using the internal DQ functions, we must have access to fQ\r
+       * (i.e. quadrature integration must be enabled) and to the problem parameters */\r
+\r
+      if (cv_mem.cv_fQSDQ) {\r
+\r
+        /* Test if quadratures are defined, so we can use fQ */\r
+        if (!cv_mem.cv_quadr) {\r
+          cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "cvInitialSetup",\r
+                         MSGCV_NULL_FQ);\r
+          return(CV_ILL_INPUT);\r
+        }\r
+\r
+        /* Test if we have the problem parameters */\r
+        if(cv_mem.cv_p == null) {\r
+          cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "cvInitialSetup",\r
+                         MSGCV_NULL_P);\r
+          return(CV_ILL_INPUT);\r
+        }\r
+\r
+      }\r
+\r
+      if (cv_mem.cv_errconQS) {\r
+        \r
+        /* Did the user specify tolerances? */\r
+        if (cv_mem.cv_itolQS == CV_NN) {\r
+          cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "cvInitialSetup",\r
+                         MSGCV_NO_TOLQS);\r
+          return(CV_ILL_INPUT);\r
+        }\r
+\r
+        /* If needed, did the user provide quadrature tolerances? */\r
+        if ( (cv_mem.cv_itolQS == CV_EE) && (cv_mem.cv_itolQ == CV_NN) ) {\r
+          cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "cvInitialSetup",\r
+                         MSGCV_NO_TOLQ);\r
+          return(CV_ILL_INPUT);\r
+        }\r
+\r
+        /* Load ewtQS */\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", "cvInitialSetup",\r
+                         MSGCV_BAD_EWTQS);\r
+          return(CV_ILL_INPUT);\r
+        }\r
+\r
+      }\r
+\r
+    } else {\r
+\r
+      cv_mem.cv_errconQS = false;\r
+\r
+    }\r
+\r
+    /* Check if lsolve function exists (if needed) and call linit function (if it exists) */\r
+    if (cv_mem.cv_iter == CV_NEWTON) {\r
+      if (cv_mem.cv_lsolve_select < 0) {\r
+        cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "cvInitialSetup",\r
+                       MSGCV_LSOLVE_NULL);\r
+        return(CV_ILL_INPUT);\r
+      }\r
+      if (cv_mem.cv_linit_select > 0) {\r
+        ier = cv_linit(cv_mem, cv_mem.cv_linit_select);\r
+        if (ier != 0) {\r
+          cvProcessError(cv_mem, CV_LINIT_FAIL, "CVODES", "cvInitialSetup",\r
+                         MSGCV_LINIT_FAIL);\r
+          return(CV_LINIT_FAIL);\r
+        }\r
+      }\r
+    }\r
+      \r
+    return(CV_SUCCESS);\r
+  }\r
+\r
+  /*\r
+   * cvEwtSet\r
+   *\r
+   * This routine is responsible for setting the error weight vector ewt,\r
+   * according to tol_type, as follows:\r
+   *\r
+   * (1) ewt[i] = 1 / (reltol * SUNRabs(ycur[i]) + *abstol), i=0,...,neq-1\r
+   *     if tol_type = CV_SS\r
+   * (2) ewt[i] = 1 / (reltol * SUNRabs(ycur[i]) + abstol[i]), i=0,...,neq-1\r
+   *     if tol_type = CV_SV\r
+   *\r
+   * cvEwtSet returns 0 if ewt is successfully set as above to a\r
+   * positive vector and -1 otherwise. In the latter case, ewt is\r
+   * considered undefined.\r
+   *\r
+   * All the real work is done in the routines cvEwtSetSS, cvEwtSetSV.\r
+   */\r
+\r
+  int cvEwtSet(NVector ycur, NVector weight, CVodeMemRec cv_mem)\r
+  {\r
+    int flag = 0;\r
+\r
+    switch(cv_mem.cv_itol) {\r
+    case CV_SS:\r
+      flag = cvEwtSetSS(cv_mem, ycur, weight);\r
+      break;\r
+    case CV_SV:\r
+      flag = cvEwtSetSV(cv_mem, ycur, weight);\r
+      break;\r
+    }\r
+\r
+    return(flag);\r
+  }\r
+  \r
+  /*\r
+   * cvEwtSetSS\r
+   *\r
+   * This routine sets ewt as decribed above in the case tol_type = CV_SS.\r
+   * It tests for non-positive components before inverting. cvEwtSetSS\r
+   * returns 0 if ewt is successfully set to a positive vector\r
+   * and -1 otherwise. In the latter case, ewt is considered undefined.\r
+   */\r
+\r
+  private int cvEwtSetSS(CVodeMemRec cv_mem, NVector ycur, NVector weight)\r
+  {\r
+    N_VAbs_Serial(ycur, cv_mem.cv_tempv);\r
+    N_VScale_Serial(cv_mem.cv_reltol, cv_mem.cv_tempv, cv_mem.cv_tempv);\r
+    N_VAddConst_Serial(cv_mem.cv_tempv, cv_mem.cv_Sabstol, cv_mem.cv_tempv);\r
+    if (N_VMin(cv_mem.cv_tempv) <= ZERO) return(-1);\r
+    N_VInv_Serial(cv_mem.cv_tempv, weight);\r
+\r
+    return(0);\r
+  }\r
+\r
+  /*\r
+   * cvEwtSetSV\r
+   *\r
+   * This routine sets ewt as decribed above in the case tol_type = CV_SV.\r
+   * It tests for non-positive components before inverting. cvEwtSetSV\r
+   * returns 0 if ewt is successfully set to a positive vector\r
+   * and -1 otherwise. In the latter case, ewt is considered undefined.\r
+   */\r
+\r
+  private int cvEwtSetSV(CVodeMemRec cv_mem, NVector ycur, NVector weight)\r
+  {\r
+    N_VAbs_Serial(ycur, cv_mem.cv_tempv);\r
+    N_VLinearSum_Serial(cv_mem.cv_reltol, cv_mem.cv_tempv, ONE,\r
+                 cv_mem.cv_Vabstol, cv_mem.cv_tempv);\r
+    if (N_VMin(cv_mem.cv_tempv) <= ZERO) return(-1);\r
+    N_VInv_Serial(cv_mem.cv_tempv, weight);\r
+    return(0);\r
+  }\r
+\r
+  private void N_VAbs_Serial(NVector x, NVector z)\r
+  {\r
+    int i, N;\r
+    double xd[], zd[];\r
+\r
+    xd = zd = null;\r
+    xd = x.data;\r
+    zd = z.data;\r
+    N = xd.length;\r
+\r
+    for (i = 0; i < N; i++)\r
+      zd[i] = Math.abs(xd[i]);\r
+\r
+    return;\r
+  }\r
+  \r
+  private void N_VAddConst_Serial(NVector x, double b, NVector z)\r
+  {\r
+    int i, N;\r
+    double xd[], zd[];\r
+\r
+    xd = zd = null;\r
+    xd = x.data;\r
+    zd = z.data;\r
+    N = xd.length;\r
+\r
+    for (i = 0; i < N; i++) \r
+      zd[i] = xd[i]+b;\r
+\r
+    return;\r
+  }\r
+\r
+  private void N_VInv_Serial(NVector x, NVector z)\r
+  {\r
+    int i, N;\r
+    double xd[], zd[];\r
+\r
+    xd = zd = null;\r
+    \r
+    xd = x.data;\r
+    zd = z.data;\r
+    N = xd.length;\r
+\r
+    for (i = 0; i < N; i++)\r
+      zd[i] = ONE/xd[i];\r
+\r
+    return;\r
+  }\r
+  \r
+  private void N_VLinearSum_Serial(double a, NVector x, double b, NVector y, NVector z)\r
+  {\r
+    int i, N;\r
+    double c, xd[], yd[], zd[];\r
+    NVector v1, v2;\r
+    boolean test;\r
+\r
+    xd = yd = zd = null;\r
+\r
+    if ((b == ONE) && (z == y)) {    /* BLAS usage: axpy y <- ax+y */\r
+      Vaxpy_Serial(a,x,y);\r
+      return;\r
+    }\r
+\r
+    if ((a == ONE) && (z == x)) {    /* BLAS usage: axpy x <- by+x */\r
+      Vaxpy_Serial(b,y,x);\r
+      return;\r
+    }\r
+\r
+    /* Case: a == b == 1.0 */\r
+\r
+    if ((a == ONE) && (b == ONE)) {\r
+      VSum_Serial(x, y, z);\r
+      return;\r
+    }\r
+\r
+    /* Cases: (1) a == 1.0, b = -1.0, (2) a == -1.0, b == 1.0 */\r
+\r
+    if ((test = ((a == ONE) && (b == -ONE))) || ((a == -ONE) && (b == ONE))) {\r
+      v1 = test ? y : x;\r
+      v2 = test ? x : y;\r
+      VDiff_Serial(v2, v1, z);\r
+      return;\r
+    }\r
+\r
+    /* Cases: (1) a == 1.0, b == other or 0.0, (2) a == other or 0.0, b == 1.0 */\r
+    /* if a or b is 0.0, then user should have called N_VScale */\r
+\r
+    if ((test = (a == ONE)) || (b == ONE)) {\r
+      c  = test ? b : a;\r
+      v1 = test ? y : x;\r
+      v2 = test ? x : y;\r
+      VLin1_Serial(c, v1, v2, z);\r
+      return;\r
+    }\r
+\r
+    /* Cases: (1) a == -1.0, b != 1.0, (2) a != 1.0, b == -1.0 */\r
+\r
+    if ((test = (a == -ONE)) || (b == -ONE)) {\r
+      c = test ? b : a;\r
+      v1 = test ? y : x;\r
+      v2 = test ? x : y;\r
+      VLin2_Serial(c, v1, v2, z);\r
+      return;\r
+    }\r
+\r
+    /* Case: a == b */\r
+    /* catches case both a and b are 0.0 - user should have called N_VConst */\r
+\r
+    if (a == b) {\r
+      VScaleSum_Serial(a, x, y, z);\r
+      return;\r
+    }\r
+\r
+    /* Case: a == -b */\r
+\r
+    if (a == -b) {\r
+      VScaleDiff_Serial(a, x, y, z);\r
+      return;\r
+    }\r
+\r
+    /* Do all cases not handled above:\r
+       (1) a == other, b == 0.0 - user should have called N_VScale\r
+       (2) a == 0.0, b == other - user should have called N_VScale\r
+       (3) a,b == other, a !=b, a != -b */\r
+    \r
+    xd = x.data;\r
+    yd = y.data;\r
+    zd = z.data;\r
+    N = xd.length;\r
+\r
+    for (i = 0; i < N; i++)\r
+      zd[i] = (a*xd[i])+(b*yd[i]);\r
+\r
+    return;\r
+  }\r
+  \r
+  private void Vaxpy_Serial(double a, NVector x, NVector y)\r
+  {\r
+    int i, N;\r
+    double xd[], yd[];\r
+\r
+    xd = yd = null;\r
+\r
+    xd = x.data;\r
+    yd = y.data;\r
+    N = xd.length;\r
+\r
+    if (a == ONE) {\r
+      for (i = 0; i < N; i++)\r
+        yd[i] += xd[i];\r
+      return;\r
+    }\r
+\r
+    if (a == -ONE) {\r
+      for (i = 0; i < N; i++)\r
+        yd[i] -= xd[i];\r
+      return;\r
+    }    \r
+\r
+    for (i = 0; i < N; i++)\r
+      yd[i] += a*xd[i];\r
+\r
+    return;\r
+  }\r
+\r
+\r
+  private void N_VScale_Serial(double c, NVector x, NVector z)\r
+  {\r
+    int i, N;\r
+    double xd[], zd[];\r
+\r
+    xd = zd = null;\r
+\r
+    if (z == x) {  /* BLAS usage: scale x <- cx */\r
+      VScaleBy_Serial(c, x);\r
+      return;\r
+    }\r
+\r
+    if (c == ONE) {\r
+      VCopy_Serial(x, z);\r
+    } else if (c == -ONE) {\r
+      VNeg_Serial(x, z);\r
+    } else {\r
+       xd = x.data;\r
+        zd = z.data;\r
+        N = xd.length;\r
+      for (i = 0; i < N; i++) \r
+        zd[i] = c*xd[i];\r
+    }\r
+\r
+    return;\r
+  }\r
+  \r
+  private void VScaleBy_Serial(double a, NVector x)\r
+  {\r
+    int i, N;\r
+    double xd[];\r
+\r
+    xd = null;\r
+\r
+    xd = x.data;\r
+    N = xd.length;\r
+\r
+    for (i = 0; i < N; i++)\r
+      xd[i] *= a;\r
+\r
+    return;\r
+  }\r
+  \r
+  private void VCopy_Serial(NVector x, NVector z)\r
+  {\r
+    int i, N;\r
+    double xd[], zd[];\r
+\r
+    xd = zd = null;\r
+\r
+    xd = x.data;\r
+    zd = z.data;\r
+    N = xd.length;\r
+\r
+    for (i = 0; i < N; i++)\r
+      zd[i] = xd[i]; \r
+\r
+    return;\r
+  }\r
+\r
+  private void VNeg_Serial(NVector x, NVector z)\r
+  {\r
+    int i, N;\r
+    double xd[], zd[];\r
+\r
+    xd = zd = null;\r
+\r
+    xd = x.data;\r
+    zd = z.data;\r
+    N = xd.length;\r
+    \r
+    for (i = 0; i < N; i++)\r
+      zd[i] = -xd[i];\r
+\r
+    return;\r
+  }\r
+  \r
+  private void VSum_Serial(NVector x, NVector y, NVector z)\r
+  {\r
+    int i, N;\r
+    double xd[], yd[], zd[];\r
+\r
+    xd = yd = zd = null;\r
+\r
+    xd = x.data;\r
+    yd = y.data;\r
+    zd = z.data;\r
+    N = xd.length;\r
+\r
+    for (i = 0; i < N; i++)\r
+      zd[i] = xd[i]+yd[i];\r
+\r
+    return;\r
+  }\r
+\r
+  private void VDiff_Serial(NVector x, NVector y, NVector z)\r
+  {\r
+    int i, N;\r
+    double xd[], yd[], zd[];\r
+\r
+    xd = yd = zd = null;\r
+\r
+    xd = x.data;\r
+    yd = y.data;\r
+    zd = z.data;\r
+    N = xd.length;\r
+\r
+    for (i = 0; i < N; i++)\r
+      zd[i] = xd[i]-yd[i];\r
+\r
+    return;\r
+  }\r
+\r
+  private void VLin1_Serial(double a, NVector x, NVector y, NVector z)\r
+  {\r
+    int i, N;\r
+    double xd[], yd[], zd[];\r
+\r
+    xd = yd = zd = null;\r
+\r
+    xd = x.data;\r
+    yd = y.data;\r
+    zd = z.data;\r
+    N = xd.length;\r
+\r
+    for (i = 0; i < N; i++)\r
+      zd[i] = (a*xd[i])+yd[i];\r
+\r
+    return;\r
+  }\r
+\r
+  private void VLin2_Serial(double a, NVector x, NVector y, NVector z)\r
+  {\r
+    int i, N;\r
+    double xd[], yd[], zd[];\r
+\r
+    xd = yd = zd = null;\r
+\r
+    xd = x.data;\r
+    yd = y.data;\r
+    zd = z.data;\r
+    N = xd.length;\r
+\r
+    for (i = 0; i < N; i++)\r
+      zd[i] = (a*xd[i])-yd[i];\r
+\r
+    return;\r
+  }\r
+  \r
+  private void VScaleSum_Serial(double c, NVector x, NVector y, NVector z)\r
+  {\r
+    int i, N;\r
+    double xd[], yd[], zd[];\r
+\r
+    xd = yd = zd = null;\r
+\r
+    xd = x.data;\r
+    yd = y.data;\r
+    zd = z.data;\r
+    N = xd.length;\r
+\r
+    for (i = 0; i < N; i++)\r
+      zd[i] = c*(xd[i]+yd[i]);\r
+\r
+    return;\r
+  }\r
+\r
+  private void VScaleDiff_Serial(double c, NVector x, NVector y, NVector z)\r
+  {\r
+    int i, N;\r
+    double xd[], yd[], zd[];\r
+\r
+    xd = yd = zd = null;\r
+\r
+    xd = x.data;\r
+    yd = y.data;\r
+    zd = z.data;\r
+    N = xd.length;\r
+\r
+    for (i = 0; i < N; i++)\r
+      zd[i] = c*(xd[i]-yd[i]);\r
+\r
+    return;\r
+  }\r
+\r
+  private int cv_efun(NVector ycur, NVector weight, CVodeMemRec cv_mem, int cv_efun_select) {\r
+         int flag = 0;\r
+         switch (cv_efun_select) {\r
+         case cvEwtSet_select:\r
+                 flag = cvEwtSet(ycur, weight, cv_mem);  \r
+                 break;\r
+         }\r
+         return flag;\r
+  }\r
+  \r
+  /*\r
+   * cvQuadEwtSet\r
+   *\r
+   */\r
+\r
+  private int cvQuadEwtSet(CVodeMemRec cv_mem, NVector qcur, NVector weightQ)\r
+  {\r
+    int flag=0;\r
+\r
+    switch (cv_mem.cv_itolQ) {\r
+    case CV_SS:\r
+      flag = cvQuadEwtSetSS(cv_mem, qcur, weightQ);\r
+      break;\r
+    case CV_SV:\r
+      flag = cvQuadEwtSetSV(cv_mem, qcur, weightQ);\r
+      break;\r
+    }\r
+\r
+    return(flag);\r
+\r
+  }\r
+\r
+  /*\r
+   * cvQuadEwtSetSS\r
+   *\r
+   */\r
+\r
+  private int cvQuadEwtSetSS(CVodeMemRec cv_mem, NVector qcur, NVector weightQ)\r
+  {\r
+    N_VAbs_Serial(qcur, cv_mem.cv_tempvQ);\r
+    N_VScale_Serial(cv_mem.cv_reltolQ, cv_mem.cv_tempvQ, cv_mem.cv_tempvQ);\r
+    N_VAddConst_Serial(cv_mem.cv_tempvQ, cv_mem.cv_SabstolQ, cv_mem.cv_tempvQ);\r
+    if (N_VMin(cv_mem.cv_tempvQ) <= ZERO) return(-1);\r
+    N_VInv_Serial(cv_mem.cv_tempvQ, weightQ);\r
+\r
+    return(0);\r
+  }\r
+\r
+  /*\r
+   * cvQuadEwtSetSV\r
+   *\r
+   */\r
+\r
+  private int cvQuadEwtSetSV(CVodeMemRec cv_mem, NVector qcur, NVector weightQ)\r
+  {\r
+    N_VAbs_Serial(qcur, cv_mem.cv_tempvQ);\r
+    N_VLinearSum_Serial(cv_mem.cv_reltolQ, cv_mem.cv_tempvQ, ONE,\r
+                 cv_mem.cv_VabstolQ, cv_mem.cv_tempvQ);\r
+    if (N_VMin(cv_mem.cv_tempvQ) <= ZERO) return(-1);\r
+    N_VInv_Serial(cv_mem.cv_tempvQ, weightQ);\r
+\r
+    return(0);\r
+  }\r
+  \r
+  /*\r
+   * cvSensEwtSet\r
+   *\r
+   */\r
+\r
+  private int cvSensEwtSet(CVodeMemRec cv_mem, NVector yScur[], NVector weightS[])\r
+  {\r
+    int flag=0;\r
+\r
+    switch (cv_mem.cv_itolS) {\r
+    case CV_EE:\r
+      flag = cvSensEwtSetEE(cv_mem, yScur, weightS);\r
+      break;\r
+    case CV_SS:\r
+      flag = cvSensEwtSetSS(cv_mem, yScur, weightS);\r
+      break;\r
+    case CV_SV:\r
+      flag = cvSensEwtSetSV(cv_mem, yScur, weightS);\r
+      break;\r
+    }\r
+\r
+    return(flag);\r
+  }\r
+\r
+  /*\r
+   * cvSensEwtSetEE\r
+   *\r
+   * In this case, the error weight vector for the i-th sensitivity is set to\r
+   *\r
+   * ewtS_i = pbar_i * efun(pbar_i*yS_i)\r
+   *\r
+   * In other words, the scaled sensitivity pbar_i * yS_i has the same error\r
+   * weight vector calculation as the solution vector.\r
+   *\r
+   */\r
+\r
+  private int cvSensEwtSetEE(CVodeMemRec cv_mem, NVector yScur[], NVector weightS[])\r
+  {\r
+    int is;\r
+    NVector pyS;\r
+    int flag;\r
+\r
+    /* Use tempvS[0] as temporary storage for the scaled sensitivity */\r
+    pyS = cv_mem.cv_tempvS[0];\r
+\r
+    for (is=0; is<cv_mem.cv_Ns; is++) {\r
+      N_VScale_Serial(cv_mem.cv_pbar[is], yScur[is], pyS);\r
+      flag = cv_efun(pyS, weightS[is], cv_mem.cv_e_data, cv_mem.cv_efun_select);\r
+      if (flag != 0) return(-1);\r
+      N_VScale_Serial(cv_mem.cv_pbar[is], weightS[is], weightS[is]);\r
+    }\r
+\r
+    return(0);\r
+  }\r
+\r
+  /*\r
+   * cvSensEwtSetSS\r
+   *\r
+   */\r
+\r
+  private int cvSensEwtSetSS(CVodeMemRec cv_mem, NVector yScur[], NVector weightS[])\r
+  {\r
+    int is;\r
+    \r
+    for (is=0; is<cv_mem.cv_Ns; is++) {\r
+      N_VAbs_Serial(yScur[is], cv_mem.cv_tempv);\r
+      N_VScale_Serial(cv_mem.cv_reltolS, cv_mem.cv_tempv, cv_mem.cv_tempv);\r
+      N_VAddConst_Serial(cv_mem.cv_tempv, cv_mem.cv_SabstolS[is], cv_mem.cv_tempv);\r
+      if (N_VMin(cv_mem.cv_tempv) <= ZERO) return(-1);\r
+      N_VInv_Serial(cv_mem.cv_tempv, weightS[is]);\r
+    }\r
+    return(0);\r
+  }\r
+  \r
+  /*\r
+   * cvSensEwtSetSV\r
+   *\r
+   */\r
+\r
+  private int cvSensEwtSetSV(CVodeMemRec cv_mem, NVector yScur[], NVector weightS[])\r
+  {\r
+    int is;\r
+    \r
+    for (is=0; is<cv_mem.cv_Ns; is++) {\r
+      N_VAbs_Serial(yScur[is], cv_mem.cv_tempv);\r
+      N_VLinearSum_Serial(cv_mem.cv_reltolS, cv_mem.cv_tempv, ONE,\r
+                   cv_mem.cv_VabstolS[is], cv_mem.cv_tempv);\r
+      if (N_VMin(cv_mem.cv_tempv) <= ZERO) return(-1);\r
+      N_VInv_Serial(cv_mem.cv_tempv, weightS[is]);\r
+    }\r
+\r
+    return(0);\r
+  }\r
+\r
+  /*\r
+   * cvQuadSensEwtSet\r
+   *\r
+   */\r
+\r
+  private int cvQuadSensEwtSet(CVodeMemRec cv_mem, NVector yQScur[], NVector weightQS[])\r
+  {\r
+    int flag=0;\r
+\r
+    switch (cv_mem.cv_itolQS) {\r
+    case CV_EE:\r
+      flag = cvQuadSensEwtSetEE(cv_mem, yQScur, weightQS);\r
+      break;\r
+    case CV_SS:\r
+      flag = cvQuadSensEwtSetSS(cv_mem, yQScur, weightQS);\r
+      break;\r
+    case CV_SV:\r
+      flag = cvQuadSensEwtSetSV(cv_mem, yQScur, weightQS);\r
+      break;\r
+    }\r
+\r
+    return(flag);\r
+  }\r
+\r
+  /*\r
+   * cvQuadSensEwtSetEE\r
+   *\r
+   * In this case, the error weight vector for the i-th quadrature sensitivity\r
+   * is set to\r
+   *\r
+   * ewtQS_i = pbar_i * cvQuadEwtSet(pbar_i*yQS_i)\r
+   *\r
+   * In other words, the scaled sensitivity pbar_i * yQS_i has the same error\r
+   * weight vector calculation as the quadrature vector.\r
+   *\r
+   */\r
+  private int cvQuadSensEwtSetEE(CVodeMemRec cv_mem, NVector yQScur[], NVector weightQS[])\r
+  {\r
+    int is;\r
+    NVector pyS;\r
+    int flag;\r
+\r
+    /* Use tempvQS[0] as temporary storage for the scaled sensitivity */\r
+    pyS = cv_mem.cv_tempvQS[0];\r
+\r
+    for (is=0; is<cv_mem.cv_Ns; is++) {\r
+      N_VScale_Serial(cv_mem.cv_pbar[is], yQScur[is], pyS);\r
+      flag = cvQuadEwtSet(cv_mem, pyS, weightQS[is]);\r
+      if (flag != 0) return(-1);\r
+      N_VScale_Serial(cv_mem.cv_pbar[is], weightQS[is], weightQS[is]);\r
+    }\r
+\r
+    return(0);\r
+  }\r
+\r
+  private int cvQuadSensEwtSetSS(CVodeMemRec cv_mem, NVector yQScur[], NVector weightQS[])\r
+  {\r
+    int is;\r
+\r
+    for (is=0; is<cv_mem.cv_Ns; is++) {\r
+      N_VAbs_Serial(yQScur[is], cv_mem.cv_tempvQ);\r
+      N_VScale_Serial(cv_mem.cv_reltolQS, cv_mem.cv_tempvQ, cv_mem.cv_tempvQ);\r
+      N_VAddConst_Serial(cv_mem.cv_tempvQ, cv_mem.cv_SabstolQS[is], cv_mem.cv_tempvQ);\r
+      if (N_VMin(cv_mem.cv_tempvQ) <= ZERO) return(-1);\r
+      N_VInv_Serial(cv_mem.cv_tempvQ, weightQS[is]);\r
+    }\r
+\r
+    return(0);\r
+  }\r
+\r
+  private int cvQuadSensEwtSetSV(CVodeMemRec cv_mem, NVector yQScur[], NVector weightQS[])\r
+  {\r
+    int is;\r
+    \r
+    for (is=0; is<cv_mem.cv_Ns; is++) {\r
+      N_VAbs_Serial(yQScur[is], cv_mem.cv_tempvQ);\r
+      N_VLinearSum_Serial(cv_mem.cv_reltolQS, cv_mem.cv_tempvQ, ONE,\r
+                   cv_mem.cv_VabstolQS[is], cv_mem.cv_tempvQ);\r
+      if (N_VMin(cv_mem.cv_tempvQ) <= ZERO) return(-1);\r
+      N_VInv_Serial(cv_mem.cv_tempvQ, weightQS[is]);\r
+    }\r
+\r
+    return(0);\r
+  }\r
+  \r
+  private int cv_linit(CVodeMemRec cv_mem, int cv_linit_select) {\r
+         int flag = 0;\r
+         switch (cv_linit_select) {\r
+         case cvDlsInitialize_select:\r
+                 flag = cvDlsInitialize(cv_mem);\r
+                 break;\r
+         }\r
+         return flag;\r
+  }\r
+  \r
+  /*-----------------------------------------------------------------\r
+  cvDlsInitialize\r
+  -----------------------------------------------------------------\r
+  This routine performs remaining initializations specific\r
+  to the direct linear solver interface (and solver itself)\r
+  -----------------------------------------------------------------*/\r
+    private int cvDlsInitialize(CVodeMemRec cv_mem)\r
+{\r
+  CVDlsMemRec cvdls_mem;\r
+\r
+  /* Return immediately if cv_mem or cv_mem->cv_lmem are NULL */\r
+  if (cv_mem == null) {\r
+    cvProcessError(null, CVDLS_MEM_NULL, "CVSDLS", \r
+                    "cvDlsInitialize", MSGD_CVMEM_NULL);\r
+    return(CVDLS_MEM_NULL);\r
+  }\r
+  if (cv_mem.cv_lmem == null) {\r
+    cvProcessError(cv_mem, CVDLS_LMEM_NULL, "CVSDLS", \r
+                    "cvDlsInitialize", MSGD_LMEM_NULL);\r
+    return(CVDLS_LMEM_NULL);\r
+  }\r
+  cvdls_mem = cv_mem.cv_lmem;\r
\r
+  cvDlsInitializeCounters(cvdls_mem);\r
+\r
+  /* Set Jacobian function and data, depending on jacDQ (in case \r
+     it has changed based on user input) */\r
+  if (cvdls_mem.jacDQ) {\r
+    cvdls_mem.jac    = cvDlsDQJac;\r
+    cvdls_mem.J_data = cv_mem;\r
+  } else {\r
+    cvdls_mem.J_data = cv_mem.cv_user_data;\r
+  }\r
+\r
+  /* Call LS initialize routine */\r
+  cvdls_mem.last_flag = SUNLinSolInitialize_Dense(cvdls_mem.LS);\r
+  return(cvdls_mem.last_flag);\r
+}\r
+    \r
+    private int SUNLinSolInitialize_Dense(SUNLinearSolver S)\r
+    {\r
+      /* all solver-specific memory has already been allocated */\r
+      S.last_flag = SUNLS_SUCCESS;\r
+      return S.last_flag;\r
+    }\r
+\r
 \r
+   private void cv_lfree(CVodeMemRec cv_mem, int cv_lfree_select) {\r
+          switch (cv_lfree_select) {\r
+          case cvDlsFree_select:\r
+                  cvDlsFree(cv_mem);\r
+          break;\r
+          }\r
+   }\r
+   \r
+   /*-----------------------------------------------------------------\r
+   cvDlsFree\r
+   -----------------------------------------------------------------\r
+   This routine frees memory associates with the CVDls solver \r
+   interface.\r
+   -----------------------------------------------------------------*/\r
+ private int cvDlsFree(CVodeMemRec cv_mem)\r
+ {\r
+   CVDlsMemRec cvdls_mem;\r
+\r
+   /* Return immediately if cv_mem or cv_mem->cv_lmem are NULL */\r
+   if (cv_mem == null)  return (CVDLS_SUCCESS);\r
+   if (cv_mem.cv_lmem == null)  return(CVDLS_SUCCESS);\r
+   cvdls_mem = cv_mem.cv_lmem;\r
+\r
+   /* Free x vector */\r
+   if (cvdls_mem.x != null) {\r
+     N_VDestroy(cvdls_mem.x);\r
+   }\r
+\r
+   /* Free savedJ memory */\r
+   if (cvdls_mem.savedJ != null) {\r
+     cvdls_mem.savedJ = null;\r
+   }\r
+\r
+   /* Nullify other SUNMatrix pointer */\r
+   cvdls_mem.A = null;\r
+\r
+   /* free CVDls interface structure */\r
+   cv_mem.cv_lmem = null;\r
+   \r
+   return(CVDLS_SUCCESS);\r
+ }\r
 \r
 \r
 }
\ No newline at end of file