topic Help for Jacobian calculations in Intel® oneAPI Math Kernel Library

Help for Jacobian calculations

sureshdeoda — Thu, 22 Sep 2011 08:31:22 GMT

Hi,

I was trying to use trust region solver for one of my code. I found that the jacobian given by the solver is calculating only half of thevalues (i.e. n*m/2)whilefor rest of thecalculation it gives zeros and brakes out with the RCI_Request == -5. ??
I have tested solver with some small eqation where i have calculated the jacobian(i.e. (f(x) - f(x'))/(x - x')) and compared the same with jacobian given by the solver which are same but not for my original code.
One more thing is i found that, though i am giving initial guess as upper bounds ,while perturbing the variable in costfun calculation goes beyond upper bounds .??

I know without code it will be difficul to find the errors, hope someone might be able to hint the possible reasons.

followingare thecalls which i use. it may help

if (dtrnlspbc_init (&handle, &n, &m, x, LW, UP, eps, &iter1, &iter2,&rs) != TR_SUCCESS) {}
if (dtrnlspbc_solve (&handle, fvec, fjac, &RCI_Request) != TR_SUCCESS) {}
if (RCI_Request == 1)
{
CostFun(&m,&n,x,fvec);
}
if (RCI_Request == 2)
{
if (djacobi (CostFun, &n, &m, fjac, x, eps) != TR_SUCCESS)
{
return 0;
}
fp_temp= fopen("solverResFjac.dat","w");
for(i=0;i fprintf(fp_temp,"%f\\n",fjac);

Jacobian(fjac1,n,m,x,fvec); // user defined jacobian
fp_temp= fopen("solverResFjac1.dat","w");
for(i=0;i fprintf(fp_temp,"%f\\n",fjac1);
}

Thanks ,
Suresh

Help for Jacobian calculations

Gennady_F_Intel — Thu, 22 Sep 2011 10:26:58 GMT

>> the RCI_Request == -5. ??

This is an unexpected output task status - please show a complete example demonstrating the claimed bug.

Help for Jacobian calculations

sureshdeoda — Fri, 30 Sep 2011 06:08:18 GMT

Hi Gennady ,

Thanks for the reply. I found there are someproblem with my code..:) .

But still if my initial guess is one of the boundary condition, in jacobian it crosses the boundary condition which i dont want to happen ,So i am calculating jacobian.

I found the values are same in solver calculated and user defined but still the next iteartion values are different with both the cases.

Please let me know if i am missing something here..

[bash]#include 
#include 
#include 
#include 

void myJacobi(double *c,double *y,double *fjec,int n,int m,double *fvec);
void myCostFun(int *m, int *n, double *x, double *f);
double y[7];

int main()
{
			
	
	//double y[7];
	double xin[7];
	double c[3] = {3 ,2 ,3};
			
	/* n - number of function variables
	m - dimension of function value */
	int n = 3, m = 7;
	/* precisions for stop-criteria (see manual for more details) */
	double eps[6];
	/* solution vector. contains values x for f(x) */
	double *x;
	/* iter1 - maximum number of iterations
	iter2 - maximum number of iterations of calculation of trial-step */
	int iter1 = 1000, iter2 = 100;
	/* initial step bound */
	double rs = 100.0;
	/* reverse communication interface parameter */
	int RCI_Request;

	/* controls of rci cycle */
	int successful;
	/* function (f(x)) value vector */
	double *fvec;
	/* jacobi matrix */
	double *fjac,*fjac1;
	/* lower and upper bounds */
	double *LW, *UP;
	/* number of iterations */
	int iter;
	/* number of stop-criterion */
	int st_cr;
	/* initial and final residuals */
	double r1, r2;
	/* TR solver handle */
	_TRNSPBC_HANDLE_t handle;
	/* cycles counter */
	int i;
	
	/* memory allocation */
	x = (double*) malloc (sizeof (double)*n);
	fvec = (double*) malloc (sizeof (double)*m);
	LW = (double*) malloc (sizeof (double)*n);
	UP = (double*) malloc (sizeof (double)*n);
	fjac = (double*) malloc (sizeof (double)*m*n);
	fjac1 = (double*) calloc (m*n,sizeof (double));
	/* set precisions for stop-criteria */
	for (i = 0; i < 6; i++)
	{
		eps  = 0.00001;
	}

	/* set the initial guess */
	for (i = 0; i < n; i++)
	{
	x  = 4.2;	
	}
	/* set the initial values */
	for (i = 0; i < m; i++)
	fvec  = 0.0;
	for (i = 0; i < m*n; i++)
	fjac  = 0.0;
	/* set bounds */
	for (i = 0; i < n; i++)
	{
	LW  = 0;	
	UP  = 4;	
	}

	for(i=0;i<7;i++)
	{
		xin = i+1;		
		y = c[0]*xin +c[1]*xin*xin +c[2]*xin*xin;
	}

	if (dtrnlspbc_init (&handle, &n, &m, x, LW, UP, eps, &iter1, &iter2,
	&rs) != TR_SUCCESS)
	{	
		printf ("| error in dtrnlspbc_initn");	
		return 0;
	}

	/* set initial rci cycle variables */
	RCI_Request = 0;
	successful = 0;
	/* rci cycle */
	while (successful == 0)
	{
	if (dtrnlspbc_solve (&handle, fvec, fjac, &RCI_Request) != TR_SUCCESS)
	{
		printf ("error in dtrnlspbc_solven");
		return 0;
	}
	/* according to rci_request value we do next step */
	if (RCI_Request == -1 ||
	RCI_Request == -2 ||
	RCI_Request == -3 ||
	RCI_Request == -4 ||
	RCI_Request == -5 ||
	RCI_Request == -6)
	/* exit rci cycle */
	successful = 1;

	if (RCI_Request == 1)
	{			
		myCostFun(&m, &n, x, fvec);	
		printf("Cost function n");
		for (i=0;i);
		}
	}
	if (RCI_Request == 2)
	{
		/* compute jacobi matrix*/	
		if (djacobi(myCostFun, &n, &m, fjac, x, eps) != TR_SUCCESS)
		{			
			printf ("error in djacobin");		
			return 0;
		}
		
		myJacobi(x,y,fjac1,n,m,fvec);
		printf("Jacobiann MKLtt User definedn");
		for (i=0;i,fjac1);
		}
	}
	}
	
	if (dtrnlspbc_get (&handle, &iter, &st_cr, &r1, &r2) != TR_SUCCESS)
	{			
		return 0;
	}
	
	
	if (dtrnlspbc_delete (&handle) != TR_SUCCESS)
	{
		return 0;
	}


	printf ("The result isn");	
	for(i=0;i);	
	//printf ("n The Cost is %fn");	

	/* free allocated memory */
	free (x);
	free (fvec);
	free (fjac);
	free (LW);
	free (UP);
	/* if final residual is less than required precision then print pass */
	if (r2 < 0.1)
	printf ("| dtrnlspbc ..........PASSn");
	/* else print failed */
	else
	printf ("| dtrnlspbc .........FAILEDn");	
	
	
	getch();

	return 0;
	
}



