LLNL / sundials

Official development repository for SUNDIALS - a SUite of Nonlinear and DIfferential/ALgebraic equation Solvers. Pull requests are welcome for bug fixes and minor changes.
https://computing.llnl.gov/projects/sundials
BSD 3-Clause "New" or "Revised" License
507 stars 123 forks source link

onemkl linear solver doesn't have a default difference quotient jacobian? #182

Open ciropom opened 2 years ago

ciropom commented 2 years ago

Hello again, I finally was able to compile link and execute my model with the onemkl dense linear solver. The problem now is that at initialization time I get this error

[CVLS ERROR]  cvLsInitialize
  No Jacobian constructor available for SUNMatrix type

[CVODE ERROR]  cvInitialSetup
  The linear solver's init routine failed.

I'm using the following initialization code:

    /* Create dense SUNMatrix for use in linear solves */
    A = SUNMatrix_OneMklDense(y_len, y_len, SUNMEMTYPE_DEVICE, memhelper, &s->qspccQueue, s->sunctx);
    if (check_flag((void *)A, "SUNMatrix_OneMklDense", 0)) return s;
    s->matrix = A; /* save for deallocation */
    // Create the SUNLinearSolver object for use by CVode
    LS = SUNLinSol_OneMklDense(s->y, A, s->sunctx);
    if (check_flag((void *)LS, "SUNLinSol_OneMklDense", 0)) return s;
    s->flag = CVodeSetLinearSolver(s->cvode_mem, LS, A);
    if(check_flag(&s->flag, "CVodeSetLinearSolver", 1)) return s;    

It looks like onemkldense solver doesn't have a default difference-quotient jacobian, like the standard dense/band/lapack solvers. Is this true? Users are forced to provide a Jacobian through CVodeSetJacFn to use this solver?

As always thank you very much for your help.

gardner48 commented 2 years ago

Yes, that's correct. A user-supplied Jacobian function is currently required when using the oneMKL interface. This function may compute the analytic Jacobian or an approximation to it e.g., using difference-quotients.

ciropom commented 2 years ago

Why there is not any pre-packaged difference-quotient jacobian as it is for the standard dense solver?

balos1 commented 2 years ago

Why there is not any pre-packaged difference-quotient jacobian as it is for the standard dense solver?

The implementation of the difference-quotients routine in CVODE currently is specific to the standard dense and banded matrices. Generalizing it to support other matrices is something that I think would be nice, but we don't currently have any plans to do so.

ciropom commented 2 years ago

Thank you for the explanation. This is blocking in my scenario, where I'm translating automatically from MATLAB to C and simulating with sundials by using QSPcc

Can you give me some hints on how to implement a difference-quotient jacobian for the MKL matrix? Thanks

drreynolds commented 2 years ago

@ciropom: you will need to emulate the code from cvLsDenseDQJac (here), but ensuring that each of the operations use device-resident vector data, and that these operations fill device-resident data within the MKL matrix.