#include "snopt.h"

typedef cint (*snFOBJ) /* Pointer to a function returning an cint */
(
 cint*   modefg,
 cint*   ierror,
 cint*   state,
 cint*   nnObj,
 void*   fgobj,
 double* fobj,
 double* gobj,
 double* x,
 char*   cu,
 cint*   lencu,
 cint*   iu,
 cint*   leniu,
 double* ru,
 cint*   lenru,
 char*   cw,
 cint*   lencw,
 cint*   iw,
 cint*   leniw,
 double* rw,
 cint*   lenrw,
 long    cu_len,
 long    cw_len
);

typedef cint (*snFCON)
(
 cint*   modefg,
 cint*   ierror,
 cint*   state,
 cint*   n,
 cint*   nnCon,
 cint*   nnJac,
 cint*   negCon,
 void*   fgcon,
 cint*   ne,
 cint*   nlocJ,
 cint*   indJ,
 cint*   locJ,
 double* fcon,
 double* gcon,
 double* x,
 char*   cu,
 cint*   lencu,
 cint*   iu,
 cint*   leniu,
 double* ru,
 cint*   lenru,
 char*   cw,
 cint*   lencw,
 cint*   iw,
 cint*   leniw,
 double* rw,
 cint*   lenrw,
 long    cu_len,
 long    cw_len
);

extern cint snmemb_
(
 cint   *iExit,
 cint   *m,
 cint   *n,
 cint   *ne,
 cint   *negCon,
 cint   *nnCon,
 cint   *nnJac,
 cint   *nnObj,
 cint   *mincw,
 cint   *miniw,
 cint   *minrw,
 char   *cw,
 cint   *lencw,
 cint   *iw,
 cint   *leniw,
 double *rw,
 cint   *lenrw,
 long    cw_len
 );

extern cint snkerb_
(
 char   *start,
 cint   *m,
 cint   *n,
 cint   *ne,
 cint   *nName,
 cint   *nnCon,
 cint   *nnObj,
 cint   *nnJac,
 cint   *iObj,
 double *ObjAdd,
 char   *prob,
 snConstraintFunc fgcon,   /* function pointers */
 snObjectiveFunc  fgobj,
 cint   (*snLog )(),
 cint   (*snLog2)(),
 cint   (*sqLog )(),
 cint   (*snSTOP)(),
 double *Jcol,
 cint   *indJ,
 cint   *locJ,
 double *bl,
 double *bu,
 char   *Names,
 cint   *hs,
 double *x,
 double *pi,
 double *rc,
 cint   *INFO,
 cint   *mincw,
 cint   *miniw,
 cint   *minrw,
 cint   *nS,
 cint   *nInf,
 double *sInf,
 double *Obj,
 char   *cu,
 cint   *lencu,
 cint   *iu,
 cint    *leniu,
 double *ru,
 cint   *lenru,
 char   *cw,
 cint   *lencw,
 cint   *iw,
 cint   *leniw,
 double *rw,
 cint   *lenrw,
 long    start_len,
 long    prob_len,
 long    Names_len,
 long    cu_len,
 long    cw_len
 );

extern cint sninit_
(
 cint   *iprint,
 cint   *isumm,
 char   *cw,
 cint   *lencw,
 cint   *iw,
 cint   *leniw,
 double *rw,
 cint   *lenrw,
 long    cw_len
);

extern cint snopen_
(
 char   *fname,
 cint   *unit,
 cint   *status,
 long    fname_clen
);

void snOpenFiles( snProblem* prob )
{
  cint       str_len;
  cint       inform;
  const char std_out[] = "stdout";

  /* open no file, standard out, or a disk file */
  if ( prob->option.print_file_name == NULL )
    prob->print_file_unit = 0;
  else if ( strcmp( prob->option.print_file_name, std_out ) == 0 )
    prob->print_file_unit = 6;
  else {
    prob->print_file_unit = 9;
    str_len = strlen( prob->option.print_file_name );
    snopen_(   prob->option.print_file_name,
	     &(prob->print_file_unit) ,
	     &inform , str_len );
  }

  if ( prob->option.sumry_file_name == NULL )
    prob->sumry_file_unit = 0;
  else if ( strcmp( prob->option.sumry_file_name , std_out ) == 0 )
    prob->sumry_file_unit = 6;
  else {
    prob->sumry_file_unit = 4;
    str_len = strlen( prob->option.sumry_file_name );
    snopen_( prob->option.sumry_file_name,
	     &(prob->sumry_file_unit) ,
	     &inform , str_len );
  }
}

/* *************************************************************** */

