Mercurial > octave

comparison libinterp/dldfcn/__ode15__.cc @ 30547:b876de975edf stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
__ode15__: Adapt for changed API of SUNDIALS version 6 (bug #61701). * m4/acinclude.m4 (SUNLINSOL_KLU), libinterp/dldfcn/__ode15__.cc: Adapt check and functions for changed API of SUNDIALS version 6.
author Markus Mützel <markus.muetzel@gmx.de>
date Sat, 25 Dec 2021 12:36:11 +0100
parents 7a4f16bff8fd
children 00d82e792b8b 796f54d4ddbf
comparison
equal deleted inserted replaced
30545:c2ee60f0d8bf 30547:b876de975edf
192 ~IDA (void) 192 ~IDA (void)
193 { 193 {
194 IDAFree (&m_mem); 194 IDAFree (&m_mem);
195 SUNLinSolFree (m_sunLinearSolver); 195 SUNLinSolFree (m_sunLinearSolver);
196 SUNMatDestroy (m_sunJacMatrix); 196 SUNMatDestroy (m_sunJacMatrix);
197 # if defined (HAVE_SUNDIALS_SUNCONTEXT)
198 SUNContext_Free (&m_sunContext);
199 # endif
197 } 200 }
198 201
199 IDA& 202 IDA&
200 set_jacobian (const octave_value& jac, DAEJacFuncDense j) 203 set_jacobian (const octave_value& jac, DAEJacFuncDense j)
201 { 204 {
251 254
252 void initialize (void); 255 void initialize (void);
253 256
254 static ColumnVector NVecToCol (N_Vector& v, octave_f77_int_type n); 257 static ColumnVector NVecToCol (N_Vector& v, octave_f77_int_type n);
255 258
259 # if defined (HAVE_SUNDIALS_SUNCONTEXT)
260 N_Vector ColToNVec (const ColumnVector& data, octave_f77_int_type n);
261 # else
256 static N_Vector ColToNVec (const ColumnVector& data, octave_f77_int_type n); 262 static N_Vector ColToNVec (const ColumnVector& data, octave_f77_int_type n);
263 # endif
257 264
258 void 265 void
259 set_up (const ColumnVector& y); 266 set_up (const ColumnVector& y);
260 267
261 void 268 void
364 DAERHSFuncIDA m_fun; 371 DAERHSFuncIDA m_fun;
365 DAEJacFuncDense m_jacfun; 372 DAEJacFuncDense m_jacfun;
366 DAEJacFuncSparse m_jacspfun; 373 DAEJacFuncSparse m_jacspfun;
367 DAEJacCellDense m_jacdcell; 374 DAEJacCellDense m_jacdcell;
368 DAEJacCellSparse m_jacspcell; 375 DAEJacCellSparse m_jacspcell;
376 # if defined (HAVE_SUNDIALS_SUNCONTEXT)
377 SUNContext m_sunContext;
378 # endif
369 SUNMatrix m_sunJacMatrix; 379 SUNMatrix m_sunJacMatrix;
370 SUNLinearSolver m_sunLinearSolver; 380 SUNLinearSolver m_sunLinearSolver;
371 }; 381 };
372 382
373 int 383 int
392 realtype *puntrr = nv_data_s (rr); 402 realtype *puntrr = nv_data_s (rr);
393 403
394 for (octave_idx_type i = 0; i < m_num; i++) 404 for (octave_idx_type i = 0; i < m_num; i++)
395 puntrr[i] = res(i); 405 puntrr[i] = res(i);
396 } 406 }
407
408 # if defined (HAVE_SUNDIALS_SUNCONTEXT)
409 # define OCTAVE_SUNCONTEXT , m_sunContext
410 # else
411 # define OCTAVE_SUNCONTEXT
412 # endif
397 413
398 void 414 void
399 IDA::set_up (const ColumnVector& y) 415 IDA::set_up (const ColumnVector& y)
400 { 416 {
401 N_Vector yy = ColToNVec (y, m_num); 417 N_Vector yy = ColToNVec (y, m_num);
405 # if defined (HAVE_SUNDIALS_SUNLINSOL_KLU) 421 # if defined (HAVE_SUNDIALS_SUNLINSOL_KLU)
406 # if defined (HAVE_SUNSPARSEMATRIX_REALLOCATE) 422 # if defined (HAVE_SUNSPARSEMATRIX_REALLOCATE)
407 // Initially allocate memory for 0 entries. We will reallocate when we 423 // Initially allocate memory for 0 entries. We will reallocate when we
408 // get the Jacobian matrix from the user and know the actual number of 424 // get the Jacobian matrix from the user and know the actual number of
409 // entries. 425 // entries.
410 m_sunJacMatrix = SUNSparseMatrix (m_num, m_num, 0, CSC_MAT); 426 m_sunJacMatrix = SUNSparseMatrix (m_num, m_num, 0, CSC_MAT
427 OCTAVE_SUNCONTEXT);
411 # else 428 # else
412 octave_f77_int_type max_elems; 429 octave_f77_int_type max_elems;
413 if (math::int_multiply_overflow (m_num, m_num, &max_elems)) 430 if (math::int_multiply_overflow (m_num, m_num, &max_elems))
414 error ("Unable to allocate memory for sparse Jacobian"); 431 error ("Unable to allocate memory for sparse Jacobian");
415 432
416 m_sunJacMatrix = SUNSparseMatrix (m_num, m_num, max_elems, CSC_MAT); 433 m_sunJacMatrix = SUNSparseMatrix (m_num, m_num, max_elems, CSC_MAT
434 OCTAVE_SUNCONTEXT);
417 # endif 435 # endif
418 if (! m_sunJacMatrix) 436 if (! m_sunJacMatrix)
419 error ("Unable to create sparse Jacobian for Sundials"); 437 error ("Unable to create sparse Jacobian for Sundials");
420 438
421 m_sunLinearSolver = SUNLinSol_KLU (yy, m_sunJacMatrix); 439 m_sunLinearSolver = SUNLinSol_KLU (yy, m_sunJacMatrix
440 OCTAVE_SUNCONTEXT);
422 if (! m_sunLinearSolver) 441 if (! m_sunLinearSolver)
423 error ("Unable to create KLU sparse solver"); 442 error ("Unable to create KLU sparse solver");
424 443
425 if (IDASetLinearSolver (m_mem, m_sunLinearSolver, m_sunJacMatrix)) 444 if (IDASetLinearSolver (m_mem, m_sunLinearSolver, m_sunJacMatrix))
426 error ("Unable to set sparse linear solver"); 445 error ("Unable to set sparse linear solver");
435 454
436 } 455 }
437 else 456 else
438 { 457 {
439 458
440 m_sunJacMatrix = SUNDenseMatrix (m_num, m_num); 459 m_sunJacMatrix = SUNDenseMatrix (m_num, m_num OCTAVE_SUNCONTEXT);
441 if (! m_sunJacMatrix) 460 if (! m_sunJacMatrix)
442 error ("Unable to create dense Jacobian for Sundials"); 461 error ("Unable to create dense Jacobian for Sundials");
443 462
444 m_sunLinearSolver = SUNLinSol_Dense (yy, m_sunJacMatrix); 463 m_sunLinearSolver = SUNLinSol_Dense (yy, m_sunJacMatrix
464 OCTAVE_SUNCONTEXT);
445 if (! m_sunLinearSolver) 465 if (! m_sunLinearSolver)
446 error ("Unable to create dense linear solver"); 466 error ("Unable to create dense linear solver");
447 467
448 if (IDASetLinearSolver (m_mem, m_sunLinearSolver, m_sunJacMatrix)) 468 if (IDASetLinearSolver (m_mem, m_sunLinearSolver, m_sunJacMatrix))
449 error ("Unable to set dense linear solver"); 469 error ("Unable to set dense linear solver");
538 } 558 }
539 559
540 N_Vector 560 N_Vector
541 IDA::ColToNVec (const ColumnVector& data, octave_f77_int_type n) 561 IDA::ColToNVec (const ColumnVector& data, octave_f77_int_type n)
542 { 562 {
543 N_Vector v = N_VNew_Serial (n); 563 N_Vector v = N_VNew_Serial (n OCTAVE_SUNCONTEXT);
544 564
545 realtype *punt = nv_data_s (v); 565 realtype *punt = nv_data_s (v);
546 566
547 for (octave_f77_int_type i = 0; i < n; i++) 567 for (octave_f77_int_type i = 0; i < n; i++)
548 punt[i] = data(i); 568 punt[i] = data(i);
561 581
562 void 582 void
563 IDA::initialize (void) 583 IDA::initialize (void)
564 { 584 {
565 m_num = to_f77_int (m_y0.numel ()); 585 m_num = to_f77_int (m_y0.numel ());
586 # if defined (HAVE_SUNDIALS_SUNCONTEXT)
587 if (SUNContext_Create (nullptr, &m_sunContext) < 0)
588 error ("__ode15__: unable to create context for SUNDIALS");
589 m_mem = IDACreate (m_sunContext);
590 # else
566 m_mem = IDACreate (); 591 m_mem = IDACreate ();
592 # endif
567 593
568 N_Vector yy = ColToNVec (m_y0, m_num); 594 N_Vector yy = ColToNVec (m_y0, m_num);
569 595
570 N_Vector yyp = ColToNVec (m_yp0, m_num); 596 N_Vector yyp = ColToNVec (m_yp0, m_num);
571 597
724 } 750 }
725 751
726 if (status == 0) 752 if (status == 0)
727 { 753 {
728 // Interpolate in tend 754 // Interpolate in tend
729 N_Vector dky = N_VNew_Serial (m_num); 755 N_Vector dky = N_VNew_Serial (m_num OCTAVE_SUNCONTEXT);
730 756
731 if (IDAGetDky (m_mem, tend, 0, dky) != 0) 757 if (IDAGetDky (m_mem, tend, 0, dky) != 0)
732 error ("IDA failed to interpolate y"); 758 error ("IDA failed to interpolate y");
733 759
734 tout(cont) = tend; 760 tout(cont) = tend;
888 const octave_idx_type num_event_args) 914 const octave_idx_type num_event_args)
889 { 915 {
890 realtype h = 0, tcur = 0; 916 realtype h = 0, tcur = 0;
891 bool status = 0; 917 bool status = 0;
892 918
893 N_Vector dky = N_VNew_Serial (m_num); 919 N_Vector dky = N_VNew_Serial (m_num OCTAVE_SUNCONTEXT);
894 920
895 N_Vector dkyp = N_VNew_Serial (m_num); 921 N_Vector dkyp = N_VNew_Serial (m_num OCTAVE_SUNCONTEXT);
896 922
897 ColumnVector yout (m_num); 923 ColumnVector yout (m_num);
898 ColumnVector ypout (m_num); 924 ColumnVector ypout (m_num);
899 std::string string = ""; 925 std::string string = "";
900 926