# HG changeset patch # User Bill Greene # Date 1550229868 -3600 # Fri Feb 15 12:24:28 2019 +0100 # Node ID 4bf27c090f5695bcf545fc4af15e2a61a3941d46 # Parent f034b29320ad5034ad5c66480f64411e9d773440 Update DAE/IDE solvers to work with SUNDIALS 3 (bug #52475). * libinterp/dldfcn/__ode15__.cc : use SUNDIALS API version 3.x diff --git a/libinterp/dldfcn/__ode15__.cc b/libinterp/dldfcn/__ode15__.cc --- a/libinterp/dldfcn/__ode15__.cc +++ b/libinterp/dldfcn/__ode15__.cc @@ -1,6 +1,7 @@ /* Copyright (C) 2016-2019 Francesco Faccio +Copyright (C) 2019 William Greene This file is part of Octave. @@ -112,7 +113,8 @@ havejacsparse (false), mem (nullptr), num (), ida_fun (nullptr), ida_jac (nullptr), dfdy (nullptr), dfdyp (nullptr), spdfdy (nullptr), spdfdyp (nullptr), fun (nullptr), jacfun (nullptr), jacspfun (nullptr), - jacdcell (nullptr), jacspcell (nullptr) + jacdcell (nullptr), jacspcell (nullptr), + sunJacMatrix (nullptr), sunLinearSolver (nullptr) { } @@ -122,11 +124,17 @@ havejacsparse (false), mem (nullptr), num (), ida_fun (ida_fcn), ida_jac (nullptr), dfdy (nullptr), dfdyp (nullptr), spdfdy (nullptr), spdfdyp (nullptr), fun (daefun), jacfun (nullptr), jacspfun (nullptr), - jacdcell (nullptr), jacspcell (nullptr) + jacdcell (nullptr), jacspcell (nullptr), + sunJacMatrix (nullptr), sunLinearSolver (nullptr) { } - ~IDA (void) { IDAFree (&mem); } + ~IDA (void) + { + IDAFree (&mem); + SUNLinSolFree(sunLinearSolver); + SUNMatDestroy(sunJacMatrix); + } IDA& set_jacobian (octave_function *jac, DAEJacFuncDense j) @@ -184,7 +192,7 @@ static N_Vector ColToNVec (const ColumnVector& data, long int n); void - set_up (void); + set_up (const ColumnVector& y); void set_tolerance (ColumnVector& abstol, realtype reltol); @@ -199,25 +207,24 @@ void resfun_impl (realtype t, N_Vector& yy, N_Vector& yyp, N_Vector& rr); - static int - jacdense (long int Neq, realtype t, realtype cj, N_Vector yy, - N_Vector yyp, N_Vector, DlsMat JJ, void *user_data, + jacdense (realtype t, realtype cj, N_Vector yy, + N_Vector yyp, N_Vector, SUNMatrix JJ, void *user_data, N_Vector, N_Vector, N_Vector) { IDA *self = static_cast (user_data); - self->jacdense_impl (Neq, t, cj, yy, yyp, JJ); + self->jacdense_impl (t, cj, yy, yyp, JJ); return 0; } void - jacdense_impl (long int Neq, realtype t, realtype cj, - N_Vector& yy, N_Vector& yyp, DlsMat& JJ); + jacdense_impl (realtype t, realtype cj, + N_Vector& yy, N_Vector& yyp, SUNMatrix& JJ); # if defined (HAVE_SUNDIALS_IDAKLU) static int jacsparse (realtype t, realtype cj, N_Vector yy, N_Vector yyp, - N_Vector, SlsMat Jac, void *user_data, N_Vector, + N_Vector, SUNMatrix Jac, void *user_data, N_Vector, N_Vector, N_Vector) { IDA *self = static_cast (user_data); @@ -227,7 +234,7 @@ void jacsparse_impl (realtype t, realtype cj, N_Vector& yy, - N_Vector& yyp, SlsMat& Jac); + N_Vector& yyp, SUNMatrix& Jac); #endif void set_maxstep (realtype maxstep); @@ -291,6 +298,8 @@ DAEJacFuncSparse jacspfun; DAEJacCellDense jacdcell; DAEJacCellSparse jacspcell; + SUNMatrix sunJacMatrix; + SUNLinearSolver sunLinearSolver; }; int @@ -323,36 +332,61 @@ } void - IDA::set_up (void) + IDA::set_up (const ColumnVector& y) { + N_Vector yy = ColToNVec(y, num); + if (havejacsparse) { # if defined (HAVE_SUNDIALS_IDAKLU) - if (IDAKLU (mem, num, num*num, CSC_MAT) != 0) - error ("IDAKLU solver not initialized"); + + sunJacMatrix = SUNSparseMatrix (num, num, num*num, CSC_MAT); + if (! sunJacMatrix) + error ("Unable to create sparse Jacobian for Sundials"); + + sunLinearSolver = SUNKLU (yy, sunJacMatrix); + if (! sunLinearSolver) + error ("Unable to create KLU sparse solver"); - IDASlsSetSparseJacFn (mem, IDA::jacsparse); + if (IDADlsSetLinearSolver (mem, sunLinearSolver, sunJacMatrix)) + error ("Unable to set sparse linear solver"); + + IDADlsSetJacFn (mem, IDA::jacsparse); + # else - error ("IDAKLU is not available in this version of Octave"); + error ("SUNDIALS SUNLINSOL KLU is not available in this version of Octave"); # endif + } else { - if (IDADense (mem, num) != 0) - error ("IDADense solver not initialized"); + + sunJacMatrix = SUNDenseMatrix (num, num); + if (! sunJacMatrix) + error ("Unable to create dense Jacobian for Sundials"); - if (havejac && IDADlsSetDenseJacFn (mem, IDA::jacdense) != 0) - error ("Dense Jacobian not set"); + sunLinearSolver = SUNDenseLinearSolver (yy, sunJacMatrix); + if (! sunLinearSolver) + error ("Unable to create dense linear solver"); + + if (IDADlsSetLinearSolver (mem, sunLinearSolver, sunJacMatrix)) + error ("Unable to set dense linear solver"); + + if (havejac && IDADlsSetJacFn (mem, IDA::jacdense) != 0) + error ("Unable to set dense Jacobian function"); + } } void - IDA::jacdense_impl (long int Neq, realtype t, realtype cj, - N_Vector& yy, N_Vector& yyp, DlsMat& JJ) + IDA::jacdense_impl (realtype t, realtype cj, + N_Vector& yy, N_Vector& yyp, SUNMatrix& JJ) { BEGIN_INTERRUPT_WITH_EXCEPTIONS; + long int Neq = NV_LENGTH_S(yy); + ColumnVector y = NVecToCol (yy, Neq); ColumnVector yp = NVecToCol (yyp, Neq); @@ -366,7 +400,7 @@ std::copy (jac.fortran_vec (), jac.fortran_vec () + jac.numel (), - JJ->data); + SUNDenseMatrix_Data(JJ)); END_INTERRUPT_WITH_EXCEPTIONS; } @@ -374,7 +408,7 @@ # if defined (HAVE_SUNDIALS_IDAKLU) void IDA::jacsparse_impl (realtype t, realtype cj, N_Vector& yy, N_Vector& yyp, - SlsMat& Jac) + SUNMatrix& Jac) { BEGIN_INTERRUPT_WITH_EXCEPTIONS; @@ -390,17 +424,18 @@ else jac = (*jacspcell) (spdfdy, spdfdyp, cj); - SparseSetMatToZero (Jac); - int *colptrs = *(Jac->colptrs); - int *rowvals = *(Jac->rowvals); + SUNMatZero_Sparse (Jac); + octave_idx_type *colptrs = SUNSparseMatrix_IndexPointers (Jac); + octave_idx_type *rowvals = SUNSparseMatrix_IndexValues (Jac); for (int i = 0; i < num + 1; i++) colptrs[i] = jac.cidx(i); + double *d = SUNSparseMatrix_Data (Jac); for (int i = 0; i < jac.nnz (); i++) { rowvals[i] = jac.ridx(i); - Jac->data[i] = jac.data(i); + d[i] = jac.data(i); } END_INTERRUPT_WITH_EXCEPTIONS; @@ -567,7 +602,7 @@ //main loop while (((posdirection == 1 && tsol < tend) - || (posdirection == 0 && tsol > tend)) + || (posdirection == 0 && tsol > tend)) && status == 0) { if (IDASolve (mem, tend, &tsol, yy, yyp, IDA_ONE_STEP) != 0) @@ -692,7 +727,7 @@ // Linear interpolation ie(0) = index(0); te(0) = tsol - val (index(0)) * (tsol - told) - / (val (index(0)) - oldval (index(0))); + / (val (index(0)) - oldval (index(0))); ColumnVector ytemp = y - ((tsol - te(0)) * (y - yold) / (tsol - told)); @@ -717,7 +752,7 @@ // Linear interpolation ie(temp+i) = index(i); te(temp+i) = tsol - val(index(i)) * (tsol - told) - / (val(index(i)) - oldval(index(i))); + / (val(index(i)) - oldval(index(i))); ColumnVector ytemp = y - (tsol - te (temp + i)) * (y - yold) / (tsol - told); @@ -1096,7 +1131,7 @@ event_fcn = options.getfield("Events").function_value (); // Set up linear solver - dae.set_up (); + dae.set_up (y0); // Integrate retval = dae.integrate (numt, tspan, y0, yp0, refine,