cint snSolve( snProblem* prob )
{
  char   start[5];

  cint   INFO;

  cint   nName = 1;
  cint   mincw, miniw, minrw;
  char   ProbName[8];
  char   Names[8*1]; /* nName is one; */

  /*  Copy the problem name into the work array.       */
  if (prob->option.ProblemName == NULL )
    strncpy( ProbName,              "    None" , 8 );
  else {
    strncpy( ProbName, prob->option.ProblemName, 8 );
  }

  /* Set the default start parameter.                  */
  /* It may be superceded by a value set as an option. */

  strcpy( start , "Cold" );


  /*  Solve the problem. */

  snkerb_(start,
	  &(prob->dimension.m),
	  &(prob->dimension.n),
	  &(prob->dimension.ne),
	  &nName,
	  &(prob->dimension.nnCon),
	  &(prob->dimension.nnObj),
	  &(prob->dimension.nnJac),
	  &(prob->linear_objective_row),
	  &(prob->option.objective_add),
	  ProbName,
	  prob->constraint_func,
	  prob->objective_func,
	  snlog_,
	  snlog2_,
	  sqlog_,
	  snstop_,
	  prob->Jcol,
	  prob->indJ,
	  prob->locJ,
	  prob->bl,
	  prob->bu,
	  Names,
	  prob->hs,
	  prob->xs,
	  prob->pi,
	  prob->rc,
	  &INFO,
	  &mincw,
	  &miniw,
	  &minrw,
	  &(prob->stats.num_superbasics),
	  &(prob->stats.num_infeasibilities),
	  &(prob->stats.sum_infeasibilities),
	  &(prob->stats.objective_value),
	  /* Pass problem workspace as user workspace */
	  prob->cw,
	  &(prob->length_cw),
	  prob->iw,
	  &(prob->length_iw),
	  prob->rw,
	  &(prob->length_rw) ,
	  prob->cw,
	  &(prob->length_cw),
          prob->iw,
          &(prob->length_iw),
          prob->rw,
          &(prob->length_rw) ,
          sizeof(start),            /* 4L start_len */
          sizeof(ProbName),         /* 8L prob_len  */
          sizeof(Names),            /* 8L Names_len */
          prob->length_cw,          /* size of character input   */
	  prob->length_cw );

  prob->stats.num_minor_iterations              = prob->iw[420];
  prob->stats.num_major_iterations              = prob->iw[421];
  prob->stats.num_objective_function_evals      = prob->iw[403];

  prob->stats.num_constraint_function_evals     = prob->iw[404];
  prob->stats.num_degenerate_steps              = prob->iw[413];
  prob->stats.i_max_scaled_primal_infeasibility = prob->iw[414];
  prob->stats.i_max_scaled_dual_infeasibility   = prob->iw[415];
  prob->stats.i_max_primal_infeasibility        = prob->iw[416];
  prob->stats.i_max_dual_infeasibility          = prob->iw[417];

  prob->stats.elastic_weight                  = prob->rw[423];
  prob->stats.scaled_merit                    = prob->rw[403];
  prob->stats.linear_objective_value          = prob->rw[433];
  prob->stats.penalty                         = prob->rw[435];
  prob->stats.norm_scaled_solution            = prob->rw[422];
  prob->stats.norm_scaled_multipliers         = prob->rw[409];
  prob->stats.norm_solution                   = prob->rw[410];
  prob->stats.norm_multipliers                = prob->rw[411];
  prob->stats.max_scaled_primal_infeasibility = prob->rw[412];
  prob->stats.max_scaled_dual_infeasibility   = prob->rw[413];
  prob->stats.max_primal_infeasibility        = prob->rw[414];
  prob->stats.max_dual_infeasibility          = prob->rw[415];
  prob->stats.max_nonlinear_violation         = prob->rw[416];

  prob->set_flags |= SN_BIT_SOLVED;  /* snSolve() has been called */

  switch (INFO) {
  case 0:
    prob->solve_status = SN_SOLUTION_FOUND;
    return SN_OK;
  case 1:
    prob->solve_status = SN_INFEASIBLE;
    return SN_OK;
  case 2:
    prob->solve_status = SN_UNBOUNDED;
    return SN_OK;
    /* or  SN_VIOLATION_LIMIT_EXCEEDED; */
  case 3:
    prob->solve_status = SN_MINOR_ITERATION_LIMIT_EXCEEDED;
    return SN_OK;
    /* or  SN_MAJOR_ITERATION_LIMIT_EXCEEDED; */
  case 4:
    prob->solve_status = SN_ACCURACY_NOT_ACHIEVED;
    return SN_OK;
  case 5:
    prob->solve_status = SN_SUPERBASICS_LIMIT_EXCEEDED;
    return SN_OK;
  case 9:
    prob->solve_status = SN_POINT_CANNOT_BE_IMPROVED;
    return SN_OK;
  case 10:
    prob->solve_status = SN_CANNOT_SATISFY_GENERAL_CONSTRAINTS;
    return SN_OK;
  case 22:
    prob->solve_status = SN_SINGULAR_BASIS;
    return SN_OK;
  }

  return SN_ERROR;
}

/* *************************************************************** */

