I noticed that the nonlinear solver uses central difference to calculate the Jacobi matrix, I was wondering is it possible to pass the analytical Jacobi elements to the code or not? I mean, in my problem, I want to use analytical Jacobi instead of numerical Jacobi. How can I do that?
The nonlinear solver uses the reverse call interface. When RCI_REQUEST is equal to 2, you have to update the Jacobian matrix FJAC and call the solver again. The method used for calculating the Jacobian is entirely in your hands. You can enter formulae for the elements of the Jacobian in line, or in a separate subroutine that you write and call; you can use finite difference approximations for some or all elements, or obtain the derivatives from a finite-field extension using complex variables.
You can, for example, take the furnished testcase ex_nlsqp_f.f and modify it so that FJAC is calculated analytically instead of calling DJACOBI.
Thanks for your reply,
I am using this solver in C++, I used the "ex_nlsqp_bc_c.c" example and tried to implement what you mentioned as the following.
I defined the function as
extern void jacobian(MKL_INT *, MKL_INT *, double *, double *);
I called the Jacobian function as
if (RCI_Request == 2)
jacobian(&m, &n, x, fjac);
and defined the Jacobian like the following:
void jacobian(MKL_INT * m, MKL_INT * n, double *x, double *fj)
fj = 2 * x;
fj = 2 * x;
fj = 1;
fj = -1;
The function was defined as:
void extended_powell(MKL_INT * m, MKL_INT * n, double *x, double *f)
f = pow(x,2) + pow(x,2) - 9;
f = x - x + 3;
When I call the jacobian, unfortunately it does not work, but if I call djacobi it works.
I tracked the iteration of jacobian, it does around 20 iterations and becomes close to the solution, but I do not know why it fails after several iteration while it is getting closer to the correct solution.
The correct solution for this problem is "zero" and "3", and in the following I printed the two last iterations of the solution when I am using jacobian.
x1= 0.10554969128064530 x2= 2.9935339376683032
x1= 0.10554965354656229 x2= 2.9935339407828785
Can anyone help me?
I wrote a detailed reply, but the forum software marked it as spam and removed it. Here is a short replacement:
Your Jacobian vector arranges the elements of the 2-dimensional matrix in row-major order. The MKL routine may expect it to be in column-major (as in Fortran) order. Try printing out the values of fjac from the finite-difference calculation and compare those with the analytically calculated values for accuracy and matching order.