Porting continues.
authorilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Thu, 1 Mar 2018 23:09:33 +0000 (23:09 +0000)
committerilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Thu, 1 Mar 2018 23:09:33 +0000 (23:09 +0000)
git-svn-id: https://citdcbmipav.cit.nih.gov/repos-pub/mipav/trunk@15395 ba61647d-9d00-f842-95cd-605cb4296b96

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

index b7b5af2..97031e8 100644 (file)
@@ -222,6 +222,16 @@ public abstract class CVODES {
        final double HLB_FACTOR = 100.0;\r
        final double FUZZ_FACTOR = 100.0;\r
        final int MAX_ITERS = 4;\r
+       // CRDOWN constant used in the estimation of the convergence rate (crate)\r
+       //        of the iterates for the nonlinear equation\r
+       final double CRDOWN = 0.3;\r
+       // DGMAX  iter == CV_NEWTON, |gamma/gammap-1| > DGMAX => call lsetup\r
+    final double DGMAX = 0.3;\r
+    // RDIV declare divergence if ratio del/delp > RDIV\r
+    final double RDIV = TWO;\r
+    // MSBP  max no. of steps between lsetup calls\r
+    final int MSBP = 20;\r
+\r
        \r
        final int RTFOUND = +1;\r
        final int CLOSERT = +3;\r
@@ -418,8 +428,53 @@ public abstract class CVODES {
        final int SRHSFUNC_RECVR = +12;\r
        final int QSRHSFUNC_RECVR = +13;\r
 \r
+       /*\r
+        * -----------------------------------------------------------------\r
+        * Communication between CVODE and a CVODE Linear Solver\r
+        * -----------------------------------------------------------------\r
+        * convfail (input to cv_lsetup)\r
+        *\r
+        * CV_NO_FAILURES : Either this is the first cv_setup call for this\r
+        *                  step, or the local error test failed on the\r
+        *                  previous attempt at this step (but the Newton\r
+        *                  iteration converged).\r
+        *\r
+        * CV_FAIL_BAD_J  : This value is passed to cv_lsetup if\r
+        *\r
+        *                  (a) The previous Newton corrector iteration\r
+        *                      did not converge and the linear solver's\r
+        *                      setup routine indicated that its Jacobian-\r
+        *                      related data is not current\r
+        *                                   or\r
+        *                  (b) During the previous Newton corrector\r
+        *                      iteration, the linear solver's solve routine\r
+        *                      failed in a recoverable manner and the\r
+        *                      linear solver's setup routine indicated that\r
+        *                      its Jacobian-related data is not current.\r
+        *\r
+        * CV_FAIL_OTHER  : During the current internal step try, the\r
+        *                  previous Newton iteration failed to converge\r
+        *                  even though the linear solver was using current\r
+        *                  Jacobian-related data.\r
+        * -----------------------------------------------------------------\r
+        */\r
+       /* Constants for convfail (input to cv_lsetup) */\r
+       final int CV_NO_FAILURES = 0;\r
+    final int CV_FAIL_BAD_J = 1;\r
+       final int CV_FAIL_OTHER = 2;\r
        \r
        /*-----------------------------------------------------------------\r
+         CVDLS solver constants\r
+         -----------------------------------------------------------------\r
+         CVD_MSBJ   maximum number of steps between Jacobian evaluations\r
+         CVD_DGMAX  maximum change in gamma between Jacobian evaluations\r
+         -----------------------------------------------------------------*/\r
+\r
+       final int CVD_MSBJ = 50;\r
+       final double CVD_DGMAX = 0.2;\r
+\r
+\r
+       /*-----------------------------------------------------------------\r
         * IV. SUNLinearSolver error codes\r
         * ---------------------------------------------------------------\r
         */\r
@@ -686,7 +741,7 @@ public abstract class CVODES {
         * @param tmp2\r
         * @return\r
         */\r
-       private int Jac(double t, NVector yv, NVector fy, double J[][], NVector tmp1, NVector tmp2) {\r
+       private int Jac(double t, NVector yv, NVector fy, double J[][], CVodeMemRec data, NVector tmp1, NVector tmp2, NVector tmp3) {\r
                double y[] = yv.getData();\r
            J[0][0] = -0.04;\r
            J[0][1] = 1.0E4 * y[2];\r
@@ -1045,7 +1100,7 @@ public abstract class CVODES {
          double cv_h0u;             /* actual initial stepsize                     */\r
          double cv_hu;              /* last successful h value used                */\r
          double cv_saved_tq5;       /* saved value of tq[5]                        */\r
-         boolean cv_jcur;         /* is Jacobian info for linear solver current? */\r
+         boolean cv_jcur[] = new boolean[1];         /* is Jacobian info for linear solver current? */\r
          int cv_convfail;             /* flag storing previous solver failure mode   */\r
          double cv_tolsf;           /* tolerance scale factor                      */\r
          int cv_qmax_alloc;           /* qmax used when allocating mem               */\r
@@ -5715,10 +5770,10 @@ else                return(snrm);
      if (do_sensi_sim)\r
        delS = cvSensUpdateNorm(cv_mem, del, cv_mem.cv_acorS, cv_mem.cv_ewtS);\r
 \r
-     /*N_VScale(ONE, cv_mem->cv_tempv, cv_mem->cv_acor);\r
+     N_VScale_Serial(ONE, cv_mem.cv_tempv, cv_mem.cv_acor);\r
      if (do_sensi_sim) \r
-       for (is=0; is<cv_mem->cv_Ns; is++)\r
-         N_VScale(ONE, cv_mem->cv_tempvS[is], cv_mem->cv_acorS[is]);*/\r
+       for (is=0; is<cv_mem.cv_Ns; is++)\r
+         N_VScale_Serial(ONE, cv_mem.cv_tempvS[is], cv_mem.cv_acorS[is]);\r
      \r
      /* Test for convergence.  If m > 0, an estimate of the convergence\r
         rate constant is stored in crate, and used in the test. \r
@@ -5730,50 +5785,50 @@ else                return(snrm);
         del and delS)\r
      */\r
      \r
-     /*Del = (do_sensi_sim) ? delS : del;\r
-     if (m > 0) cv_mem->cv_crate = SUNMAX(CRDOWN * cv_mem->cv_crate, Del / Delp);\r
-     dcon = Del * SUNMIN(ONE, cv_mem->cv_crate) / cv_mem->cv_tq[4];\r
+     Del = (do_sensi_sim) ? delS : del;\r
+     if (m > 0) cv_mem.cv_crate = Math.max(CRDOWN * cv_mem.cv_crate, Del / Delp);\r
+     dcon = Del * Math.min(ONE, cv_mem.cv_crate) / cv_mem.cv_tq[4];\r
 \r
      if (dcon <= ONE) {\r
        if (m == 0)\r
-         if (do_sensi_sim && cv_mem->cv_errconS)\r
-           cv_mem->cv_acnrm = delS;\r
+         if (do_sensi_sim && cv_mem.cv_errconS)\r
+           cv_mem.cv_acnrm = delS;\r
          else\r
-           cv_mem->cv_acnrm = del;\r
+           cv_mem.cv_acnrm = del;\r
        else {\r
-         cv_mem->cv_acnrm = N_VWrmsNorm(cv_mem->cv_acor, cv_mem->cv_ewt);\r
-         if (do_sensi_sim && cv_mem->cv_errconS)\r
-           cv_mem->cv_acnrm = cvSensUpdateNorm(cv_mem, cv_mem->cv_acnrm,\r
-                                               cv_mem->cv_acorS, cv_mem->cv_ewtS);\r
+         cv_mem.cv_acnrm = N_VWrmsNorm_Serial(cv_mem.cv_acor, cv_mem.cv_ewt);\r
+         if (do_sensi_sim && cv_mem.cv_errconS)\r
+           cv_mem.cv_acnrm = cvSensUpdateNorm(cv_mem, cv_mem.cv_acnrm,\r
+                                               cv_mem.cv_acorS, cv_mem.cv_ewtS);\r
        }\r
-       return(CV_SUCCESS);*/  /* Convergence achieved */\r
-     /*}*/\r
+       return(CV_SUCCESS);  /* Convergence achieved */\r
+     }\r
 \r
      /* Stop at maxcor iterations or if iter. seems to be diverging */\r
 \r
-     /*m++;\r
-     if ((m==cv_mem->cv_maxcor) || ((m >= 2) && (Del > RDIV * Delp)))\r
-       return(CONV_FAIL);*/\r
+     m++;\r
+     if ((m==cv_mem.cv_maxcor) || ((m >= 2) && (Del > RDIV * Delp)))\r
+       return(CONV_FAIL);\r
 \r
      /* Save norm of correction, evaluate f, and loop again */\r
 \r
-    /* Delp = Del;\r
+     Delp = Del;\r
 \r
-     retval = cv_mem->cv_f(cv_mem->cv_tn, cv_mem->cv_y,\r
-                           cv_mem->cv_tempv, cv_mem->cv_user_data);\r
-     cv_mem->cv_nfe++;\r
+     retval = f(cv_mem.cv_tn, cv_mem.cv_y,\r
+                           cv_mem.cv_tempv, cv_mem.cv_user_data);\r
+     cv_mem.cv_nfe++;\r
      if (retval < 0) return(CV_RHSFUNC_FAIL);\r
      if (retval > 0) return(RHSFUNC_RECVR);\r
 \r
      if (do_sensi_sim) {\r
-       wrk1 = cv_mem->cv_ftemp;\r
-       wrk2 = cv_mem->cv_ftempS[0];\r
-       retval = cvSensRhsWrapper(cv_mem, cv_mem->cv_tn, cv_mem->cv_y,\r
-                                 cv_mem->cv_tempv, cv_mem->cv_yS,\r
-                                 cv_mem->cv_tempvS, wrk1, wrk2);\r
+       wrk1 = cv_mem.cv_ftemp;\r
+       wrk2 = cv_mem.cv_ftempS[0];\r
+       retval = cvSensRhsWrapper(cv_mem, cv_mem.cv_tn, cv_mem.cv_y,\r
+                                 cv_mem.cv_tempv, cv_mem.cv_yS,\r
+                                 cv_mem.cv_tempvS, wrk1, wrk2);\r
        if (retval < 0) return(CV_SRHSFUNC_FAIL);\r
        if (retval > 0) return(SRHSFUNC_RECVR);\r
-     }  */    \r
+     }  \r
 \r
    } /* end loop */\r
 \r
@@ -5813,103 +5868,553 @@ else                return(snrm);
 \r
  private int cvNlsNewton(CVodeMemRec cv_mem, int nflag)\r
  {\r
-  /* N_Vector vtemp1, vtemp2, vtemp3, wrk1, wrk2;\r
+   NVector vtemp1, vtemp2, vtemp3, wrk1, wrk2;\r
    int convfail, ier;\r
-   booleantype callSetup, do_sensi_sim;\r
-   int retval, is;*/\r
+   boolean callSetup, do_sensi_sim;\r
+   int retval, is;\r
    \r
    /* Are we computing sensitivities with the CV_SIMULTANEOUS approach? */\r
-  /* do_sensi_sim = (cv_mem->cv_sensi && (cv_mem->cv_ism==CV_SIMULTANEOUS));\r
+   do_sensi_sim = (cv_mem.cv_sensi && (cv_mem.cv_ism==CV_SIMULTANEOUS));\r
 \r
-   vtemp1 = cv_mem->cv_acor; */ /* rename acor as vtemp1 for readability  */\r
-  /* vtemp2 = cv_mem->cv_y;  */   /* rename y as vtemp2 for readability     */\r
-  /* vtemp3 = cv_mem->cv_tempv; *//* rename tempv as vtemp3 for readability */\r
+   vtemp1 = cv_mem.cv_acor;  /* rename acor as vtemp1 for readability  */\r
+   vtemp2 = cv_mem.cv_y;     /* rename y as vtemp2 for readability     */\r
+   vtemp3 = cv_mem.cv_tempv; /* rename tempv as vtemp3 for readability */\r
    \r
    /* Set flag convfail, input to lsetup for its evaluation decision */\r
-  /* convfail = ((nflag == FIRST_CALL) || (nflag == PREV_ERR_FAIL)) ?\r
-     CV_NO_FAILURES : CV_FAIL_OTHER;*/\r
+   convfail = ((nflag == FIRST_CALL) || (nflag == PREV_ERR_FAIL)) ?\r
+     CV_NO_FAILURES : CV_FAIL_OTHER;\r
 \r
    /* Decide whether or not to call setup routine (if one exists) */\r
-   /*if (cv_mem->cv_lsetup) {      \r
+   if (cv_mem.cv_lsetup_select > 0) {      \r
      callSetup = (nflag == PREV_CONV_FAIL) || (nflag == PREV_ERR_FAIL) ||\r
-       (cv_mem->cv_nst == 0) ||\r
-       (cv_mem->cv_nst >= cv_mem->cv_nstlp + MSBP) ||\r
-       (SUNRabs(cv_mem->cv_gamrat-ONE) > DGMAX);*/\r
+       (cv_mem.cv_nst == 0) ||\r
+       (cv_mem.cv_nst >= cv_mem.cv_nstlp + MSBP) ||\r
+       (Math.abs(cv_mem.cv_gamrat-ONE) > DGMAX);\r
 \r
      /* Decide whether to force a call to setup */\r
-    /* if (cv_mem->cv_forceSetup) {\r
-       callSetup = SUNTRUE;\r
+    if (cv_mem.cv_forceSetup) {\r
+       callSetup = true;\r
        convfail = CV_FAIL_OTHER;\r
      }\r
 \r
    } else {  \r
-     cv_mem->cv_crate = ONE;\r
-     cv_mem->cv_crateS = ONE; */ /* if NO lsetup all conv. rates are set to ONE */\r
-   /*  callSetup = SUNFALSE;\r
-   }*/\r
+     cv_mem.cv_crate = ONE;\r
+     cv_mem.cv_crateS = ONE; /* if NO lsetup all conv. rates are set to ONE */\r
+     callSetup = false;\r
+   }\r
    \r
    /* Looping point for the solution of the nonlinear system.\r
       Evaluate f at the predicted y, call lsetup if indicated, and\r
       call cvNewtonIteration for the Newton iteration itself.      */\r
 \r
-   /*for(;;) {\r
+   for(;;) {\r
 \r
-     retval = cv_mem->cv_f(cv_mem->cv_tn, cv_mem->cv_zn[0],\r
-                           cv_mem->cv_ftemp, 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_ftemp, cv_mem.cv_user_data);\r
+     cv_mem.cv_nfe++; \r
      if (retval < 0) return(CV_RHSFUNC_FAIL);\r
      if (retval > 0) return(RHSFUNC_RECVR);\r
 \r
      if (do_sensi_sim) {\r
-       wrk1 = cv_mem->cv_tempv;\r
-       wrk2 = cv_mem->cv_tempvS[0];\r
-       retval = cvSensRhsWrapper(cv_mem, cv_mem->cv_tn, cv_mem->cv_zn[0],\r
-                                 cv_mem->cv_ftemp, cv_mem->cv_znS[0],\r
-                                 cv_mem->cv_ftempS, wrk1, wrk2);\r
+       wrk1 = cv_mem.cv_tempv;\r
+       wrk2 = cv_mem.cv_tempvS[0];\r
+       retval = cvSensRhsWrapper(cv_mem, cv_mem.cv_tn, cv_mem.cv_zn[0],\r
+                                 cv_mem.cv_ftemp, cv_mem.cv_znS[0],\r
+                                 cv_mem.cv_ftempS, wrk1, wrk2);\r
        if (retval < 0) return(CV_SRHSFUNC_FAIL);\r
        if (retval > 0) return(SRHSFUNC_RECVR);\r
      }\r
 \r
      if (callSetup) {\r
-       ier = cv_mem->cv_lsetup(cv_mem, convfail, cv_mem->cv_zn[0],\r
-                               cv_mem->cv_ftemp, &(cv_mem->cv_jcur), \r
-                               vtemp1, vtemp2, vtemp3);\r
-       cv_mem->cv_nsetups++;\r
-       callSetup = SUNFALSE;\r
-       cv_mem->cv_forceSetup = SUNFALSE;\r
-       cv_mem->cv_gamrat = ONE; \r
-       cv_mem->cv_gammap = cv_mem->cv_gamma;\r
-       cv_mem->cv_crate = ONE;\r
-       cv_mem->cv_crateS = ONE;*/ /* after lsetup all conv. rates are reset to ONE */\r
-       /*cv_mem->cv_nstlp = cv_mem->cv_nst;*/\r
+       ier = cv_lsetup(cv_mem, convfail, cv_mem.cv_zn[0],\r
+                               cv_mem.cv_ftemp, cv_mem.cv_jcur,\r
+                               vtemp1, vtemp2, vtemp3, cv_mem.cv_lsetup_select);\r
+       cv_mem.cv_nsetups++;\r
+       callSetup = false;\r
+       cv_mem.cv_forceSetup = false;\r
+       cv_mem.cv_gamrat = ONE; \r
+       cv_mem.cv_gammap = cv_mem.cv_gamma;\r
+       cv_mem.cv_crate = ONE;\r
+       cv_mem.cv_crateS = ONE; /* after lsetup all conv. rates are reset to ONE */\r
+       cv_mem.cv_nstlp = cv_mem.cv_nst;\r
        /* Return if lsetup failed */\r
-       /*if (ier < 0) return(CV_LSETUP_FAIL);\r
+       if (ier < 0) return(CV_LSETUP_FAIL);\r
        if (ier > 0) return(CONV_FAIL);\r
-     }*/\r
+     }\r
 \r
      /* Set acor to zero and load prediction into y vector */\r
-    /* N_VConst(ZERO, cv_mem->cv_acor);\r
-     N_VScale(ONE, cv_mem->cv_zn[0], cv_mem->cv_y);\r
+     N_VConst_Serial(ZERO, cv_mem.cv_acor);\r
+     N_VScale_Serial(ONE, cv_mem.cv_zn[0], cv_mem.cv_y);\r
 \r
      if (do_sensi_sim)\r
-       for (is=0; is<cv_mem->cv_Ns; is++) {\r
-         N_VConst(ZERO, cv_mem->cv_acorS[is]);\r
-         N_VScale(ONE, cv_mem->cv_znS[0][is], cv_mem->cv_yS[is]);\r
-       }*/\r
+       for (is=0; is<cv_mem.cv_Ns; is++) {\r
+         N_VConst_Serial(ZERO, cv_mem.cv_acorS[is]);\r
+         N_VScale_Serial(ONE, cv_mem.cv_znS[0][is], cv_mem.cv_yS[is]);\r
+       }\r
 \r
      /* Do the Newton iteration */\r
-    /* ier = cvNewtonIteration(cv_mem);*/\r
+       ier = cvNewtonIteration(cv_mem);\r
 \r
      /* If there is a convergence failure and the Jacobian-related \r
         data appears not to be current, loop again with a call to lsetup\r
         in which convfail=CV_FAIL_BAD_J.  Otherwise return.                 */\r
-    /* if (ier != TRY_AGAIN) return(ier);\r
+     if (ier != TRY_AGAIN) return(ier);\r
      \r
-     callSetup = SUNTRUE;\r
+     callSetup = true;\r
      convfail = CV_FAIL_BAD_J;\r
-   }*/\r
-        return -100;\r
+   }\r
  }\r
