# what is the structure of fjac vector for dtrnlsp_solve?

## what is the structure of fjac vector for dtrnlsp_solve?

Hello. I am trying to use dtrnlsp_solve function, but I do not know the format of fjac. I calculate fjac myself. According to the manual, it is a 1D array of size m*n where m is the number of function values and n is the number of variables. But, the manual does not say how the fjack elements are arranged. There are two possible ways:

for(i=0;i
or

for(j=0;j
Which one is correct?

13 posts / 0 new
For more complete information about compiler optimizations, see our Optimization Notice.

Quoting - gpwr9k5
Hello. I am trying to use dtrnlsp_solve function, but I do not know the format of fjac. I calculate fjac myself. According to the manual, it is a 1D array of size m*n where m is the number of function values and n is the number of variables. But, the manual does not say how the fjack elements are arranged. There are two possible ways:

for(i=0;i
or

for(j=0;j
Which one is correct?

Hi,

The Fortran version is an array FJac(m,n), so I believe it goes the same as with other two-dimensional matrix A of size m x n stored densely in a one-dimensional array:
A[i][j] = B[i*n+j] in C (i=0, ... , m-1, j=0, ... , n-1)
A(i,j) = B(j*m+i) in Fortran (i=1, ... , m, j=1, ... , n).

so this would be your first version. See Chapter 7 of User's Manual.

A.

Quoting - ArturGuzik

Hi,

The Fortran version is an array FJac(m,n), so I believe it goes the same as with other two-dimensional matrix A of size m x n stored densely in a one-dimensional array:
A[i][j] = B[i*n+j] in C (i=0, ... , m-1, j=0, ... , n-1)
A(i,j) = B(j*m+i) in Fortran (i=1, ... , m, j=1, ... , n).

so this would be your first version. See Chapter 7 of User's Manual.

A.

Thanks. That is what I have been using. Something does not work. The solver keeps exiting with 0 iterations and st_cr=6. I tried different values of rs and even changed the sign of jacobian, but still the same result. If I set eps[5]=0 to avoid early exit, then initialization fails. If you can advice on how to avoid early exit with st_cr=6, that would be great.
Thanks again.

Quoting - gpwr9k5

Thanks. That is what I have been using. Something does not work. The solver keeps exiting with 0 iterations and st_cr=6. I tried different values of rs and even changed the sign of jacobian, but still the same result. If I set eps[5]=0 to avoid early exit, then initialization fails. If you can advice on how to avoid early exit with st_cr=6, that would be great.
Thanks again.

Well, sounds as problems with initialization/settings.

(1) does MKL example (from disk) work for you?

(2) show a bit more of _init call, so we can take a look.

str_cr = 6 means that solution converged, so maybe you have a FJac zero(s) or ..... Do you call it properly (handle, pointers etc)?

A.

Quoting - ArturGuzik
Well, sounds as problems with initialization/settings.

(1) does MKL example (from disk) work for you?

(2) show a bit more of _init call, so we can take a look.

str_cr = 6 means that solution converged, so maybe you have a FJac zero(s) or ..... Do you call it properly (handle, pointers etc)?

A.

The example from disk works, saying "PASS". I am using this solver to train a neural network. I just copied the example code and changed everything that was related to my task. From what I see, all handles and variables get set correctly. Jacobian is not zero. dtrnlsp_solve computes the starting and ending residuals and they are the same. I compute the same rezidual myself (as mySSE in the code below) and it is way different.

```double NN::train(double **in,double *tgt,int nep,const double maxSSE)
{
int nit1=1000;			// maximum number of iterations
int nit2=100;			// maximum number of iterations of calculation of trial-step
double rs=1.0;			// initial step bound
double eps[6];			// precisions for stop-criteria
int RCI_Request;		// reverse communication interface variable
int successful;			// rci cycle variable
int nit;				// number of iterations performed
int st_cr;				// number of stop-criterion
double init_sse,sse;	// initial and final sse's
_TRNSP_HANDLE_t handle;	// TR solver handle

//	set precisions for stop-criteria
for(int i=0;i<6;i++) eps[i]=1.0e-30;
cout << np << endl;

//	initialize solver (allocate memory, set initial values)
if(dtrnlsp_init(&handle,&nw,&np,wv,eps,&nit1,&nit2,&rs) != TR_SUCCESS)
{
cout << "error in dtrnlsp_initn";
cin.get();
return 0;
}

//	set initial rci cycle variables
RCI_Request = 0;
successful = 0;
double mySSE=0.0;

//	main training cycle
while(successful == 0)
{
//	call tr solver, which reurns RCI_request telling us the next step
if (dtrnlsp_solve(&handle,err,jac,&RCI_Request) != TR_SUCCESS)
{
cout << "error in dtrnlsp_solven";
cin.get();
return 0;
}

//	act on RCI_request
if (RCI_Request == -1 || RCI_Request == -2 || RCI_Request == -3 ||
RCI_Request == -4 || RCI_Request == -5 || RCI_Request == -6)
successful = 1; // one of the eps is below maxSSE; exit rci cycle
if (RCI_Request == 1)
{
//	recalculate errors
for(int p=0;p ```

Quoting - gpwr9k5

```double NN::train(double **in,double *tgt,int nep,const double maxSSE)
{

//	initialize solver (allocate memory, set initial values)
if(dtrnlsp_init(&handle,&nw,&np,wv,eps,&nit1,&nit2,&rs) != TR_SUCCESS)

```

Quick look -- nothing suspicious. Where do you get/set nw, np, that is sizes of the problem/arrays (used in _init).

A.

Hi all,

I would like to recommend to change the precisions: for(int i=0;i<6;i++) eps[i]=1.0e-30;
I suggest to use 1.0e-12 or large.

Thank you!
--Nikita

Quoting - ArturGuzik

Quick look -- nothing suspicious. Where do you get/set nw, np, that is sizes of the problem/arrays (used in _init).

A.

I passthem when I call my train() function. nw=71 (number of weoghts/variables), np=1500 (number of patterns/functions). So, my fjac array is huge, 71x1500=106,500. Could that be a problem? I tried to increase esp[] to 1e-12, but now it quits with st_cr=2. My biggest puzzle is why the initial and final residuals computed by the solver are the same and so different from mine. I pass the solver thevalues err[p]=tgt[p]-out[p] as the function values. I compute an SSE as sum(err[p]^2,p=0..np-1). The solver probably uses something else. For instance, solver computes 42.9975. I compute 77431.3. Even if I divide my sse by (np) or (2*np), I still get a different residual that the solver's. How does the solver computes the residual from the passed function values?

Quoting - gpwr9k5

I passthem when I call my train() function. nw=71 (number of weoghts/variables), np=1500 (number of patterns/functions). So, my fjac array is huge, 71x1500=106,500. Could that be a problem? I tried to increase esp[] to 1e-12, but now it quits with st_cr=2. My biggest puzzle is why the initial and final residuals computed by the solver are the same and so different from mine. I pass the solver thevalues err[p]=tgt[p]-out[p] as the function values. I compute an SSE as sum(err[p]^2,p=0..np-1). The solver probably uses something else. For instance, solver computes 42.9975. I compute 77431.3. Even if I divide my sse by (np) or (2*np), I still get a different residual that the solver's. How does the solver computes the residual from the passed function values?

Quoting - gpwr9k5

Hi,

A.

Quoting - gpwr9k5

And one more, your Jac is 1500 x 71, not 71 x 1500, right?

A.

Quoting - ArturGuzik
And one more, your Jac is 1500 x 71, not 71 x 1500, right?

A.

Number of function values m=1500
Number of variables n=71.
My jac is a 1D array oflength 21,300 computed as jac[i*n+j]=dFunction[j]/dVariable[j], i=0..m-1, j=0..n-1
I got the residuals to agree with my calculations: r1/2=sqrt(Function[j]^2,j=0..m-1).
I changed the m from 1500 to 300 and still got stopped with st_cr=6. All eps[]=1e-30. I checked my jac values and their magnituderanges from 1e-6 to 1.
Here is the whole code of the function, that is doing all calculations (nw=n - number of variables, np=m - number of function values, err[p], p=0,..,np-1 - evaluated function, which = tgt[p]-Function_Output)

```double NN::train(double **in,double *tgt,int nep,const double maxSSE)
{
int nit1=1000;			// maximum number of iterations
int nit2=100;			// maximum number of iterations of calculation of trial-step
double rs=1.0;			// initial step bound
double eps[6];			// precisions for stop-criteria
int RCI_Request;		// reverse communication interface variable
int successful;			// rci cycle variable
int nit;				// number of iterations performed
int st_cr;				// number of stop-criterion
double init_sse,sse;	// initial and final sse's
_TRNSP_HANDLE_t handle;	// TR solver handle

//	set precisions for stop-criteria
for(int i=0;i<6;i++) eps[i]=1.0e-50;

//	initialize solver (allocate memory, set initial values)
if(dtrnlsp_init(&handle,&nw,&np,wv,eps,&nit1,&nit2,&rs) != TR_SUCCESS)
{
cout << "error in dtrnlsp_initn";
cin.get();
return 0;
}

//	set initial rci cycle variables
RCI_Request = 0;
successful = 0;
double mySSE;

//	main training cycle
while(successful == 0)
{
//	call tr solver, which reurns RCI_request telling us the next step
if (dtrnlsp_solve(&handle,err,jac,&RCI_Request) != TR_SUCCESS)
{
cout << "error in dtrnlsp_solven";
cin.get();
return 0;
}

//	act on RCI_request
if (RCI_Request == -1 || RCI_Request == -2 || RCI_Request == -3 ||
RCI_Request == -4 || RCI_Request == -5 || RCI_Request == -6)
successful = 1; // one of the eps is below maxSSE; exit rci cycle
if (RCI_Request == 1)
{
//	recalculate errors
mySSE=0.0;
for(int p=0;p < np;p++)					// for each training pattern
{
ffwd(in[p]);
err[p]=tgt[p]-out[nl-1][0];	// only one neuron in the output layer
mySSE+=err[p]*err[p];
}
}
if (RCI_Request == 2)
{
//	compute Jacobi matrix
cout << "compute Jacobiann";
for(int p=0;p < np;p++)					// for each training pattern
{
//	update outputs of each neuron
ffwd(in[p]);

//	find sensitivity for the output layer i=nl-1
s[nl-1][0]=-1.0;
if(oaf==1)
s[nl-1][0]*=afDeriv(out[nl-1][0]);

//	propagate sensitivities from output layer to hidden layers
for(int i=nl-2;i > 0;i--)				// for each layer except input & output
for(int j=0;j < ls[i];j++)		// for each neuron in current layer
{
double sum=0.0;
for(int k=0;k < ls[i+1];k++)	// for each neuron in later layer
sum+=s[i+1][k]*w[i+1][k][j];
s[i][j]=afDeriv(out[i][j])*sum;
}

int iw=0;
for(int i=1;i < nl;i++)					// for each layer except input
for(int j=0;j < ls[i];j++)			// for each neuron in current layer
for(int k=0;k <= ls[i-1];k++)		// for each input of curr neuron incl bias
{
if(k==ls[i-1]) jac[p*nw+iw]=s[i][j];		// bias input
else jac[p*nw+iw]=s[i][j]*out[i-1][k];	// data input
if(p*nw+iw>=21100)
cout << i << "," << j << "," << k << ": " << jac[p*nw+iw] << endl;
iw++;
}
}
}
}

// get solution statuse
if (dtrnlsp_get(&handle,&nit,&st_cr,&init_sse,&sse) != TR_SUCCESS)
{
cout << "error in dtrnlsp_getn";
cin.get();
return 0;
}

//	free handle memory
if (dtrnlsp_delete(&handle) != TR_SUCCESS)
{
cout << "error in dtrnlsp_deleten";
cin.get();
return 0;
}
cout << "finished after " << nit << " iterationsn";
cout << "stop criterion is " << st_cr << endl;
cout << "sse: " << init_sse << " -> " << sse << endl;
cout << "mySSE: " << sqrt(mySSE) << endl;
return(sse);
}```

I think I found my mistake