Now call cvDlsDenseDQJac if a difference quotient Jacobian is used instead of a user...
authorilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Mon, 12 Mar 2018 20:45:20 +0000 (20:45 +0000)
committerilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Mon, 12 Mar 2018 20:45:20 +0000 (20:45 +0000)
git-svn-id: https://citdcbmipav.cit.nih.gov/repos-pub/mipav/trunk@15412 ba61647d-9d00-f842-95cd-605cb4296b96

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

index b728fb9..337f934 100644 (file)
@@ -2,6 +2,8 @@ package gov.nih.mipav.model.algorithms;
 \r
 import java.io.RandomAccessFile;\r
 \r
+import gov.nih.mipav.model.algorithms.CVODES.CVodeMemRec;\r
+import gov.nih.mipav.model.algorithms.CVODES.NVector;\r
 import gov.nih.mipav.view.MipavUtil;\r
 import gov.nih.mipav.view.Preferences;\r
 \r
@@ -259,6 +261,8 @@ public abstract class CVODES {
     // LONG_WAIT   number of steps to wait before considering an order change when\r
     //             q==1 and MXNEF1 error test failures have occurred\r
     final int LONG_WAIT = 10;\r
+    /* Constant for DQ Jacobian approximation */\r
+    final double MIN_INC_MULT = 1000.0;\r
 \r
 \r
        final int RTFOUND = +1;\r
@@ -559,7 +563,7 @@ public abstract class CVODES {
 \r
     final int cvsRoberts_dns = 1;\r
     final int cvSDirectDemo_ls_Problem_1 = 2;\r
-    int problem = cvsRoberts_dns;\r
+    int problem = cvSDirectDemo_ls_Problem_1;\r
     boolean testMode = true;\r
        \r
     // Linear Solver options for runcvsDirectDemo\r
@@ -578,6 +582,34 @@ public abstract class CVODES {
        final double P1_DTOUT = 2.214773875;\r
        final double P1_TOL_FACTOR = 1.0E4;\r
        \r
+       // USe the following code to call CVODES from another module:\r
+       /*boolean testme = true;\r
+    class CVODEStest extends CVODES {\r
+         public CVODEStest() {\r
+                 super();\r
+         }\r
+         public int f(double t, NVector yv, NVector ydotv, CVodeMemRec user_data) {\r
+                 return 0;\r
+         }\r
+         \r
+         public int g(double t, NVector yv, double gout[], CVodeMemRec user_data) {\r
+                 return 0;\r
+         }\r
+         \r
+         public int fQ(double t, NVector x, NVector y, CVodeMemRec user_data) {\r
+                 return 0;\r
+         }\r
+         \r
+         public int Jac(double t, NVector yv, NVector fy, double J[][], CVodeMemRec data, NVector tmp1, \r
+                         NVector tmp2, NVector tmp3) {\r
+                 return 0;\r
+         }\r
+       \r
+      }\r
+    if (testme) {\r
+       new CVODEStest();\r
+       return;\r
+    } */\r
        public CVODES() {\r
                // eps returns the distance from 1.0 to the next larger double-precision\r
                // number, that is, eps = 2^-52.\r
@@ -604,6 +636,9 @@ public abstract class CVODES {
            if (problem == cvsRoberts_dns) {\r
                runcvsRoberts_dns();\r
            }\r
+           else if (problem == cvSDirectDemo_ls_Problem_1) {\r
+               runcvsDirectDemo_Problem_1();\r
+           }\r
            return;\r
        }\r
        \r
@@ -1086,6 +1121,78 @@ public abstract class CVODES {
                    } // for (iout = 1, tout = P1_T1; iout <= P1_NOUT; iout++, tout += P1_DTOUT)\r
                    PrintFinalStats(cvode_mem, miter, ero);\r
                } // for (miter = FUNC; miter <= DENSE_DQ; miter++)\r
+               CVodeFree(cvode_mem);\r
+               \r
+               cvode_mem = CVodeCreate(CV_BDF, CV_FUNCTIONAL);\r
+               if (cvode_mem == null) {\r
+                       return;\r
+               }\r
+               for (miter = FUNC; miter <= DENSE_DQ; miter++) {\r
+                   ero = ZERO;\r
+                   y.data[0] = TWO;\r
+                   y.data[1] = ZERO;\r
+                   \r
+                   firstrun = (miter == FUNC);\r
+                   if (firstrun) {\r
+                       flag = CVodeInit(cvode_mem, f1, P1_T0, y);\r
+                       if (flag < 0) {\r
+                               return;\r
+                       }\r
+                       flag = CVodeSStolerances(cvode_mem, reltol, abstol);\r
+                       if (flag < 0) {\r
+                               return;\r
+                       }\r
+                   } // if (firstrun)\r
+                   else {\r
+                       flag = CVodeSetIterType(cvode_mem, CV_NEWTON);\r
+                       if (flag < 0) {\r
+                               return;\r
+                       }\r
+                       flag = CVodeReInit(cvode_mem, P1_T0, y);\r
+                       if (flag < 0) {\r
+                               return;\r
+                       }\r
+                   } // else\r
+                   \r
+                   flag = PrepareNextRun(cvode_mem, CV_BDF, miter, y, A, 0, 0, LS);\r
+                   if (flag < 0) {\r
+                       return;\r
+                   }\r
+                   \r
+                   for (iout = 1, tout = P1_T1; iout <= P1_NOUT; iout++, tout += P1_DTOUT) {\r
+                       flag = CVode(cvode_mem, tout, y, t, CV_NORMAL);\r
+                       if (flag < 0) {\r
+                               return;\r
+                       }\r
+                       // Returns the order on the last successful step\r
+                       qu = cvode_mem.cv_qu;\r
+                       // Return the step size used on the last successful step\r
+                       hu = cvode_mem.cv_hu;\r
+                       System.out.println("time = " + t[0]);\r
+                       System.out.println("y[0] = "+ y.data[0]);\r
+                       System.out.println("y[1] = " + y.data[1]);\r
+                       System.out.println("Step order qu = " + qu);\r
+                       System.out.println("Step size hu = " + hu);\r
+                       if (flag != CV_SUCCESS) {\r
+                               nerr++;\r
+                               break;\r
+                       }\r
+                       if ((iout %2) == 0) {\r
+                           er = Math.abs(y.data[0])/abstol;\r
+                           if (er > ero) ero = er;\r
+                           if (er > P1_TOL_FACTOR) {\r
+                               nerr++;\r
+                               System.out.println("Error exceeds " + P1_TOL_FACTOR + " * tolerance");\r
+                           }\r
+                       } // if ((iout %2) == 0)\r
+                       \r
+                   } // for (iout = 1, tout = P1_T1; iout <= P1_NOUT; iout++, tout += P1_DTOUT)\r
+                   PrintFinalStats(cvode_mem, miter, ero);\r
+               }\r
+               CVodeFree(cvode_mem);\r
+               N_VDestroy(y);\r
+               return;\r
+               \r
        }\r
        \r
        private int PrepareNextRun(CVodeMemRec cvode_mem, int lmm, int miter, NVector y, double A[][],\r
@@ -1157,7 +1264,12 @@ public abstract class CVODES {
        \r
        private void PrintFinalStats(CVodeMemRec cvode_mem, int miter, double ero) {\r
            long lenrw, leniw;\r
-           long nst, nfe, nsetups, netf, nni, ncfn, nje, nfeLS;\r
+           long nje = 0L;\r
+           long nfeLS = 0L;\r
+           long nfe,nst, nsetups, netf, nni, ncfn;\r
+           long lenrwLS[] = new long[1];\r
+           long leniwLS[] = new long[1];\r
+           int flag;\r
            // Integrator work space requirements\r
            leniw = cvode_mem.cv_liw;\r
            lenrw = cvode_mem.cv_lrw;\r
@@ -1190,10 +1302,19 @@ public abstract class CVODES {
                        nje = cvode_mem.cv_lmem.nje;    \r
                        // the number of calls to the ODE function needed for the DQ Jacobian approximation\r
                        nfeLS = cvode_mem.cv_lmem.nfeDQ;\r
+                       flag = CVDlsGetWorkSpace(cvode_mem, lenrwLS, leniwLS);\r
+                       if (flag != CVDLS_SUCCESS) {\r
+                               return;\r
+                       }\r
                } // if (miter != DIAG)\r
                else {\r
                } // else\r
+               System.out.println("Linear solver real workspace length = " + lenrwLS[0]);\r
+               System.out.println("Linear solver integer workspace length = " + leniwLS[0]);\r
+               System.out.println("Number of Jacobian evaluations = " + nje);\r
+               System.out.println("Number of f evaluations in linear solver = " + nfeLS);\r
            } // if (miter != FUNC)\r
+           System.out.println("Error overrun = " + ero);\r
        }\r
        \r
        private void N_VNew_Serial(NVector y, int length) {\r
@@ -1219,6 +1340,10 @@ public abstract class CVODES {
                    ydot[2] = 3.0E7*y[1]*y[1];\r
                    ydot[1] = -ydot[0] - ydot[2];\r
                }\r
+               else if (problem == cvSDirectDemo_ls_Problem_1) {\r
+                       ydot[0] = y[1];\r
+                       ydot[1] = (ONE - y[0]*y[0]) * P1_ETA * y[1] - y[0];\r
+               }\r
                return 0;\r
        }\r
        \r
@@ -1244,6 +1369,7 @@ public abstract class CVODES {
                    gout[0] = y[0] - 0.0001;\r
                    gout[1] = y[2] - 0.01;\r
                }\r
+               \r
                return 0;\r
        }\r
        \r
@@ -1274,6 +1400,12 @@ public abstract class CVODES {
                    J[2][1] = 6.0E7 * y[1];\r
                    J[2][2] = ZERO;\r
                }\r
+               else if (problem == cvSDirectDemo_ls_Problem_1) {\r
+                       J[0][0] = ZERO;\r
+                       J[0][1] = ONE;\r
+                       J[1][0] = -TWO * P1_ETA * y[0] * y[1] - ONE;\r
+                       J[1][1] = P1_ETA * (ONE - y[0]*y[0]);\r
+               }\r
            \r
            return 0;\r
        }\r
@@ -1287,7 +1419,7 @@ public abstract class CVODES {
        // 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
+       public class NVector {\r
                double data[];\r
                boolean own_data;\r
                \r
@@ -1300,7 +1432,7 @@ public abstract class CVODES {
        }\r
        \r
          \r
-       private class CVodeMemRec {\r
+       public class CVodeMemRec {\r
            \r
          double cv_uround;         /* machine unit roundoff                         */   \r
 \r
@@ -6631,7 +6763,10 @@ else                return(snrm);
          return(-1);\r
        }\r
 \r
-       if (testMode) {\r
+       if (cvdls_mem.jacDQ) {\r
+          cvDlsDenseDQJac(cv_mem.cv_tn, y, fy, cvdls_mem.A, cv_mem, tmp1);  \r
+       }\r
+       else if (testMode) {\r
           retval = JacTestMode(cv_mem.cv_tn, y, fy, cvdls_mem.A, \r
                    cvdls_mem.J_data, tmp1, tmp2, tmp3);   \r
        }\r
@@ -9799,7 +9934,7 @@ else                return(snrm);
 \r
          return(CVDLS_SUCCESS);\r
        }\r
-       \r
+\r
        private void N_VSpace_Serial(NVector v, int lrw[], int liw[])\r
        {\r
          lrw[0] = v.data.length;\r
@@ -9807,6 +9942,90 @@ else                return(snrm);
 \r
          return;\r
        }\r
+       \r
+       /*-----------------------------------------------------------------\r
+         cvDlsDenseDQJac \r
+         -----------------------------------------------------------------\r
+         This routine generates a dense difference quotient approximation \r
+         to the Jacobian of f(t,y). It assumes that a dense SUNMatrix is \r
+         stored column-wise, and that elements within each column are \r
+         contiguous. The address of the jth column of J is obtained via\r
+         the accessor function SUNDenseMatrix_Column, and this pointer \r
+         is associated with an N_Vector using the N_VSetArrayPointer\r
+         function.  Finally, the actual computation of the jth column of \r
+         the Jacobian is done with a call to N_VLinearSum.\r
+         -----------------------------------------------------------------*/ \r
+       private int cvDlsDenseDQJac(double t, NVector y, NVector fy, \r
+                           double Jac[][], CVodeMemRec cv_mem, NVector tmp1)\r
+       {\r
+         double fnorm, minInc, inc, inc_inv, yjsaved, srur;\r
+         double y_data[];\r
+         double ewt_data[];\r
+         NVector ftemp, jthCol;\r
+         int j, N, i;\r
+         int retval = 0;\r
+         CVDlsMemRec cvdls_mem;\r
+\r
+         /* access DlsMem interface structure */\r
+         cvdls_mem = cv_mem.cv_lmem;\r
+\r
+         /* access matrix dimension */\r
+         N = Jac.length;\r
+\r
+         /* Rename work vector for readibility */\r
+         ftemp = tmp1;\r
+\r
+         /* Create an empty vector for matrix column calculations */\r
+         jthCol = N_VCloneEmpty_Serial(tmp1);\r
+\r
+         /* Obtain pointers to the data for ewt, y */\r
+         ewt_data = cv_mem.cv_ewt.data;\r
+         y_data   = y.data;\r
+\r
+         /* Set minimum increment based on uround and norm of f */\r
+         srur = Math.sqrt(cv_mem.cv_uround);\r
+         fnorm = N_VWrmsNorm_Serial(fy, cv_mem.cv_ewt);\r
+         minInc = (fnorm != ZERO) ?\r
+           (MIN_INC_MULT * Math.abs(cv_mem.cv_h) * cv_mem.cv_uround * N * fnorm) : ONE;\r
+\r
+         for (j = 0; j < N; j++) {\r
+\r
+           /* Generate the jth col of J(tn,y) */\r
+\r
+          for (i = 0; i < Jac.length; i++) {\r
+                  jthCol.data[i] = Jac[i][j];\r
+          }\r
+\r
+           yjsaved = y_data[j];\r
+           inc = Math.max(srur*Math.abs(yjsaved), minInc/ewt_data[j]);\r
+           y_data[j] += inc;\r
+\r
+           retval = f(t, y, ftemp, cv_mem.cv_user_data);\r
+           cvdls_mem.nfeDQ++;\r
+           if (retval != 0) break;\r
+           \r
+           y_data[j] = yjsaved;\r
+\r
+           inc_inv = ONE/inc;\r
+           N_VLinearSum_Serial(inc_inv, ftemp, -inc_inv, fy, jthCol);\r
+\r
+           /* DENSE_COL(Jac,j) = N_VGetArrayPointer(jthCol); */  /* UNNECESSARY? */\r
+         }\r
+\r
+         N_VDestroy(jthCol);\r
+\r
+         return(retval);\r
+       }\r
+       \r
+       private NVector N_VCloneEmpty_Serial(NVector w)\r
+       {\r
+         NVector v = new NVector();\r
+         v.data = new double[w.data.length];\r
+         v.own_data = false;\r
+         return(v);\r
+       }\r
+\r
+\r
 \r
 \r
 }
\ No newline at end of file