cint snInit
( snProblem *prob,
  cint       n,
  cint       m,
  cint       ne,
  cint       negCon,
  cint       nnCon,
  cint       nnObj,
  cint       nnJac,
  cint       linear_objective_row,
  char      *PrintFileName,
  char      *SummaryFileName )
{
  cint       pre_alloc_len = 500;
  cint       i_sumry       = 0;
  cint       i_print       = 0;
  cint       nb;
  cint       iExit;

  prob->set_flags = 0;

  prob->dimension.n          =  n;
  prob->dimension.m          =  m;
  prob->dimension.ne         = ne;
  prob->dimension.negCon     = ( negCon < 1 ) ? 1 : negCon;
  prob->dimension.nnCon      = nnCon;
  prob->dimension.nnObj      = nnObj;
  prob->dimension.nnJac      = nnJac;
  prob->linear_objective_row = linear_objective_row;

  prob->error_msg = NULL;
  snClearLastError( prob );

  prob->option.print_file_name = NULL;
  prob->option.sumry_file_name = NULL;
  prob->option.ProblemName     = NULL;

  snSetPrintFileName  ( prob, PrintFileName   );
  snSetSummaryFileName( prob, SummaryFileName );
  snOpenFiles( prob );

  prob->rw = (double*)malloc( sizeof(double)*pre_alloc_len   );
  prob->iw = (cint*  )malloc( sizeof(cint  )*pre_alloc_len   );
  prob->cw = (char*  )malloc( sizeof(char  )*pre_alloc_len*8 );

  sninit_( &(prob->print_file_unit), &(prob->sumry_file_unit),
	   prob->cw, &pre_alloc_len,
	   prob->iw, &pre_alloc_len,
	   prob->rw, &pre_alloc_len,
	   pre_alloc_len ); /* size of character input passed to
				f77 subroutines */

  prob->length_cw = prob->length_iw = prob->length_rw = 500;

  /*  Get min{iw,rw,cw} size requirements
      Suppress printing                   */

  i_print =  prob->iw[11]; prob->iw[11] = 0;
  i_sumry =  prob->iw[12]; prob->iw[12] = 0;

  snmemb_
    ( &iExit,
      &(prob->dimension.m),
      &(prob->dimension.n),
      &(prob->dimension.ne),
      &(prob->dimension.negCon),
      &(prob->dimension.nnCon),
      &(prob->dimension.nnJac),
      &(prob->dimension.nnObj),
      &(prob->length_cw),
      &(prob->length_iw),
      &(prob->length_rw),
      prob->cw,
      &pre_alloc_len,
      prob->iw,
      &pre_alloc_len,
      prob->rw,
      &pre_alloc_len,
      pre_alloc_len );

  prob->iw[11] = i_print;
  prob->iw[12] = i_sumry;

  /* resize the workspace to the required value, copying the original
  content */

  prob->cw = (char  *)realloc( prob->cw , sizeof(char  )*prob->length_cw*8 );
  prob->iw = (cint  *)realloc( prob->iw , sizeof(cint  )*prob->length_iw   );
  prob->rw = (double*)realloc( prob->rw , sizeof(double)*prob->length_rw   );

  prob->iw[6] = prob->length_cw;
  prob->iw[4] = prob->length_iw;
  prob->iw[2] = prob->length_rw;

  nb         = prob->dimension.n + prob->dimension.m;

  prob->Jcol = (double*)malloc( sizeof(double) * prob->dimension.ne  );
  prob->indJ = (cint  *)malloc( sizeof(cint  ) * prob->dimension.ne  );
  prob->locJ = (cint  *)malloc( sizeof(cint  ) *(prob->dimension.n + 1));
  prob->bl   = (double*)malloc( sizeof(double) * nb );
  prob->bu   = (double*)malloc( sizeof(double) * nb );

  prob->pi   = (double*)malloc( sizeof(double) * prob->dimension.m );
  prob->hs   = (cint*  )malloc( sizeof(cint  ) * nb );
  prob->xs   = (double*)malloc( sizeof(double) * nb );
  prob->rc   = (double*)malloc( sizeof(double) * nb );

  return SN_OK;
}

cint snDelete( snProblem* prob )
{
  if ( prob->iw == NULL )
    /* in the event that snDelete gets called after it has already
       been invoked */
    return SN_OK;

  snClearLastError( prob );

  free(prob->iw);      prob->iw = NULL;
  free(prob->rw);      prob->rw = NULL;
  free(prob->cw);      prob->cw = NULL;

  free(prob->Jcol);    prob->Jcol = NULL;
  free(prob->indJ);    prob->indJ = NULL;
  free(prob->locJ);    prob->locJ = NULL;

  free(prob->pi);      prob->pi = NULL;

  free(prob->hs);      prob->hs = NULL;

  free(prob->bl);      prob->bl = NULL;
  free(prob->bu);      prob->bu = NULL;
  free(prob->xs);      prob->xs = NULL;
  free(prob->rc);      prob->rc = NULL;

  return SN_OK;
}

cint snCheck( snProblem* prob )
{
  /* use this subroutine to check the consistency among all the options.
     This routine is typically called once before snSolve().
   */

  if ( !( prob->set_flags & SN_BIT_MAJOR_OPTIMALITY_TOLERANCE ) )
    ; /* major_optimality_tolerance was not explicitly set */

  return SN_OK;
}