+     \r
+     int cv_lsetup(CVodeMemRec cv_mem, int convfail, NVector y, \r
+             NVector fy, boolean jcurPtr[], \r
+             NVector tmp1, NVector tmp2, NVector tmp3, int select) {\r
+        int flag = 0;\r
+        switch(select) {\r
+        case cvDlsSetup_select:\r
+                flag = cvDlsSetup(cv_mem, convfail, y, fy, jcurPtr,\r
+                                tmp1, tmp2, tmp3);\r
+        }\r
+        return flag;\r
+     }\r
+     \r
+     /*-----------------------------------------------------------------\r
+     cvDlsSetup\r
+     -----------------------------------------------------------------\r
+     This routine determines whether to update a Jacobian matrix (or\r
+     use a stored version), based on heuristics regarding previous \r
+     convergence issues, the number of time steps since it was last\r
+     updated, etc.; it then creates the system matrix from this, the\r
+     'gamma' factor and the identity matrix, \r
+       A = I-gamma*J.\r
+     This routine then calls the LS 'setup' routine with A.\r
+     -----------------------------------------------------------------*/\r
+   int cvDlsSetup(CVodeMemRec cv_mem, int convfail, NVector y, \r
+                  NVector fy, boolean jcurPtr[], \r
+                  NVector tmp1, NVector tmp2, NVector tmp3)\r
+   {\r
+     boolean jbad, jok;\r
+     double dgamma;\r
+     CVDlsMemRec cvdls_mem;\r
+     int retval;\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
+                       "cvDlsSetup", 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
+                       "cvDlsSetup", MSGD_LMEM_NULL);\r
+       return(CVDLS_LMEM_NULL);\r
+     }\r
+     cvdls_mem = cv_mem.cv_lmem;\r
+\r
+     /* Use nst, gamma/gammap, and convfail to set J eval. flag jok */\r
+     dgamma = Math.abs((cv_mem.cv_gamma/cv_mem.cv_gammap) - ONE);\r
+     jbad = (cv_mem.cv_nst == 0) || \r
+       (cv_mem.cv_nst > cvdls_mem.nstlj + CVD_MSBJ) ||\r
+       ((convfail == CV_FAIL_BAD_J) && (dgamma < CVD_DGMAX)) ||\r
+       (convfail == CV_FAIL_OTHER);\r
+     jok = !jbad;\r
+    \r
+     /* If jok = SUNTRUE, use saved copy of J */\r
+     if (jok) {\r
+       jcurPtr[0] = false;\r
+       // Copy savedJ to A\r
+       retval = SUNMatCopy(cvdls_mem.savedJ, cvdls_mem.A);\r
+       if (retval != 0) {\r
+         cvProcessError(cv_mem, CVDLS_SUNMAT_FAIL, "CVSDLS", \r
+                         "cvDlsSetup",  MSGD_MATCOPY_FAILED);\r
+         cvdls_mem.last_flag = CVDLS_SUNMAT_FAIL;\r
+         return(-1);\r
+       }\r
+\r
+     /* If jok = SUNFALSE, call jac routine for new J value */\r
+     } else {\r
+       cvdls_mem.nje++;\r
+       cvdls_mem.nstlj = cv_mem.cv_nst;\r
+       jcurPtr[0] = true;\r
+       retval = SUNMatZero(cvdls_mem.A);\r
+       if (retval != 0) {\r
+         cvProcessError(cv_mem, CVDLS_SUNMAT_FAIL, "CVSDLS", \r
+                         "cvDlsSetup",  MSGD_MATZERO_FAILED);\r
+         cvdls_mem.last_flag = CVDLS_SUNMAT_FAIL;\r
+         return(-1);\r
+       }\r
+\r
+       retval = Jac(cv_mem.cv_tn, y, fy, cvdls_mem.A, \r
+                               cvdls_mem.J_data, tmp1, tmp2, tmp3);\r
+       if (retval < 0) {\r
+         cvProcessError(cv_mem, CVDLS_JACFUNC_UNRECVR, "CVSDLS", \r
+                         "cvDlsSetup",  MSGD_JACFUNC_FAILED);\r
+         cvdls_mem.last_flag = CVDLS_JACFUNC_UNRECVR;\r
+         return(-1);\r
+       }\r
+       if (retval > 0) {\r
+         cvdls_mem.last_flag = CVDLS_JACFUNC_RECVR;\r
+         return(1);\r
+       }\r
+\r
+       retval = SUNMatCopy(cvdls_mem.A, cvdls_mem.savedJ);\r
+       if (retval != 0) {\r
+         cvProcessError(cv_mem, CVDLS_SUNMAT_FAIL, "CVSDLS", \r
+                         "cvDlsSetup",  MSGD_MATCOPY_FAILED);\r
+         cvdls_mem.last_flag = CVDLS_SUNMAT_FAIL;\r
+         return(-1);\r
+       }\r
+\r
+     }\r
+     \r
+     /* Scale and add I to get A = I - gamma*J */\r
+     retval = SUNMatScaleAddI(-cv_mem.cv_gamma, cvdls_mem.A);\r
+     if (retval != 0) {\r
+       cvProcessError(cv_mem, CVDLS_SUNMAT_FAIL, "CVSDLS", \r
+                      "cvDlsSetup",  MSGD_MATSCALEADDI_FAILED);\r
+       cvdls_mem.last_flag = CVDLS_SUNMAT_FAIL;\r
+       return(-1);\r
+     }\r
+\r
+     /* Call generic linear solver 'setup' with this system matrix, and\r
+        return success/failure flag */\r
+     cvdls_mem.last_flag = SUNLinSolSetup_Dense(cvdls_mem.LS, cvdls_mem.A);\r
+     return(cvdls_mem.last_flag);\r
+   }\r
+   \r
+   private int SUNLinSolSetup_Dense(SUNLinearSolver S, double A[][])\r
+   {\r
+     double A_cols[][];\r
+     int pivots[];\r
+     int i,j;\r
+     \r
+     /* check for valid inputs */\r
+     if ( (A == null) || (S == null) ) \r
+       return(SUNLS_MEM_NULL);\r
+     \r
+     \r
+     /* access data pointers (return with failure on NULL) */\r
+     A_cols = null;\r
+     pivots = null;\r
+     //A_cols = SUNDenseMatrix_Cols(A);\r
+     // return SM_COLS_D(A)\r
+     // // #define SM_COLS_D(A)        ( SM_CONTENT_D(A)->cols )\r
+     A_cols = new double[A[0].length][A.length];\r
+     for (i = 0; i < A.length; i++) {\r
+        for (j = 0; j < A[0].length; j++) {\r
+                A_cols[j][i] = A[i][j];\r
+        }\r
+     }\r
+     pivots = S.pivots;\r
+     if ( (A_cols == null) || (pivots == null) ) {\r
+       S.last_flag = SUNLS_MEM_FAIL;\r
+       return(S.last_flag);\r
+     }\r
+     \r
+     /* perform LU factorization of input matrix */\r
+     S.last_flag = denseGETRF(A_cols, A.length,\r
+                              A[0].length, pivots);\r
+     for (i = 0; i < A.length; i++) {\r
+        for (j = 0; j < A[0].length; j++) {\r
+                A[i][j] = A_cols[j][i];\r
+        }\r
+     }\r
+\r
+     /* store error flag (if nonzero, this row encountered zero-valued pivod) */\r
+     if (S.last_flag > 0)\r
+       return(SUNLS_LUFACT_FAIL);\r
+     return(SUNLS_SUCCESS);\r
+   }\r
+   \r
+   private int denseGETRF(double a[][], int m, int n, int p[])\r
+   {\r
+     int i, j, k, l, q;\r
+     // a is [n][m]\r
+     // n column each with m row numbers\r
+     double col_j[] = new double[m];\r
+     double col_k[] = new double[m];\r
+     double temp, mult, a_kj;\r
+\r
+     /* k-th elimination step number */\r
+     for (k=0; k < n; k++) {\r
+\r
+         for (q = 0; q < m; q++) {\r
+            col_k[q]  = a[k][q];\r
+         }\r
+\r
+       /* find l = pivot row number */\r
+       l=k;\r
+       for (i=k+1; i < m; i++)\r
+         if (Math.abs(col_k[i]) > Math.abs(col_k[l])) l=i;\r
+       p[k] = l;\r
+\r
+       /* check for zero pivot element */\r
+       if (col_k[l] == ZERO) return(k+1);\r
+       \r
+       /* swap a(k,1:n) and a(l,1:n) if necessary */    \r
+       if ( l!= k ) {\r
+         for (i=0; i<n; i++) {\r
+           temp = a[i][l];\r
+           a[i][l] = a[i][k];\r
+           a[i][k] = temp;\r
+         }\r
+       }\r
+\r
+       /* Scale the elements below the diagonal in\r
+        * column k by 1.0/a(k,k). After the above swap\r
+        * a(k,k) holds the pivot element. This scaling\r
+        * stores the pivot row multipliers a(i,k)/a(k,k)\r
+        * in a(i,k), i=k+1, ..., m-1.                      \r
+        */\r
+       mult = ONE/col_k[k];\r
+       for(i=k+1; i < m; i++) col_k[i] *= mult;\r
+\r
+       /* row_i = row_i - [a(i,k)/a(k,k)] row_k, i=k+1, ..., m-1 */\r
+       /* row k is the pivot row after swapping with row l.      */\r
+       /* The computation is done one column at a time,          */\r
+       /* column j=k+1, ..., n-1.                                */\r
+\r
+       for (j=k+1; j < n; j++) {\r
+         for (q = 0; q < m; q++) {\r
+             col_j[q] = a[j][q];\r
+         }\r
+         a_kj = col_j[k];\r
+\r
+         /* a(i,j) = a(i,j) - [a(i,k)/a(k,k)]*a(k,j)  */\r
+         /* a_kj = a(k,j), col_k[i] = - a(i,k)/a(k,k) */\r
+\r
+         if (a_kj != ZERO) {\r
+       for (i=k+1; i < m; i++)\r
+         col_j[i] -= a_kj * col_k[i];\r
+         }\r
+       }\r
+     }\r
+\r
+     /* return 0 to indicate success */\r
+\r
+     return(0);\r
+   }\r
+\r
+   \r
+   private int SUNMatCopy(double A[][], double B[][]) {\r
+          int i,j;\r
+          if ((A == null) || (B == null) || (A.length != B.length) || (A[0].length != B[0].length)) {\r
+                  return 1;\r
+          }\r
+          for (i = 0; i < A.length; i++) {\r
+                  for (j = 0; j < A[0].length; j++) {\r
+                          B[i][j] = A[i][j];\r
+                  }\r
+          }\r
+          return 0;\r
+   }\r
+   \r
+   private int SUNMatZero(double A[][]) {\r
+          int i, j;\r
+          if (A == null) {\r
+                  return 1;\r
+          }\r
+          for (i = 0; i < A.length; i++) {\r
+                  for (j = 0; j < A[0].length; j++) {\r
+                          A[i][j] = 0.0;\r
+                  }\r
+          }\r
+          return 0;\r
+   }\r
+   \r
+   private int SUNMatScaleAddI(double c, double A[][]) {\r
+          int i, j;\r
+          if ((A == null) || (A.length != A[0].length)) {\r
+                  return 1;\r
+          }\r
+          for (i = 0; i < A.length; i++) {\r
+                  for (j = 0; j < A[0].length; j++) {\r
+                          A[i][j] = c * A[i][j];\r
+                          if (i == j) {\r
+                                  A[i][j] = A[i][j] + 1.0;\r
+                          }\r
+                  }\r
+          }\r
+          return 0;\r
+   }\r
+\r
+   /*\r
+    * cvNewtonIteration\r
+    *\r
+    * This routine performs the Newton iteration. If the iteration succeeds,\r
+    * it returns the value CV_SUCCESS. If not, it may signal the cvNlsNewton \r
+    * routine to call lsetup again and reattempt the iteration, by\r
+    * returning the value TRY_AGAIN. (In this case, cvNlsNewton must set \r
+    * convfail to CV_FAIL_BAD_J before calling setup again). \r
+    * Otherwise, this routine returns one of the appropriate values \r
+    * CV_LSOLVE_FAIL, CV_RHSFUNC_FAIL, CV_SRHSFUNC_FAIL, CONV_FAIL,\r
+    * RHSFUNC_RECVR, or SRHSFUNC_RECVR back to cvNlsNewton.\r
+    *\r
+    * If sensitivities are computed using the CV_SIMULTANEOUS approach, this\r
+    * routine performs a quasi-Newton on the combined nonlinear system.\r
+    * The iteration matrix of the combined system is block diagonal with\r
+    * each block being the iteration matrix of the original system. Thus\r
+    * we solve linear systems with the same matrix but different right\r
+    * hand sides.\r
+    */\r
+\r
+   private int cvNewtonIteration(CVodeMemRec cv_mem)\r
+   {\r
+     int m;\r
+     double del, delS, Del, Delp, dcon;\r
+     NVector bS[] = null;\r
+     NVector b, wrk1, wrk2;\r
+     boolean do_sensi_sim;\r
+     int retval, is;\r
+     \r
+     /* Are we computing sensitivities with the CV_SIMULTANEOUS approach? */\r
+     do_sensi_sim = (cv_mem.cv_sensi && (cv_mem.cv_ism==CV_SIMULTANEOUS));\r
+\r
+     cv_mem.cv_mnewt = m = 0;\r
+\r
+     /* Initialize delS and Delp to avoid compiler warning message */\r
+     delS = Delp = ZERO;\r
+\r
+     /* At this point, ftemp  <- f(t_n, y_n(0))\r
+        ftempS <- fS(t_n, y_n(0), s_n(0))\r
+        acor   <- 0\r
+        acorS  <- 0\r
+        y      <- y_n(0)\r
+        yS     <- yS_n(0)                 */\r
+\r
+     /* Looping point for Newton iteration */\r
+     for(;;) {\r
+       \r
+       /* Evaluate the residual of the nonlinear system */\r
+       N_VLinearSum_Serial(cv_mem.cv_rl1, cv_mem.cv_zn[1], ONE,\r
+                    cv_mem.cv_acor, cv_mem.cv_tempv);\r
+       N_VLinearSum_Serial(cv_mem.cv_gamma, cv_mem.cv_ftemp, -ONE,\r
+                    cv_mem.cv_tempv, cv_mem.cv_tempv);\r
+\r
+       /* Call the lsolve function */\r
+       b = cv_mem.cv_tempv;\r
+       /*retval = cv_mem->cv_lsolve(cv_mem, b, cv_mem.cv_ewt,\r
+                                  cv_mem.cv_y, cv_mem.cv_ftemp); \r
+       cv_mem->cv_nni++;\r
+\r
+       if (retval < 0) return(CV_LSOLVE_FAIL);*/\r
+       \r
+       /* If lsolve had a recoverable failure and Jacobian data is\r
+          not current, signal to try the solution again            */\r
+      /* if (retval > 0) { \r
+         if ((!cv_mem->cv_jcur) && (cv_mem->cv_lsetup))\r
+           return(TRY_AGAIN);\r
+         else\r
+           return(CONV_FAIL);\r
+       }*/\r
+\r
+       /* Solve the sensitivity linear systems and do the same \r
+          tests on the return value of lsolve. */\r
+    \r
+      /* if (do_sensi_sim) {\r
+\r
+         for (is=0; is<cv_mem->cv_Ns; is++) {\r
+           N_VLinearSum(cv_mem->cv_rl1, cv_mem->cv_znS[1][is], ONE,\r
+                        cv_mem->cv_acorS[is], cv_mem->cv_tempvS[is]);\r
+           N_VLinearSum(cv_mem->cv_gamma, cv_mem->cv_ftempS[is], -ONE,\r
+                        cv_mem->cv_tempvS[is], cv_mem->cv_tempvS[is]);\r
+         }\r
+         bS = cv_mem->cv_tempvS;\r
+         for (is=0; is<cv_mem->cv_Ns; is++) {\r
+           retval = cv_mem->cv_lsolve(cv_mem, bS[is], cv_mem->cv_ewtS[is],\r
+                                      cv_mem->cv_y, cv_mem->cv_ftemp);\r
+           if (retval < 0) return(CV_LSOLVE_FAIL);\r
+           if (retval > 0) { \r
+             if ((!cv_mem->cv_jcur) && (cv_mem->cv_lsetup))\r
+               return(TRY_AGAIN);\r
+             else\r
+               return(CONV_FAIL);\r
+           }\r
+         }\r
+       }*/\r
+       \r
+       /* Get WRMS norm of correction; add correction to acor and y */\r
+\r
+       /*del = N_VWrmsNorm(b, cv_mem->cv_ewt);\r
+       N_VLinearSum(ONE, cv_mem->cv_acor, ONE, b, cv_mem->cv_acor);\r
+       N_VLinearSum(ONE, cv_mem->cv_zn[0], ONE, cv_mem->cv_acor, cv_mem->cv_y);\r
+\r
+       if (do_sensi_sim) {\r
+         delS = cvSensUpdateNorm(cv_mem, del, bS, cv_mem->cv_ewtS);\r
+         for (is=0; is<cv_mem->cv_Ns; is++) {\r
+           N_VLinearSum(ONE, cv_mem->cv_acorS[is], ONE,\r
+                        bS[is], cv_mem->cv_acorS[is]);\r
+           N_VLinearSum(ONE, cv_mem->cv_znS[0][is], ONE,\r
+                        cv_mem->cv_acorS[is], cv_mem->cv_yS[is]);\r
+         }\r
+       }*/\r
+\r
+       /* Test for convergence.  If m > 0, an estimate of the convergence\r
+          rate constant is stored in crate, and used in the test.        */\r
+\r
+       /*Del = (do_sensi_sim) ? delS : del;\r
+       if (m > 0)\r
+         cv_mem->cv_crate = SUNMAX(CRDOWN * cv_mem->cv_crate, Del/Delp);\r
+       dcon = Del * SUNMIN(ONE, cv_mem->cv_crate) / cv_mem->cv_tq[4];\r
+       \r
+       if (dcon <= ONE) {\r
+         if (m == 0)\r
+           if (do_sensi_sim && cv_mem->cv_errconS)\r
+             cv_mem->cv_acnrm = delS;\r
+           else\r
+             cv_mem->cv_acnrm = del;\r
+         else {\r
+           cv_mem->cv_acnrm = N_VWrmsNorm(cv_mem->cv_acor, cv_mem->cv_ewt);\r
+           if (do_sensi_sim && cv_mem->cv_errconS)\r
+             cv_mem->cv_acnrm = cvSensUpdateNorm(cv_mem, cv_mem->cv_acnrm,\r
+                                                 cv_mem->cv_acorS, cv_mem->cv_ewtS);\r
+         }\r
+         cv_mem->cv_jcur = SUNFALSE;\r
+         return(CV_SUCCESS); */ /* Convergence achieved */\r
+       /*}\r
+\r
+       cv_mem->cv_mnewt = ++m;*/\r
+       \r
+       /* Stop at maxcor iterations or if iter. seems to be diverging.\r
+          If still not converged and Jacobian data is not current, \r
+          signal to try the solution again                            */\r
+       /*if ((m == cv_mem->cv_maxcor) || ((m >= 2) && (Del > RDIV * Delp))) {\r
+         if ((!cv_mem->cv_jcur) && (cv_mem->cv_lsetup))\r
+           return(TRY_AGAIN);\r
+         else\r
+           return(CONV_FAIL);\r
+       }*/\r
+       \r
+       /* Save norm of correction, evaluate f, and loop again */\r
+      /* Delp = Del;\r
+       retval = cv_mem->cv_f(cv_mem->cv_tn, cv_mem->cv_y,\r
+                             cv_mem->cv_ftemp, cv_mem->cv_user_data);\r
+       cv_mem->cv_nfe++;\r
+       if (retval < 0) return(CV_RHSFUNC_FAIL);\r
+       if (retval > 0) {\r
+         if ((!cv_mem->cv_jcur) && (cv_mem->cv_lsetup))\r
+           return(TRY_AGAIN);\r
+         else\r
+           return(RHSFUNC_RECVR);\r
+       }\r
+\r
+       if (do_sensi_sim) {\r
+         wrk1 = cv_mem->cv_tempv;\r
+         wrk2 = cv_mem->cv_tempvS[0];\r
+         retval = cvSensRhsWrapper(cv_mem, cv_mem->cv_tn, cv_mem->cv_y,\r
+                                   cv_mem->cv_ftemp, cv_mem->cv_yS,\r
+                                   cv_mem->cv_ftempS, wrk1, wrk2);\r
+         if (retval < 0) return(CV_SRHSFUNC_FAIL);\r
+         if (retval > 0) {\r
+           if ((!cv_mem->cv_jcur) && (cv_mem->cv_lsetup))\r
+             return(TRY_AGAIN);\r
+           else\r
+             return(SRHSFUNC_RECVR);\r
+         }\r
+       }*/\r
+\r
+     } /* end loop */\r
+\r
+   }\r
 \r
 \r
 }
\ No newline at end of file