cvsRoberts_dnsL example added which uses LAPACK dgetrf and dgetrs instead of denseGET...
authorilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Wed, 28 Mar 2018 18:54:25 +0000 (18:54 +0000)
committerilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Wed, 28 Mar 2018 18:54:25 +0000 (18:54 +0000)
git-svn-id: https://citdcbmipav.cit.nih.gov/repos-pub/mipav/trunk@15432 ba61647d-9d00-f842-95cd-605cb4296b96

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

index 1f18483..c1cc127 100644 (file)
@@ -2,6 +2,7 @@ package gov.nih.mipav.model.algorithms;
 \r
 import java.io.RandomAccessFile;\r
 \r
+import gov.nih.mipav.model.structures.jama.LinearEquations2;\r
 import gov.nih.mipav.view.MipavUtil;\r
 \r
 public abstract class CVODES {\r
@@ -545,7 +546,8 @@ public abstract class CVODES {
        public enum SUNLinearSolver_Type{\r
                  SUNLINEARSOLVER_DIRECT,\r
                  SUNLINEARSOLVER_ITERATIVE,\r
-                 SUNLINEARSOLVER_CUSTOM\r
+                 SUNLINEARSOLVER_CUSTOM,\r
+                 SUNLINEARLAPACKSOLVER_DIRECT\r
                };\r
 \r
     final int cvDlsDQJac = -1;\r
@@ -557,10 +559,12 @@ public abstract class CVODES {
     final int cvDlsInitialize_select = 1;\r
     // For cv_mem.cv_lsolve_select\r
     final int cvDlsSolve_select = 1;\r
+    final int cvDlsLapackSolve_select = 2;\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
+    final int cvDlsLapackSetup_select = 2;\r
     // For cv_mem.cv_fS1_select\r
     final int cvSensRhs1InternalDQ_select = -1;\r
     // For cv_mem.cv_f_select\r
@@ -569,8 +573,9 @@ public abstract class CVODES {
     final int cvsRoberts_dns = 1;\r
     final int cvsDirectDemo_ls_Problem_1 = 2;\r
     final int cvsRoberts_dns_uw = 3;\r
-    final int cvsRoberts_FSA_dns = 4;\r
-    int problem = cvsDirectDemo_ls_Problem_1;\r
+    final int cvsRoberts_dnsL = 4;\r
+    final int cvsRoberts_FSA_dns = 5;\r
+    int problem = cvsRoberts_dnsL;\r
     boolean testMode = true;\r
        \r
     // Linear Solver options for runcvsDirectDemo\r
@@ -655,6 +660,9 @@ public abstract class CVODES {
            else if (problem == cvsRoberts_dns_uw) {\r
                runcvsRoberts_dns_uw();\r
            }\r
+           else if (problem == cvsRoberts_dnsL) {\r
+               runcvsRoberts_dnsL();   \r
+           }\r
            else if (problem == cvsRoberts_FSA_dns) {\r
                runcvsRoberts_FSA_dns();\r
            }\r
@@ -1169,6 +1177,237 @@ public abstract class CVODES {
                return;\r
        }\r
        \r
+       private void runcvsRoberts_dnsL() {\r
+               /* -----------------------------------------------------------------\r
+                * Programmer(s): Radu Serban @ LLNL\r
+                * -----------------------------------------------------------------\r
+                * Example problem:\r
+                * \r
+                * The following is a simple example problem, with the coding\r
+                * needed for its solution by CVODE. The problem is from\r
+                * chemical kinetics, and consists of the following three rate\r
+                * equations:         \r
+                *    dy1/dt = -.04*y1 + 1.e4*y2*y3\r
+                *    dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e7*(y2)^2\r
+                *    dy3/dt = 3.e7*(y2)^2\r
+                * on the interval from t = 0.0 to t = 4.e10, with initial\r
+                * conditions: y1 = 1.0, y2 = y3 = 0. The problem is stiff.\r
+                * While integrating the system, we also use the rootfinding\r
+                * feature to find the points at which y1 = 1e-4 or at which\r
+                * y3 = 0.01. This program solves the problem with the BDF method,\r
+                * Newton iteration with the LAPACK dense linear solver, and a\r
+                * user-supplied Jacobian routine.\r
+                * It uses a scalar relative tolerance and a vector absolute\r
+                * tolerance. Output is printed in decades from t = .4 to t = 4.e10.\r
+                * Run statistics (optional outputs) are printed at the end.\r
+                * -----------------------------------------------------------------*/\r
+               /** Problem Constants */\r
+               final int NEQ = 3; // Number of equations\r
+               final double Y1 = 1.0; // Initial y components\r
+               final double Y2 = 0.0;\r
+               final double Y3 = 0.0;\r
+               //final double RTOL = 1.0E-4; // scalar relative tolerance\r
+               final double RTOL = 1.0E-12;\r
+               //final double ATOL1 = 1.0E-8; // vector absolute tolerance components\r
+               final double ATOL1 = 1.0E-12;\r
+               //final double ATOL2 = 1.0E-14;\r
+               final double ATOL2 = 1.0E-15;\r
+               //final double ATOL3 = 1.0E-6;\r
+               final double ATOL3 = 1.0E-12;\r
+               final double T0 = 0.0; // initial time\r
+               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
+               final int NOUT = 12;\r
+               int f = cvsRoberts_dnsL;\r
+               int g = cvsRoberts_dnsL;\r
+               int Jac = cvsRoberts_dnsL;\r
+               \r
+               // Set initial conditions\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 abstolr[] = new double[]{ATOL1,ATOL2,ATOL3};\r
+               N_VNew_Serial(abstol,NEQ);\r
+               abstol.setData(abstolr);\r
+               CVodeMemRec cvode_mem;\r
+               int flag;\r
+               int flagr;\r
+               double A[][];\r
+               SUNLinearSolver LS;\r
+               double tout;\r
+               int iout;\r
+               double t[] = new double[1];\r
+               int rootsfound[] = new int[2];\r
+               int i;\r
+               \r
+               // Call CVodeCreate to create the solver memory and specify the\r
+               // Backward Differentiation Formula and the use of a Newton\r
+               // iteration\r
+               cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);\r
+               if (cvode_mem == null) {\r
+                   return;     \r
+               }\r
+               // Allow unlimited steps in reaching tout\r
+               flag = CVodeSetMaxNumSteps(cvode_mem, -1);\r
+               if (flag != CV_SUCCESS) {\r
+                       return;\r
+               }\r
+               // Allow many error test failures\r
+               flag = CVodeSetMaxErrTestFails(cvode_mem, Integer.MAX_VALUE);\r
+               if (flag != CV_SUCCESS) {\r
+                       return;\r
+               }\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
+               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 = SUNDenseLinearLapackSolver(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
+               if (flag != CVDLS_SUCCESS) {\r
+                       return;\r
+               }\r
+               \r
+               // Set the user-supplied Jacobian routine Jac\r
+               flag = CVDlsSetJacFn(cvode_mem, Jac);\r
+               if (flag != CVDLS_SUCCESS) {\r
+                       return;\r
+               }\r
+               \r
+               // In loop, call CVode, print results, and test for error.\r
+               // Break out of loop when NOUT preset output times have been reached.\r
+               System.out.println(" \n3-species kinetics problem\n");\r
+               \r
+               iout = 0;\r
+               tout = T1;\r
+               \r
+               while (true) {\r
+                       flag = CVode(cvode_mem, tout, y, t, CV_NORMAL);\r
+                       switch (iout) {\r
+                       case 0:\r
+                               System.out.println("Example manual: at t = 0.4 y[0] = 0.98517 y[1] = 3.3864E-5 y[2] = 1.4794E-2");\r
+                               break;\r
+                       case 1:\r
+                               System.out.println("Example manual: at t = 4.0 y[0] = 0.90552 y[1] = 2.2405E-5 y[2] = 9.4459E-2");\r
+                               break;\r
+                       case 2:\r
+                               System.out.println("Example manual: at t = 40.0 y[0] = 0.71583 y[1] = 9.1856E-6 y[2] = 0.28416");\r
+                               break;\r
+                       case 3:\r
+                               System.out.println("Example manual: at t = 400.0 y[0] = 0.45052 y[1] = 3.2229E-6 y[2] = 0.54947");\r
+                               break;\r
+                       case 4:\r
+                               System.out.println("Example manual: at t = 4.0E3 y[0] = 0.18317 y[1] = 8.9403E-7 y[2] = 0.81683");\r
+                               break;\r
+                       case 5:\r
+                               System.out.println("Example manual: at t = 4.0E4 y[0] = 3.8977E-2 y[1] = 1.6215E-7 y[2] = 0.96102");\r
+                               break;\r
+                       case 6:\r
+                               System.out.println("Example manual: at t = 4.0E5 y[0] = 4.9387E-3 y[1] = 1.9852E-8 y[2] = 0.99506");\r
+                               break;\r
+                       case 7:\r
+                               System.out.println("Example manual: at t = 4.0E6 y[0] = 5.1684E-4 y[1] = 2.0684E-9 y[2] = 0.99948");\r
+                               break;\r
+                       case 8:\r
+                               System.out.println("Example manual: at t = 4.0E7 y[0] = 5.2039E-5 y[1] = 2.0817E-10 y[2] = 0.99995");\r
+                               break;\r
+                       case 9:\r
+                               System.out.println("Example manual at t = 4.0E8 y[0] = 5.2106E-6 y[1] = 2.0842E-11 y[2] = 0.99999");\r
+                               break;\r
+                       case 10:\r
+                               System.out.println("Example manual at t = 4.0E9 y[0] = 5.1881E-7 y[1] = 2.0752E-12 y[2] = 1.0000");\r
+                               break;\r
+                       case 11:\r
+                               System.out.println("Example manual at t = 4.0E10 y[0] = 6.5181E-8 y[1] = 2.6072E-13 y[2] = 1.0000");\r
+                               break;\r
+                       }\r
+                       System.out.println("MIPAV: at t = " + t[0] + " y[0] = " + y.data[0] + " y[1] = " + y.data[1] + " y[2] = " + y.data[2]);\r
+                       \r
+                       if (flag == CV_ROOT_RETURN) {\r
+                           flagr = CVodeGetRootInfo(cvode_mem, rootsfound)     ;\r
+                           if (flagr == CV_MEM_NULL) {\r
+                               return;\r
+                           }\r
+                           System.out.println("Roots found are " + rootsfound[0] + " and " + rootsfound[1]);\r
+                       }\r
+                       \r
+                       if (flag < 0) {\r
+                               System.out.println("CVode failed with flag = " + flag);\r
+                               break;\r
+                       }\r
+                       \r
+                       if (flag == CV_SUCCESS) {\r
+                               iout++;\r
+                               tout *= TMULT;\r
+                       }\r
+                       \r
+                       if (iout == NOUT) {\r
+                               break;\r
+                       }\r
+               } // while (true)\r
+               \r
+               // Print some final statistics\r
+               PrintFinalStats(cvode_mem);\r
+               \r
+               // Check the solution error\r
+               flag = check_ans(y, reltol, abstol);\r
+               \r
+               // Free y and abstol vectors\r
+               N_VDestroy(y);\r
+               N_VDestroy(abstol);\r
+               \r
+               // Free the integrator memory\r
+               CVodeFree(cvode_mem);\r
+               \r
+               // Free the linear solver memory\r
+               SUNLinSolFree_Dense(LS);\r
+               \r
+               // Free the matrix memory\r
+               for (i = 0; i < NEQ; i++) {\r
+                       A[i] = null;\r
+               }\r
+               A = null;\r
+               return;\r
+       }\r
+       \r
        private void runcvsRoberts_FSA_dns() {\r
                /* -----------------------------------------------------------------\r
                 * Programmer(s): Scott D. Cohen, Alan C. Hindmarsh, and\r
@@ -1950,7 +2189,7 @@ public abstract class CVODES {
        private int fTestMode(double t, NVector yv, NVector ydotv, UserData user_data) {\r
                double y[] = yv.getData();\r
                double ydot[] = ydotv.getData();\r
-               if ((problem == cvsRoberts_dns) || (problem == cvsRoberts_dns_uw)) {\r
+               if ((problem == cvsRoberts_dns) || (problem == cvsRoberts_dns_uw) || (problem == cvsRoberts_dnsL)) {\r
                    ydot[0] = -0.04*y[0] + 1.0E4*y[1]*y[2];\r
                    ydot[2] = 3.0E7*y[1]*y[1];\r
                    ydot[1] = -ydot[0] - ydot[2];\r
@@ -2019,7 +2258,7 @@ public abstract class CVODES {
         */\r
        private int gTestMode(double t, NVector yv, double gout[], UserData user_data) {\r
                double y[] = yv.getData();\r
-               if ((problem == cvsRoberts_dns) || (problem == cvsRoberts_dns_uw) || (problem == cvsRoberts_FSA_dns)) {\r
+               if ((problem == cvsRoberts_dns) || (problem == cvsRoberts_dns_uw) || (problem == cvsRoberts_dnsL) || (problem == cvsRoberts_FSA_dns)) {\r
                    gout[0] = y[0] - 0.0001;\r
                    gout[1] = y[2] - 0.01;\r
                }\r
@@ -2041,7 +2280,7 @@ public abstract class CVODES {
         */\r
        private int JacTestMode(double t, NVector yv, NVector fy, double J[][], UserData data, NVector tmp1, NVector tmp2, NVector tmp3) {\r
                double y[] = yv.getData();\r
-               if ((problem == cvsRoberts_dns) || (problem == cvsRoberts_dns_uw)) {\r
+               if ((problem == cvsRoberts_dns) || (problem == cvsRoberts_dns_uw) || (problem == cvsRoberts_dnsL)) {\r
                    J[0][0] = -0.04;\r
                    J[0][1] = 1.0E4 * y[2];\r
                    J[0][2] = 1.0E4 * y[1];\r
@@ -3279,6 +3518,26 @@ public abstract class CVODES {
         return S;\r
      }\r
      \r
+     private SUNLinearSolver SUNDenseLinearLapackSolver(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 int[matrixRows];\r
+        S.type = SUNLinearSolver_Type.SUNLINEARLAPACKSOLVER_DIRECT;\r
+        return S;\r
+     }\r
+     \r
      /*---------------------------------------------------------------\r
      CVDlsSetLinearSolver specifies the direct linear solver.\r
     ---------------------------------------------------------------*/\r
@@ -3301,7 +3560,8 @@ public abstract class CVODES {
       }\r
 \r
       /* Test if solver and vector are compatible with DLS */\r
-      if (LS.type != SUNLinearSolver_Type.SUNLINEARSOLVER_DIRECT) {\r
+      if ((LS.type != SUNLinearSolver_Type.SUNLINEARSOLVER_DIRECT) &&\r
+        (LS.type != SUNLinearSolver_Type.SUNLINEARLAPACKSOLVER_DIRECT)) {\r
         cvProcessError(cv_mem, CVDLS_ILL_INPUT, "CVSDLS", \r
                        "CVDlsSetLinearSolver", \r
                        "Non-direct LS supplied to CVDls interface");\r
@@ -3321,9 +3581,19 @@ public abstract class CVODES {
       //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
+      if (LS.type == SUNLinearSolver_Type.SUNLINEARSOLVER_DIRECT) {\r
+          cv_mem.cv_lsetup_select = cvDlsSetup_select;\r
+      }\r
+      else {\r
+         cv_mem.cv_lsetup_select = cvDlsLapackSetup_select;\r
+      }\r
       //cv_mem.cv_lsolve = cvDlsSolve;\r
-      cv_mem.cv_lsolve_select = cvDlsSolve_select;\r
+      if (LS.type == SUNLinearSolver_Type.SUNLINEARSOLVER_DIRECT) {\r
+          cv_mem.cv_lsolve_select = cvDlsSolve_select;\r
+      }\r
+      else {\r
+         cv_mem.cv_lsolve_select = cvDlsLapackSolve_select;\r
+      }\r
       //cv_mem.cv_lfree  = cvDlsFree;\r
       cv_mem.cv_lfree_select = cvDlsFree_select;\r
       \r
@@ -7465,8 +7735,9 @@ else                return(snrm);
         int flag = 0;\r
         switch(select) {\r
         case cvDlsSetup_select:\r
+        case cvDlsLapackSetup_select:\r
                 flag = cvDlsSetup(cv_mem, convfail, y, fy, jcurPtr,\r
-                                tmp1, tmp2, tmp3);\r
+                                tmp1, tmp2, tmp3,select);\r
         }\r
         return flag;\r
      }\r
@@ -7484,7 +7755,7 @@ else                return(snrm);
      -----------------------------------------------------------------*/\r
    int cvDlsSetup(CVodeMemRec cv_mem, int convfail, NVector y, \r
                   NVector fy, boolean jcurPtr[], \r
-                  NVector tmp1, NVector tmp2, NVector tmp3)\r
+                  NVector tmp1, NVector tmp2, NVector tmp3, int select)\r
    {\r
      boolean jbad, jok;\r
      double dgamma;\r
@@ -7580,7 +7851,12 @@ else                return(snrm);
 \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
+     if (select == cvDlsSetup_select) {\r
+         cvdls_mem.last_flag = SUNLinSolSetup_Dense(cvdls_mem.LS, cvdls_mem.A);\r
+     }\r
+     else {\r
+        cvdls_mem.last_flag = SUNLinSolSetup_LapackDense(cvdls_mem.LS, cvdls_mem.A);\r
+     }\r
      return(cvdls_mem.last_flag);\r
    }\r
    \r
@@ -7739,6 +8015,33 @@ else                return(snrm);
           }\r
           return 0;\r
    }\r
+   \r
+   private int SUNLinSolSetup_LapackDense(SUNLinearSolver S, double A[][]) \r
+   {\r
+     int pivots[];\r
+     int n;\r
+     LinearEquations2 le2 = new LinearEquations2();\r
+     \r
+     /* check for valid inputs */\r
+     if ( (A == null) || (S == null) ) \r
+       return(SUNLS_MEM_NULL);\r
+     \r
+     n = A.length;\r
+     int ier[] = new int[1];\r
+     pivots = S.pivots;\r
+     le2.dgetrf(n, n, A, n, pivots, ier);\r
+        //if (ier[0] > 0) {\r
+           //MipavUtil.displayError("In dgetrf factor U is exactly singular");\r
+        //}\r
+        S.last_flag = ier[0];\r
+        if (ier[0] > 0) {\r
+                return(SUNLS_LUFACT_FAIL);\r
+        }\r
+        if (ier[0] < 0) {\r
+                return(SUNLS_PACKAGE_FAIL_UNREC);\r
+        }\r
+        return(SUNLS_SUCCESS);\r
+   }\r
 \r
    /*\r
     * cvNewtonIteration\r
@@ -7932,7 +8235,8 @@ else                return(snrm);
           int flag = 0;\r
           switch (select) {\r
           case cvDlsSolve_select:\r
-                  flag = cvDlsSolve(cv_mem, b, weight, ycur, fcur);\r
+          case cvDlsLapackSolve_select:\r
+                  flag = cvDlsSolve(cv_mem, b, weight, ycur, fcur, select);\r
                   break;\r
           }\r
           return flag;\r
@@ -7946,7 +8250,7 @@ else                return(snrm);
    the solution appropriately when gamrat != 1.\r
    -----------------------------------------------------------------*/\r
  private int cvDlsSolve(CVodeMemRec cv_mem, NVector b, NVector weight,\r
-                NVector ycur, NVector fcur)\r
+                NVector ycur, NVector fcur, int select)\r
  {\r
    int retval;\r
    CVDlsMemRec cvdls_mem;\r
@@ -7965,7 +8269,12 @@ else                return(snrm);
    cvdls_mem = cv_mem.cv_lmem;\r
 \r
    /* call the generic linear system solver, and copy b to x */\r
-   retval = SUNLinSolSolve_Dense(cvdls_mem.LS, cvdls_mem.A, cvdls_mem.x, b, ZERO);\r
+   if (select == cvDlsSolve_select) {\r
+       retval = SUNLinSolSolve_Dense(cvdls_mem.LS, cvdls_mem.A, cvdls_mem.x, b, ZERO);\r
+   }\r
+   else {\r
+          retval = SUNLinSolSolve_LapackDense(cvdls_mem.LS, cvdls_mem.A, cvdls_mem.x, b, ZERO);\r
+   }\r
    N_VScale_Serial(ONE, cvdls_mem.x, b);\r
    \r
    /* scale the correction to account for change in gamma */\r
@@ -8055,6 +8364,45 @@ else                return(snrm);
           b[0] /= a[0][0];\r
 \r
         }\r
+        \r
+        private int SUNLinSolSolve_LapackDense(SUNLinearSolver S, double A[][], NVector x, \r
+                NVector b, double tol)\r
+       {\r
+    double xdata[];\r
+       int pivots[];\r
+       int n;\r
+       int ier[] = new int[1];\r
+       double xarr[][];\r
+       LinearEquations2 le2 = new LinearEquations2();\r
+       int i;\r
+       \r
+       if ( (A == null) || (S == null) || (x == null) || (b == null) ) \r
+       return(SUNLS_MEM_NULL);\r
+       \r
+       /* copy b into x */\r
+       N_VScale_Serial(ONE, b, x);\r
+       \r
+       xdata = x.data;\r
+       xarr = new double[x.data.length][1];\r
+       for (i = 0; i < x.data.length; i++) {\r
+               xarr[i][0] = xdata[i];\r
+       }\r
+       pivots = S.pivots;\r
+       if ( (xdata == null)  || (pivots == null) ) {\r
+       S.last_flag = SUNLS_MEM_FAIL;\r
+       return(S.last_flag);\r
+       }\r
+       \r
+       // Call LAPACK to solve the linear system\r
+       n = A.length;\r
+       le2.dgetrs('N', n, 1, A, n, pivots, xarr, n, ier);\r
+       S.last_flag = ier[0];\r
+       if (ier[0] < 0) {\r
+               return(SUNLS_PACKAGE_FAIL_UNREC);\r
+       }\r
+       S.last_flag = SUNLS_SUCCESS;\r
+       return(S.last_flag);\r
+       }\r
 \r
         /*\r
          * cvHandleNFlag\r