void myCostFun(int *m, int *n, double *x, double *f) //  (double* c,double *y,double *fvec)
{	
	int  i= 0;
	double fx[7];
	double xin[7];

	for(i=0;i<7;i++)
	{
		xin = i+1;
		fx = x[0]*xin +x[1]*xin*xin +x[2]*xin*xin;		
		f = (y-fx);
	}	
	//return resCost;
}



void myJacobi(double *c1,double *y,double *fjec,int n,int m,double *fvec)
{
	int i =0,j=0,k=0;
	float x[7];
	double fx[7];
	double fv[7];
	float var = 0.01;
	float del;
	double c[3];
	
	for(j = 0;j = c1;
		}
		del = (c+var*c)- c;
		c= c+var*c;
		for(i=0;i = i+1;
						
			fx = c[0]*x +c[1]*x*x +c[2]*x*x;			
			fv = (y-fx);

			fjec[j*m+i] = (fv - fvec)/del;

		}
	}


}[/bash]

Thanks,
Suresh

Help for Jacobian calculations

Gennady_F_Intel — Fri, 30 Sep 2011 10:42:20 GMT

well, sounds better :), but somehow I don't see when thenext iteration values are different with both the cases.

please show these values which you received.

Help for Jacobian calculations

sureshdeoda — Fri, 30 Sep 2011 11:22:56 GMT

Hi,

The output is as follows..

Cost function in 1st iteration -

-4.600000

-16.000000

-34.200000

-59.200000

-91.000000

-129.600000

-175.000000

Jacobian -

MKL User Defined

-1.000000 -1.000000

-2.000000 -2.000000

-3.000000 -3.000000

-4.000000 -4.000000

-5.000000 -5.000000

-6.000000 -6.000000

-7.000000 -7.000000

-1.000000 -1.000000

-4.000000 -4.000000

-9.000000 -9.000000

-16.000000 -16.000000

-25.000000 -25.000000

-36.000000 -36.000000

-49.000000 -49.000000

-1.000000 -1.000000

-4.000000 -4.000000

-9.000000 -9.000000

-16.000000 -16.000000

-25.000000 -25.000000

-36.000000 -36.000000

-49.000000 -49.000000

These are the cost function and jacobian for 1st iteration which are same..

But the cost function( new values of variable ) with djacobi is ...
Cost function in 2nd iteration

-0.022444

-0.036071

-0.040883

-0.036880

-0.024060

-0.002425

0.028026

whereas afterusing myjacobi the output is,

Cost function in 2nd iteartion

0.780000

3.120000

7.020000

12.480000

19.500000

28.080000

38.220000

if i feed the same data as cost and jacobian why there should be difference in next iteration?

Thanks,
Suresh