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

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

index dffd090..8b08cad 100644 (file)
@@ -125,58 +125,53 @@ public abstract class CVODES {
         */\r
 \r
        final int  CV_SUCCESS = 0;\r
-       /*#define CV_TSTOP_RETURN          1\r
-       #define CV_ROOT_RETURN           2\r
-\r
-       #define CV_WARNING              99\r
-\r
-       #define CV_TOO_MUCH_WORK        -1\r
-       #define CV_TOO_MUCH_ACC         -2\r
-       #define CV_ERR_FAILURE          -3\r
-       #define CV_CONV_FAILURE         -4\r
-\r
-       #define CV_LINIT_FAIL           -5\r
-       #define CV_LSETUP_FAIL          -6\r
-       #define CV_LSOLVE_FAIL          -7\r
-       #define CV_RHSFUNC_FAIL         -8\r
-       #define CV_FIRST_RHSFUNC_ERR    -9\r
-       #define CV_REPTD_RHSFUNC_ERR    -10\r
-       #define CV_UNREC_RHSFUNC_ERR    -11\r
-       #define CV_RTFUNC_FAIL          -12\r
-\r
-       #define CV_MEM_FAIL             -20\r
-       #define CV_MEM_NULL             -21\r
-       #define CV_ILL_INPUT            -22\r
-       #define CV_NO_MALLOC            -23\r
-       #define CV_BAD_K                -24\r
-       #define CV_BAD_T                -25\r
-       #define CV_BAD_DKY              -26\r
-       #define CV_TOO_CLOSE            -27\r
-\r
-       #define CV_NO_QUAD              -30\r
-       #define CV_QRHSFUNC_FAIL        -31\r
-       #define CV_FIRST_QRHSFUNC_ERR   -32\r
-       #define CV_REPTD_QRHSFUNC_ERR   -33\r
-       #define CV_UNREC_QRHSFUNC_ERR   -34\r
-\r
-       #define CV_NO_SENS              -40\r
-       #define CV_SRHSFUNC_FAIL        -41\r
-       #define CV_FIRST_SRHSFUNC_ERR   -42\r
-       #define CV_REPTD_SRHSFUNC_ERR   -43\r
-       #define CV_UNREC_SRHSFUNC_ERR   -44\r
-\r
-       #define CV_BAD_IS               -45\r
-\r
-       #define CV_NO_QUADSENS          -50\r
-       #define CV_QSRHSFUNC_FAIL       -51\r
-       #define CV_FIRST_QSRHSFUNC_ERR  -52\r
-       #define CV_REPTD_QSRHSFUNC_ERR  -53\r
-       #define CV_UNREC_QSRHSFUNC_ERR  -54*/\r
-\r
+       final int CV_TSTOP_RETURN = 1;\r
+       final int CV_ROOT_RETURN = 2;\r
+\r
+       final int CV_WARNING = 99;\r
+\r
+       final int CV_TOO_MUCH_WORK = -1;\r
+       final int CV_TOO_MUCH_ACC = -2;\r
+       final int CV_ERR_FAILURE = -3;\r
+       final int CV_CONV_FAILURE = -4;\r
+\r
+       final int CV_LINIT_FAIL = -5;\r
+       final int CV_LSETUP_FAIL = -6;\r
+       final int CV_LSOLVE_FAIL = -7;\r
+       final int CV_RHSFUNC_FAIL = -8;\r
+       final int CV_FIRST_RHSFUNC_ERR = -9;\r
+       final int CV_REPTD_RHSFUNC_ERR = -10;\r
+       final int CV_UNREC_RHSFUNC_ERR = -11;\r
+       final int  CV_RTFUNC_FAIL = -12;\r
+\r
+       final int CV_MEM_FAIL = -20;\r
        final int CV_MEM_NULL = -21;\r
        final int CV_ILL_INPUT = -22;\r
-\r
-\r
+       final int CV_NO_MALLOC = -23;\r
+       final int CV_BAD_K = -24;\r
+       final int CV_BAD_T = -25;\r
+       final int CV_BAD_DKY = -26;\r
+       final int CV_TOO_CLOSE = -27;\r
+\r
+       final int CV_NO_QUAD = -30;\r
+       final int CV_QRHSFUNC_FAIL = -31;\r
+       final int CV_FIRST_QRHSFUNC_ERR = -32;\r
+       final int CV_REPTD_QRHSFUNC_ERR = -33;\r
+       final int CV_UNREC_QRHSFUNC_ERR = -34;\r
+\r
+       final int CV_NO_SENS = -40;\r
+       final int CV_SRHSFUNC_FAIL = -41;\r
+       final int CV_FIRST_SRHSFUNC_ERR = -42;\r
+       final int CV_REPTD_SRHSFUNC_ERR = -43;\r
+       final int CV_UNREC_SRHSFUNC_ERR = -44;\r
+\r
+       final int CV_BAD_IS = -45;\r
+\r
+       final int CV_NO_QUADSENS = -50;\r
+       final int CV_QSRHSFUNC_FAIL = -51;\r
+       final int CV_FIRST_QSRHSFUNC_ERR = -52;\r
+       final int CV_REPTD_QSRHSFUNC_ERR = -53;\r
+       final int CV_UNREC_QSRHSFUNC_ERR = -54;\r
 \r
        final double HMIN_DEFAULT = 0.0;    /* hmin default value     */\r
        final double HMAX_INV_DEFAULT = 0.0;    /* hmax_inv default value */\r
@@ -187,6 +182,9 @@ public abstract class CVODES {
        final int MXNEF = 7;\r
        final double CORTES = 0.1;\r
        final double ZERO = 0.0;\r
+       final double ONE = 1.0;\r
+       final double ETAMX1 = 10000.0; \r
+\r
        \r
        final String MSG_TIME = "t = %lg";\r
        final String MSG_TIME_H = "t = %lg and h = %lg";\r
@@ -403,14 +401,24 @@ public abstract class CVODES {
                final double T1 = 0.4; // first output time\r
                final double TMULT = 10.0; // output time factor\r
                final int NOUT = 12; // number of output times\r
+               int f = cvsRoberts_dns;\r
+               int g = cvsRoberts_dns;\r
                \r
                // Set initial conditions\r
-               double y[] = new double[]{Y1,Y2,Y3};\r
+               NVector y = new NVector();\r
+               NVector abstol = new NVector();\r
+               double yr[] = new double[]{Y1,Y2,Y3};\r
+               N_VNew_Serial(y, NEQ);\r
+               y.setData(yr);\r
                double reltol = RTOL; // Set the scalar relative tolerance\r
                // Set the vector absolute tolerance\r
-               double abstol[] = new double[]{ATOL1,ATOL2,ATOL3};\r
+               double abstolr[] = new double[]{ATOL1,ATOL2,ATOL3};\r
+               N_VNew_Serial(abstol,NEQ);\r
+               abstol.setData(abstolr);\r
                CVodeMemRec cvode_mem;\r
                int flag;\r
+               double A[][];\r
+               SUNLinearSolver LS;\r
                \r
                // Call CVodeCreate to create the solver memory and specify the\r
                // Backward Differentiation Formula and the use of a Newton\r
@@ -419,15 +427,65 @@ public abstract class CVODES {
                if (cvode_mem == null) {\r
                    return;     \r
                }\r
-               flag = CVodeInit(cvode_mem, T0, y);\r
                \r
                // Call CVodeInit to initialize the integrator memory and specify the\r
                // user's right hand side function in y' = f(t,y), the initial time T0, and\r
-               // the initial dependent variable vector y\r
+           // the initial dependent variable vector y\r
+               flag = CVodeInit(cvode_mem, f, T0, y);\r
+               if (flag != CV_SUCCESS) {\r
+                       return;\r
+               }\r
                \r
+               // Call CVodeSVtolerances to specify the scalar relative tolerance\r
+               // and vector absolute tolerances\r
+               flag = CVodeSVtolerances(cvode_mem, reltol, abstol);\r
+               if (flag != CV_SUCCESS) {\r
+                       return;\r
+               }\r
+               \r
+               // Call CVRootInit to specify the root function g with 2 components\r
+               flag = CVodeRootInit(cvode_mem, 2, g);\r
+               if (flag != CV_SUCCESS) {\r
+                       return;\r
+               }\r
+               \r
+               // Create dense SUNMATRIX for use in linear solver\r
+               // indexed by columns stacked on top of each other\r
+               try {\r
+                   A = new double[NEQ][NEQ];\r
+               }\r
+               catch (OutOfMemoryError e) {\r
+                   MipavUtil.displayError("Out of memory error trying to create A");\r
+                   return;\r
+               }\r
+               \r
+               // Create dense SUNLinearSolver object for use by CVode\r
+               LS = SUNDenseLinearSolver(y, A);\r
+               if (LS == null) {\r
+                       return;\r
+               }\r
+               \r
+               // Call CVDlsSetLinearSolver to attach the matrix and linear solver to CVode\r
+               //flag = CVDlsSetLinearSolver(cvode_mem, LS, A);\r
        }\r
        \r
-       private int f(double t, double y[], double ydot[]) {\r
+       private void N_VNew_Serial(NVector y, int length) {\r
+           if ((y != null) && (length > 0)) {\r
+               y.data = new double[length];\r
+               y.own_data = true;\r
+           }\r
+       }\r
+       \r
+       /**\r
+        * f routine.  Compte function f(t,y).\r
+        * @param t\r
+        * @param yv\r
+        * @param ydotv\r
+        * @return\r
+        */\r
+       private int f(double t, NVector yv, NVector ydotv) {\r
+               double y[] = yv.getData();\r
+               double ydot[] = ydotv.getData();\r
                if (problem == 1) {\r
                    ydot[0] = -0.04*y[0] + 1.0E4*y[1]*y[2];\r
                    ydot[2] = 3.0E7*y[1]*y[1];\r
@@ -436,11 +494,39 @@ public abstract class CVODES {
                return 0;\r
        }\r
        \r
+       /**\r
+        * g routine.  Compute functions g_i(t,y) for i = 0,1.\r
+        * @param t\r
+        * @param yv\r
+        * @param gout\r
+        * @return\r
+        */\r
+       private int g(double t, NVector yv, double gout[]) {\r
+               double y[] = yv.getData();\r
+               if (problem == 1) {\r
+                   gout[0] = y[0] - 0.0001;\r
+                   gout[1] = y[2] - 0.01;\r
+               }\r
+               return 0;\r
+       }\r
+       \r
        // Types: struct CVodeMemRec, CVodeMem\r
        // -----------------------------------------------------------------\r
        // The type CVodeMem is type pointer to struct CVodeMemRec.\r
        // This structure contains fields to keep track of problem state.\r
        \r
+       private class NVector {\r
+               double data[];\r
+               boolean own_data;\r
+               \r
+               void setData(double data[]) {\r
+                       this.data = data;\r
+               }\r
+               double[] getData() {\r
+                   return data;        \r
+               }\r
+       }\r
+       \r
          \r
        private class CVodeMemRec {\r
            \r
@@ -451,6 +537,7 @@ public abstract class CVODES {
            --------------------------*/\r
 \r
          //CVRhsFn cv_f;               /* y' = f(t,y(t))                                */\r
+         int cv_f;\r
          //void *cv_user_data;         /* user pointer passed to f                      */\r
 \r
          int cv_lmm;                 /* lmm = ADAMS or BDF                            */\r
@@ -459,7 +546,7 @@ public abstract class CVODES {
          int cv_itol;                /* itol = CV_SS, CV_SV, or CV_WF, or CV_NN       */\r
          double cv_reltol;         /* relative tolerance                            */\r
          double cv_Sabstol;        /* scalar absolute tolerance                     */\r
-         double cv_Vabstol[];        /* vector absolute tolerance                     */\r
+         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
          //void *cv_e_data;            /* user pointer passed to efun                   */\r
@@ -477,7 +564,7 @@ public abstract class CVODES {
          int cv_itolQ;               /* itolQ = CV_SS or CV_SV                        */\r
          double cv_reltolQ;        /* relative tolerance for quadratures            */\r
          double cv_SabstolQ;       /* scalar absolute tolerance for quadratures     */\r
-         double cv_VabstolQ[];       /* vector absolute tolerance for quadratures     */\r
+         NVector cv_VabstolQ = new NVector();       /* vector absolute tolerance for quadratures     */\r
 \r
          /*------------------------\r
            Sensitivity Related Data \r
@@ -506,7 +593,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
-         double cv_VabstolS[];      /* vector absolute tolerances for sensi.        */\r
+         NVector cv_VabstolS = new NVector();      /* vector absolute tolerances for sensi.        */\r
 \r
          /*-----------------------------------\r
            Quadrature Sensitivity Related Data \r
@@ -523,13 +610,13 @@ 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
-         double cv_VabstolQS[];     /* vector absolute tolerances for yQS           */\r
+         NVector cv_VabstolQS = new NVector();     /* vector absolute tolerances for yQS           */\r
 \r
          /*-----------------------\r
            Nordsieck History Array \r
            -----------------------*/\r
 \r
-         double cv_zn[] = new double[L_MAX];      /* Nordsieck array, of size N x (q+1).\r
+         NVector cv_zn[] = new NVector[L_MAX];   /* Nordsieck array, of size N x (q+1).\r
                                         zn[j] is a vector of length N (j=0,...,q)\r
                                         zn[j] = [1/factorial(j)] * h^j * \r
                                         (jth derivative of the interpolating poly.)  */\r
@@ -538,37 +625,37 @@ public abstract class CVODES {
            Vectors of length N \r
            -------------------*/\r
 \r
-         double cv_ewt[];            /* error weight vector                          */\r
-         double cv_y[];              /* y is used as temporary storage by the solver.\r
-                                        The memory is provided by the user to CVode \r
+         NVector cv_ewt = new NVector();            /* error weight vector                          */\r
+         NVector cv_y = new NVector();              // y is used as temporary storage by the solver.\r
+                                     /*   The memory is provided by the user to CVode \r
                                         where the vector is named yout.              */\r
-         double cv_acor[];           /* In the context of the solution of the\r
-                                        nonlinear equation, acor = y_n(m) - y_n(0).\r
+         NVector cv_acor = new NVector();           // In the context of the solution of the\r
+                                      /*  nonlinear equation, acor = y_n(m) - y_n(0).\r
                                         On return, this vector is scaled to give\r
                                         the estimated local error in y.              */\r
-         double cv_tempv[];          /* temporary storage vector                     */\r
-         double cv_ftemp[];          /* temporary storage vector                     */\r
+         NVector cv_tempv = new NVector();          /* temporary storage vector                     */\r
+         NVector cv_ftemp = new NVector();          /* temporary storage vector                     */\r
          \r
          /*--------------------------\r
            Quadrature Related Vectors \r
            --------------------------*/\r
 \r
-         double cv_znQ[] = new double[L_MAX];     /* Nordsieck arrays for quadratures             */\r
-         double cv_ewtQ[];           /* error weight vector for quadratures          */\r
-         double cv_yQ[];             /* Unlike y, yQ is not allocated by the user    */\r
-         double cv_acorQ[];          /* acorQ = yQ_n(m) - yQ_n(0)                    */\r
-         double cv_tempvQ[];         /* temporary storage vector (~ tempv)           */\r
+         NVector cv_znQ[] = new NVector[L_MAX];     /* Nordsieck arrays for quadratures             */\r
+         NVector cv_ewtQ = new NVector();           /* error weight vector for quadratures          */\r
+         NVector cv_yQ = new NVector();             /* Unlike y, yQ is not allocated by the user    */\r
+         NVector cv_acorQ = new NVector();          /* acorQ = yQ_n(m) - yQ_n(0)                    */\r
+         NVector cv_tempvQ = new NVector();         /* temporary storage vector (~ tempv)           */\r
 \r
          /*---------------------------\r
            Sensitivity Related Vectors \r
            ---------------------------*/\r
 \r
-         double cv_znS[][] = new double[L_MAX][];    /* Nordsieck arrays for sensitivities           */\r
-         double cv_ewtS[][];          /* error weight vectors for sensitivities       */\r
-         double cv_yS[][];            /* yS=yS0 (allocated by the user)               */\r
-         double cv_acorS[][];         /* acorS = yS_n(m) - yS_n(0)                    */\r
-         double cv_tempvS[][];        /* temporary storage vector (~ tempv)           */\r
-         double cv_ftempS[][];        /* temporary storage vector (~ ftemp)           */\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
+         NVector cv_tempvS[];        /* temporary storage vector (~ tempv)           */\r
+         NVector cv_ftempS[];        /* temporary storage vector (~ ftemp)           */\r
 \r
          boolean cv_stgr1alloc;  /* Did we allocate ncfS1, ncfnS1, and nniS1?    */\r
 \r
@@ -576,12 +663,12 @@ public abstract class CVODES {
            Quadrature Sensitivity Related Vectors \r
            --------------------------------------*/\r
 \r
-         double cv_znQS[][] = new double[L_MAX][];   /* Nordsieck arrays for quadr. sensitivities    */\r
-         double cv_ewtQS[][];         /* error weight vectors for sensitivities       */\r
-         double cv_yQS[][];           /* Unlike yS, yQS is not allocated by the user  */\r
-         double cv_acorQS[][];        /* acorQS = yQS_n(m) - yQS_n(0)                 */\r
-         double cv_tempvQS[][];       /* temporary storage vector (~ tempv)           */\r
-         double cv_ftempQ[];         /* temporary storage vector (~ ftemp)           */\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
+         NVector cv_tempvQS[];       /* temporary storage vector (~ tempv)           */\r
+         NVector cv_ftempQ = new NVector();         /* temporary storage vector (~ ftemp)           */\r
          \r
          /*-----------------\r
            Tstop information\r
@@ -788,6 +875,7 @@ public abstract class CVODES {
            ----------------*/\r
 \r
          //CVRootFn cv_gfun;        /* Function g for roots sought                     */\r
+         int cv_gfun;\r
          int cv_nrtfn;            /* number of components of g                       */\r
          int cv_iroots[];          /* array for root information                      */\r
          int cv_rootdir[];         /* array specifying direction of zero-crossing     */\r
@@ -857,10 +945,6 @@ public abstract class CVODES {
 \r
          cv_mem = null;\r
          cv_mem = new CVodeMemRec();\r
-         if (cv_mem == null) {\r
-           cvProcessError(null, 0, "CVODES", "CVodeCreate", MSGCV_CVMEM_FAIL);\r
-           return null;\r
-         }\r
 \r
          /* Zero out cv_mem */\r
          //memset(cv_mem, 0, sizeof(struct CVodeMemRec));\r
@@ -879,6 +963,7 @@ public abstract class CVODES {
          /* Set default values for integrator optional inputs */\r
 \r
          //cv_mem.cv_f          = null;\r
+         cv_mem.cv_f = -1;\r
          //cv_mem.cv_user_data  = null;\r
          cv_mem.cv_itol       = CV_NN;\r
          cv_mem.cv_user_efun  = false;\r
@@ -908,6 +993,7 @@ public abstract class CVODES {
          cv_mem.cv_iroots     = null;\r
          cv_mem.cv_rootdir    = null;\r
          //cv_mem.cv_gfun       = null;\r
+         cv_mem.cv_gfun = -1;\r
          cv_mem.cv_nrtfn      = 0;  \r
          cv_mem.cv_gactive    = null;\r
          cv_mem.cv_mxgnull    = 1;\r
@@ -1039,10 +1125,11 @@ public abstract class CVODES {
       * errfp and an error flag is returned. Otherwise, it returns CV_SUCCESS\r
       */\r
 \r
-     int CVodeInit(CVodeMemRec cv_mem, double t0, double y0[])\r
+     int CVodeInit(CVodeMemRec cv_mem, int f, double t0, NVector y0)\r
      {\r
        boolean nvectorOK, allocOK;\r
-       long lrw1, liw1;\r
+       long lrw1[] = new long[1];\r
+       long liw1[] = new long[1];\r
        int i,k;\r
 \r
        /* Check cvode_mem */\r
@@ -1061,11 +1148,11 @@ public abstract class CVODES {
          return(CV_ILL_INPUT);\r
        }\r
 \r
-       //if (f == NULL) {\r
-         //cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVodeInit",\r
-                        //MSGCV_NULL_F);\r
-         //return(CV_ILL_INPUT);\r
-       //}\r
+       if (f < 0) {\r
+         cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVodeInit",\r
+                        MSGCV_NULL_F);\r
+         return(CV_ILL_INPUT);\r
+       }\r
 \r
        /* Test if all required vector operations are implemented */\r
 \r
@@ -1078,97 +1165,106 @@ public abstract class CVODES {
 \r
        /* Set space requirements for one N_Vector */\r
 \r
-       /*if (y0->ops->nvspace != NULL) {\r
-         N_VSpace(y0, &lrw1, &liw1);\r
+       if (y0.getData() != null) {\r
+         N_VSpace_Serial(y0, lrw1, liw1);\r
        } else {\r
-         lrw1 = 0;\r
-         liw1 = 0;\r
+         lrw1[0] = 0;\r
+         liw1[0] = 0;\r
        }\r
-       cv_mem->cv_lrw1 = lrw1;\r
-       cv_mem->cv_liw1 = liw1;*/\r
+       cv_mem.cv_lrw1 = lrw1[0];\r
+       cv_mem.cv_liw1 = liw1[0];\r
 \r
        /* Allocate the vectors (using y0 as a template) */\r
 \r
-      /* allocOK = cvAllocVectors(cv_mem, y0);\r
+      allocOK = cvAllocVectors(cv_mem, y0);\r
        if (!allocOK) {\r
          cvProcessError(cv_mem, CV_MEM_FAIL, "CVODES", "CVodeInit",\r
                         MSGCV_MEM_FAIL);\r
          return(CV_MEM_FAIL);\r
-       }*/\r
+       }\r
 \r
        /* All error checking is complete at this point */\r
 \r
        /* Copy the input parameters into CVODES state */\r
 \r
-       /*cv_mem->cv_f  = f;\r
-       cv_mem->cv_tn = t0;*/\r
+       //cv_mem.cv_f  = f;\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
+       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
+       cv_mem.cv_qu     = 0;\r
+       cv_mem.cv_hu     = ZERO;\r
+       cv_mem.cv_tolsf  = ONE;\r
 \r
        /* Set the linear solver addresses to NULL.\r
           (We check != NULL later, in CVode, if using CV_NEWTON.) */\r
 \r
-       /*cv_mem->cv_linit  = NULL;\r
-       cv_mem->cv_lsetup = NULL;\r
-       cv_mem->cv_lsolve = NULL;\r
-       cv_mem->cv_lfree  = NULL;\r
-       cv_mem->cv_lmem   = NULL;*/\r
+       //cv_mem.cv_linit  = null;\r
+       //cv_mem.cv_lsetup = null;\r
+       //cv_mem.cv_lsolve = null;\r
+       //cv_mem.cv_lfree  = null;\r
+       //cv_mem.cv_lmem   = null;\r
 \r
        /* Set forceSetup to SUNFALSE */\r
 \r
-       //cv_mem->cv_forceSetup = SUNFALSE;\r
+       cv_mem.cv_forceSetup = false;\r
 \r
        /* Initialize zn[0] in the history array */\r
 \r
-       //N_VScale(ONE, y0, cv_mem->cv_zn[0]);\r
+       N_VScale(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;\r
-       cv_mem->cv_netf    = 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
+       cv_mem.cv_nst     = 0;\r
+       cv_mem.cv_nfe     = 0;\r
+       cv_mem.cv_ncfn    = 0;\r
+       cv_mem.cv_netf    = 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
+       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
+       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
        /* NOTE: We do this even if stab lim det was not\r
           turned on yet. This way, the user can turn it\r
           on at any time */\r
 \r
-       /*cv_mem->cv_nor = 0;\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
+           cv_mem.cv_ssdat[i-1][k-1] = ZERO;\r
 \r
        /* Problem has been successfully initialized */\r
 \r
-       /*cv_mem->cv_MallocDone = SUNTRUE;*/\r
+       cv_mem.cv_MallocDone = true;\r
 \r
        return(CV_SUCCESS);\r
      }\r
      \r
+     void N_VSpace_Serial(NVector v, long lrw[], long liw[])\r
+     {\r
+       lrw[0] = v.getData().length;\r
+       liw[0] = 1;\r
+\r
+       return;\r
+     }\r
+\r
+     \r
      /*\r
       * cvCheckNvector\r
       * This routine checks if all required vector operations are present.\r
@@ -1181,21 +1277,373 @@ public abstract class CVODES {
           (tmpl->ops->nvdestroy   == NULL) || // set null\r
           (tmpl->ops->nvlinearsum == NULL) || // (double a, double x[], double b, double y[], double z[])\r
                                               // z[i] = a*x[i] + b*y[i];\r
-          (tmpl->ops->nvconst     == NULL) ||\r
-          (tmpl->ops->nvprod      == NULL) ||\r
-          (tmpl->ops->nvdiv       == NULL) ||\r
-          (tmpl->ops->nvscale     == NULL) ||\r
-          (tmpl->ops->nvabs       == NULL) ||\r
-          (tmpl->ops->nvinv       == NULL) ||\r
-          (tmpl->ops->nvaddconst  == NULL) ||\r
-          (tmpl->ops->nvmaxnorm   == NULL) ||\r
-          (tmpl->ops->nvwrmsnorm  == NULL) ||\r
-          (tmpl->ops->nvmin       == NULL))\r
+          (tmpl->ops->nvconst     == NULL) || // set all to constant\r
+          (tmpl->ops->nvprod      == NULL) || // (double x[], double y[], double z[])\r
+                                              // z[i] = x[i] * y[i]\r
+          (tmpl->ops->nvdiv       == NULL) || // double x[], double y[], double z[])\r
+                                              // z[i] = x[i]/y[i]\r
+          (tmpl->ops->nvscale     == NULL) || // (double c, double x[], double z[])\r
+                                              // z[i] = c * x[i]\r
+          (tmpl->ops->nvabs       == NULL) || // (double x[], double z[])\r
+                                              // z[i] = abs(x[i])\r
+          (tmpl->ops->nvinv       == NULL) || // (double z[], double z[])\r
+                                              // z[i] = 1.0/x[i]\r
+          (tmpl->ops->nvaddconst  == NULL) || // (double x[], double b, double z[])\r
+                                              // z[i] = x[i] + b\r
+          (tmpl->ops->nvmaxnorm   == NULL) || // (double x[])\r
+                                              // max(abs(x[i]))\r
+          (tmpl->ops->nvwrmsnorm  == NULL) || // (double x[], double w[])\r
+                                              // prod = x[i] * w[i]\r
+                                              // sum = sum over all i of prod[i]*prod[i]\r
+                                              // answer = sqrt(sum/N)\r
+          (tmpl->ops->nvmin       == NULL))   // (double x[])\r
+                                              // min(x[i])\r
          return(SUNFALSE);\r
        else\r
          return(SUNTRUE);\r
      }*/\r
+     \r
+     /*\r
+      * cvAllocVectors\r
+      *\r
+      * This routine allocates the CVODES vectors ewt, acor, tempv, ftemp, and\r
+      * zn[0], ..., zn[maxord].\r
+      * If all memory allocations are successful, cvAllocVectors returns SUNTRUE. \r
+      * Otherwise all allocated memory is freed and cvAllocVectors returns SUNFALSE.\r
+      * This routine also sets the optional outputs lrw and liw, which are\r
+      * (respectively) the lengths of the real and integer work spaces\r
+      * allocated here.\r
+      */\r
+\r
+     private boolean cvAllocVectors(CVodeMemRec cv_mem, NVector tmpl)\r
+     {\r
+       int i, j;\r
+\r
+       /* Allocate ewt, acor, tempv, ftemp */\r
+       \r
+       cv_mem.cv_ewt = N_VClone(tmpl);\r
+       if (cv_mem.cv_ewt == null) return(false);\r
+\r
+       cv_mem.cv_acor = N_VClone(tmpl);\r
+       if (cv_mem.cv_acor == null) {\r
+         N_VDestroy(cv_mem.cv_ewt);\r
+         return false;\r
+       }\r
+\r
+       cv_mem.cv_tempv = N_VClone(tmpl);\r
+       if (cv_mem.cv_tempv == null) {\r
+         N_VDestroy(cv_mem.cv_ewt);\r
+         N_VDestroy(cv_mem.cv_acor);\r
+         return false;\r
+       }\r
+\r
+       cv_mem.cv_ftemp = N_VClone(tmpl);\r
+       if (cv_mem.cv_ftemp == null) {\r
+         N_VDestroy(cv_mem.cv_tempv);\r
+         N_VDestroy(cv_mem.cv_ewt);\r
+         N_VDestroy(cv_mem.cv_acor);\r
+         return false;\r
+       }\r
+\r
+       /* Allocate zn[0] ... zn[qmax] */\r
+\r
+       for (j=0; j <= cv_mem.cv_qmax; j++) {\r
+         cv_mem.cv_zn[j] = N_VClone(tmpl);\r
+         if (cv_mem.cv_zn[j] == null) {\r
+           N_VDestroy(cv_mem.cv_ewt);\r
+           N_VDestroy(cv_mem.cv_acor);\r
+           N_VDestroy(cv_mem.cv_tempv);\r
+           N_VDestroy(cv_mem.cv_ftemp);\r
+           for (i=0; i < j; i++) N_VDestroy(cv_mem.cv_zn[i]);\r
+           return false;\r
+         }\r
+       }\r
+\r
+       /* Update solver workspace lengths  */\r
+       cv_mem.cv_lrw += (cv_mem.cv_qmax + 5)*cv_mem.cv_lrw1;\r
+       cv_mem.cv_liw += (cv_mem.cv_qmax + 5)*cv_mem.cv_liw1;\r
+\r
+       /* Store the value of qmax used here */\r
+       cv_mem.cv_qmax_alloc = cv_mem.cv_qmax;\r
+\r
+       return true;\r
+     }\r
+     \r
+     private NVector N_VClone(NVector y) {\r
+         NVector x = new NVector();\r
+         x.data = y.data.clone();\r
+         x.own_data = y.own_data;\r
+         return x;\r
+     }\r
+\r
+     private void N_VDestroy(NVector x) {\r
+        x.data = null;\r
+        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
+       if (cv_mem==null) {\r
+         cvProcessError(null, CV_MEM_NULL, "CVODES",\r
+                        "CVodeSVtolerances", MSGCV_NO_MEM);\r
+         return(CV_MEM_NULL);\r
+       }\r
 \r
+       if (cv_mem.cv_MallocDone == false) {\r
+         cvProcessError(cv_mem, CV_NO_MALLOC, "CVODES",\r
+                        "CVodeSVtolerances", MSGCV_NO_MALLOC);\r
+         return(CV_NO_MALLOC);\r
+       }\r
+\r
+       /* Check inputs */\r
+\r
+       if (reltol < ZERO) {\r
+         cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES",\r
+                        "CVodeSVtolerances", MSGCV_BAD_RELTOL);\r
+         return(CV_ILL_INPUT);\r
+       }\r
+\r
+       if (N_VMin(abstol) < ZERO) {\r
+         cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES",\r
+                        "CVodeSVtolerances", MSGCV_BAD_ABSTOL);\r
+         return(CV_ILL_INPUT);\r
+       }\r
+\r
+       /* Copy tolerances into memory */\r
+       \r
+       if ( !(cv_mem.cv_VabstolMallocDone) ) {\r
+         cv_mem.cv_Vabstol = N_VClone(cv_mem.cv_ewt);\r
+         cv_mem.cv_lrw += cv_mem.cv_lrw1;\r
+         cv_mem.cv_liw += cv_mem.cv_liw1;\r
+         cv_mem.cv_VabstolMallocDone = true;\r
+       }\r
+\r
+       cv_mem.cv_reltol = reltol;\r
+       N_VScale(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
+\r
+       return(CV_SUCCESS);\r
+     }\r
+\r
+     private double N_VMin(NVector x) {\r
+        int i;\r
+        if ((x == null) || (x.data == null) || (x.data.length == 0)) {\r
+                return Double.NaN;\r
+        }\r
+        double data[] = x.data;\r
+        int n = data.length;\r
+        double min = data[0];\r
+        for (i = 1; i < n; i++) {\r
+                if (data[i] < min) {\r
+                        min = data[i];\r
+                }\r
+        }\r
+        return min;\r
+     }\r
+     \r
+     /*\r
+      * CVodeRootInit\r
+      *\r
+      * CVodeRootInit initializes a rootfinding problem to be solved\r
+      * during the integration of the ODE system.  It loads the root\r
+      * function pointer and the number of root functions, and allocates\r
+      * workspace memory.  The return value is CV_SUCCESS = 0 if no errors\r
+      * occurred, or a negative value otherwise.\r
+      */\r
+\r
+     int CVodeRootInit(CVodeMemRec cv_mem, int nrtfn, int g)\r
+     {\r
+       int i, nrt;\r
+\r
+       /* Check cvode_mem */\r
+       if (cv_mem==null) {\r
+         cvProcessError(null, CV_MEM_NULL, "CVODES", "CVodeRootInit",\r
+                        MSGCV_NO_MEM);\r
+         return(CV_MEM_NULL);\r
+       }\r
+\r
+       nrt = (nrtfn < 0) ? 0 : nrtfn;\r
+\r
+       /* If rerunning CVodeRootInit() with a different number of root\r
+          functions (changing number of gfun components), then free\r
+          currently held memory resources */\r
+       if ((nrt != cv_mem.cv_nrtfn) && (cv_mem.cv_nrtfn > 0)) {\r
+         cv_mem.cv_glo = null;\r
+         cv_mem.cv_ghi = null;\r
+         cv_mem.cv_grout = null;\r
+         cv_mem.cv_iroots = null;\r
+         cv_mem.cv_rootdir = null;\r
+         cv_mem.cv_gactive = null;\r
+\r
+         cv_mem.cv_lrw -= 3 * (cv_mem.cv_nrtfn);\r
+         cv_mem.cv_liw -= 3 * (cv_mem.cv_nrtfn);\r
+       }\r
+\r
+       /* If CVodeRootInit() was called with nrtfn == 0, then set cv_nrtfn to\r
+          zero and cv_gfun to NULL before returning */\r
+       if (nrt == 0) {\r
+         cv_mem.cv_nrtfn = nrt;\r
+         cv_mem.cv_gfun = -1;\r
+         return(CV_SUCCESS);\r
+       }\r
+\r
+       /* If rerunning CVodeRootInit() with the same number of root functions\r
+          (not changing number of gfun components), then check if the root\r
+          function argument has changed */\r
+       /* If g != NULL then return as currently reserved memory resources\r
+          will suffice */\r
+       if (nrt == cv_mem.cv_nrtfn) {\r
+         if (g != cv_mem.cv_gfun) {\r
+           if (g < 0) {\r
+            cv_mem.cv_glo = null;\r
+            cv_mem.cv_ghi = null;\r
+            cv_mem.cv_grout = null;\r
+            cv_mem.cv_iroots = null;\r
+             cv_mem.cv_rootdir = null;\r
+             cv_mem.cv_gactive = null;\r
+\r
+             cv_mem.cv_lrw -= 3*nrt;\r
+             cv_mem.cv_liw -= 3*nrt;\r
+\r
+             cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVodeRootInit",\r
+                            MSGCV_NULL_G);\r
+             return(CV_ILL_INPUT);\r
+           }\r
+           else {\r
+            cv_mem.cv_gfun = g;\r
+             return(CV_SUCCESS);\r
+           }\r
+         }\r
+         else return(CV_SUCCESS);\r
+       }\r
+\r
+       /* Set variable values in CVode memory block */\r
+       cv_mem.cv_nrtfn = nrt;\r
+       if (g < 0) {\r
+         cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVodeRootInit",\r
+                        MSGCV_NULL_G);\r
+         return(CV_ILL_INPUT);\r
+       }\r
+       else cv_mem.cv_gfun = g;\r
+\r
+       /* Allocate necessary memory and return */\r
+       cv_mem.cv_glo = null;\r
+       cv_mem.cv_glo = new double[nrt];\r
+       if (cv_mem.cv_glo == null) {\r
+         cvProcessError(cv_mem, CV_MEM_FAIL, "CVODES", "CVodeRootInit",\r
+                        MSGCV_MEM_FAIL);\r
+         return(CV_MEM_FAIL);\r
+       }\r
+         \r
+       cv_mem.cv_ghi = null;\r
+       cv_mem.cv_ghi = new double[nrt];\r
+       if (cv_mem.cv_ghi == null) {\r
+         cv_mem.cv_glo = null;\r
+         cvProcessError(cv_mem, CV_MEM_FAIL, "CVODES", "CVodeRootInit",\r
+                        MSGCV_MEM_FAIL);\r
+         return(CV_MEM_FAIL);\r
+       }\r
+         \r
+       cv_mem.cv_grout = null;\r
+       cv_mem.cv_grout = new double[nrt];\r
+       if (cv_mem.cv_grout == null) {\r
+         cv_mem.cv_glo = null;\r
+         cv_mem.cv_ghi = null;\r
+         cvProcessError(cv_mem, CV_MEM_FAIL, "CVODES", "CVodeRootInit",\r
+                        MSGCV_MEM_FAIL);\r
+         return(CV_MEM_FAIL);\r
+       }\r
+\r
+       cv_mem.cv_iroots = null;\r
+       cv_mem.cv_iroots = new int[nrt];\r
+       if (cv_mem.cv_iroots == null) {\r
+         cv_mem.cv_glo = null;\r
+         cv_mem.cv_ghi = null;\r
+         cv_mem.cv_grout = null;\r
+         cvProcessError(cv_mem, CV_MEM_FAIL, "CVODES", "CVodeRootInit",\r
+                        MSGCV_MEM_FAIL);\r
+         return(CV_MEM_FAIL);\r
+       }\r
+\r
+       cv_mem.cv_rootdir = null;\r
+       cv_mem.cv_rootdir = new int[nrt];\r
+       if (cv_mem.cv_rootdir == null) {\r
+         cv_mem.cv_glo = null; \r
+         cv_mem.cv_ghi = null;\r
+         cv_mem.cv_grout = null;\r
+         cv_mem.cv_iroots = null;\r
+         cvProcessError(cv_mem, CV_MEM_FAIL, "CVODES", "CVodeRootInit",\r
+                        MSGCV_MEM_FAIL);\r
+         return(CV_MEM_FAIL);\r
+       }\r
+\r
+\r
+       cv_mem.cv_gactive = null;\r
+       cv_mem.cv_gactive = new boolean[nrt];\r
+       if (cv_mem.cv_gactive == null) {\r
+         cv_mem.cv_glo = null; \r
+         cv_mem.cv_ghi = null;\r
+         cv_mem.cv_grout = null;\r
+         cv_mem.cv_iroots = null;\r
+         cv_mem.cv_rootdir = null;\r
+         cvProcessError(cv_mem, CV_MEM_FAIL, "CVODES", "CVodeRootInit",\r
+                        MSGCV_MEM_FAIL);\r
+         return(CV_MEM_FAIL);\r
+       }\r
+\r
+\r
+       /* Set default values for rootdir (both directions) */\r
+       for(i=0; i<nrt; i++) cv_mem.cv_rootdir[i] = 0;\r
+\r
+       /* Set default values for gactive (all active) */\r
+       for(i=0; i<nrt; i++) cv_mem.cv_gactive[i] = true;\r
+\r
+       cv_mem.cv_lrw += 3*nrt;\r
+       cv_mem.cv_liw += 3*nrt;\r
+\r
+       return(CV_SUCCESS);\r
+     }\r
+     \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
+     }\r
+     \r
+     private SUNLinearSolver SUNDenseLinearSolver(NVector y, double A[][]) {\r
+        int matrixRows = A.length;\r
+        int matrixColumns = A[0].length;\r
+        if (matrixRows != matrixColumns) {\r
+                MipavUtil.displayError("Did not have the matrix rows = matrix columns required for a LinearSolver");\r
+                return null;\r
+        }\r
+        int vecLength = y.getData().length;\r
+        if (matrixRows != vecLength) {\r
+                MipavUtil.displayError("Did not have the matrix rows equal to vector length required for a LinearSolver");\r
+                return null;\r
+        }\r
+        SUNLinearSolver S = new SUNLinearSolver();\r
+        S.N = matrixRows;\r
+        S.last_flag = 0;\r
+        S.pivots = new long[matrixRows];\r
+        return S;\r
+     }\r
 \r
 \r
 }
\ No newline at end of file