#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "Standalone.h"

static int fill_structure;	/* Do we need to fill in the structure of    */
				/* the Jacobian?                             */

static double mat[5][5] = {
  {  2.0000,       1.0000,        1.0000,        1.0000,      -1.0000 },
  {  1.0000,       3.0000,        1.0000,        1.0000,       1.0000 }, 
  {  1.0000,       1.0000,        2.0000,        1.0000,       1.0000 }, 
  {  1.0000,       1.0000,        1.0000,        3.0000,      -1.0000 }, 
  { -1.0000,       1.0000,        1.0000,       -1.0000,       4.0000 } 
};


static double q[5] = { 
  8.000000000000000000,
 -2.000000000000000000,
  4.000000000000000000,
 -2.000000000000000000,
  0.000000000000000000
};


int funcEval(int n, double *z, double *f)
{
  /************************************************************************/
  /* Evaluate the function at z, placing the result in f.                 */
  /************************************************************************/

  int i, j;
  
  for (i=0; i<n; i++) {
    f[i] = q[i];
    for (j=0; j<n; j++) {
      f[i] += mat[i][j]*z[j];
    }
  }
  return 0;
}

int jacEval(int n, int nnz, double *z, int *col_start, int *col_len,
	    int *row, double *data)
{
      int i, j;

      if (fill_structure) {
	/*******************************************************************/
	/* This is the first call to jacEval, fill in the structure.       */
	/*******************************************************************/
	for (j=0; j<n; j++) {
	  col_start[j] = j*n + 1;
	  col_len[j] = n;
	  for (i=0; i<n; i++) {
	    row[i + j*n] = i+1;
	  }
	}
	fill_structure = 0;
      }
      
      /***********************************************************************/
      /* Fill in the data values for the jacobian (may change per iteration) */
      /***********************************************************************/

      for (i=0; i<n; i++) {
	for (j=0; j<n; j++) {
	  data[n*j+i] = mat[i][j];
	}
      }
      
     return 0;
}

int main()
{
     double *lb;		/* Lower bounds on the variables             */
     double *ub;		/* Upper bounds on the variables             */
     double *z;			/* Solution vector                           */
     double *f;			/* Function evaluation                       */
     
     int n;			/* Number of variable (N*N)                  */
     int nnz;			/* Number of nonzeros (5*N*N - 4*N)          */
     int status;		/* Termination status from PATH              */
     int i, j;
     
     Output_SetLog(stdout);
     n = 5;

     nnz = n*n;
     
     /************************************************************************/
     /* Allocate space for the bounds, starting point, and function value.   */
     /************************************************************************/

     lb = (double *)malloc(sizeof(double)*n);
     ub = (double *)malloc(sizeof(double)*n);
     z = (double *)malloc(sizeof(double)*n);
     f = (double *)malloc(sizeof(double)*n);

     /************************************************************************/
     /* Fill in the lower and upper bounds and a starting point.             */
     /************************************************************************/

     for (j = 0; j < n; j++)
       {
	 lb[j] = 0.0;
	 ub[j] = 1e20;
	 z[j] = 0.1;
       }

     /************************************************************************/
     /* Set fill_structure to true (we ALWAYS need to fill in the structure  */
     /* of the jacobian at least once per call to PATH).                     */
     /************************************************************************/

     fill_structure = 1;

     /************************************************************************/
     /* Call PATH.                                                           */
     /************************************************************************/

     pathMain(n, nnz, &status, z, f, lb, ub);

     /************************************************************************/
     /* Deallocate memory.                                                   */
     /************************************************************************/

     free(lb);
     free(ub);
     free(z);
     free(f);

     return 0;
}