octave-nkf: src/DLD-FUNCTIONS/__glpk_

comparison src/DLD-FUNCTIONS/glpk.cc @ 5240:2bb31f40c339

[project @ 2005-03-24 02:34:13 by jwe]

author	jwe
date	Thu, 24 Mar 2005 02:34:13 +0000
parents	d432b7809fe5
children	ab89f95de831

comparison

equal deleted inserted replaced

-:8a7224914e16
+:2bb31f40c339
 extern "C" {
 #include <glpk.h>
 }
-#define OCTOUT octave_stdout
-#define OCTERR octave_stdout
 #define NIntP 17
 #define NRealP 10
 int lpxIntParam[NIntP] = {
-1,
+0,
 1,
 0,
 1,
 0,
 -1,
 int
 glpk_fault_hook (void * /* info */, char *msg)
 {
-OCTERR << "*** SEVERE CRITICAL ERROR *** from GLPK !\n\n"<<msg<<" %s\n";
+error ("CRITICAL ERROR in GLPK: %s", msg);
 longjmp (mark, -1);
 }
 int
 glpk_print_hook (void * /* info */, char *msg)
 {
-OCTERR << msg << "\n";
+message (0, "%s\n", msg);
 return 1;
 }
 int
 double *a, double *b, char *ctype, int *freeLB, double *lb,
 int *freeUB, double *ub, int *vartype, int isMIP, int lpsolver,
 int save_pb, double *xmin, double *fmin, double *status,
 double *lambda, double *redcosts, double *time, double *mem)
 {
-int error;
+int errnum;
 int typx = 0;
 int method;
 clock_t t_start = clock();
 if (save_pb)
 {
 if (lpx_write_cpxlp (lp, "outpb.lp") != 0)
 	{
-	  OCTERR << "Unable to write problem\n";
+	  error ("__glpk__: unable to write problem");
 	  longjmp (mark, -1);
 	}
 }
 //-- scale the problem data (if required)
 case 'S':
 {
 	if (isMIP)
 	  {
 	    method = 'I';
-	    error = lpx_simplex (lp);
+	    errnum = lpx_simplex (lp);
-	    error = lpx_integer (lp);
+	    errnum = lpx_integer (lp);
 	  }
 	else
-	  error = lpx_simplex(lp);
+	  errnum = lpx_simplex(lp);
 }
 break;
 case 'T':
-error = lpx_interior(lp);
+errnum = lpx_interior(lp);
 break;
 default:
 insist (method != method);
 }
-/*  error assumes the following results:
+/*  errnum assumes the following results:
-error=0 <=> No errors
+errnum = 0 <=> No errors
-error=1 <=> Iteration limit exceeded.
+errnum = 1 <=> Iteration limit exceeded.
-error=2 <=> Numerical problems with basis matrix.
+errnum = 2 <=> Numerical problems with basis matrix.
 */
-if (error == LPX_E_OK)
+if (errnum == LPX_E_OK)
 {
 if (isMIP)
 	{
 	  *status = static_cast<double> (lpx_mip_status (lp));
 	  *fmin = lpx_mip_obj_val (lp);
 return 0;
 }
 lpx_delete_prob (lp);
-*status= static_cast<double> (error);
+*status= static_cast<double> (errnum);
-return error;
+return errnum;
 }
 #endif
+#define OCTAVE_GLPK_GET_REAL_PARAM(NAME, IDX) \
+do \
+{ \
+if (PARAM.contains (NAME)) \
+	{ \
+	  Cell tmp = PARAM.contents (NAME); \
+\
+if (! tmp.is_empty ()) \
+	    { \
+	      lpxRealParam[IDX] = tmp(0).scalar_value (); \
+\
+if (error_state) \
+		{ \
+		  error ("glpk: invalid value in param." NAME); \
+		  return retval; \
+		} \
+	    } \
+	  else \
+	    { \
+	      error ("glpk: invalid value in param." NAME); \
+	      return retval; \
+	    } \
+	} \
+} \
+while (0)
+#define OCTAVE_GLPK_GET_INT_PARAM(NAME, VAL) \
+do \
+{ \
+if (PARAM.contains (NAME)) \
+	{ \
+	  Cell tmp = PARAM.contents (NAME); \
+\
+if (! tmp.is_empty ()) \
+	    { \
+	      VAL = tmp(0).int_value (); \
+\
+if (error_state) \
+		{ \
+		  error ("glpk: invalid value in param." NAME); \
+		  return retval; \
+		} \
+	    } \
+	  else \
+	    { \
+	      error ("glpk: invalid value in param." NAME); \
+	      return retval; \
+	    } \
+	} \
+} \
+while (0)
 DEFUN_DLD (__glpk__, args, ,
 "__glpk__: internal interface for the GLPK library.\n\
 You should be using using glpk instead")
 {
 #if defined (HAVE_GLPK)
 int nrhs = args.length ();
-if (nrhs < 1)
+if (nrhs != 9)
 {
-OCTERR<<"Use the script glpk for the optimization\n";
+print_usage ("__glpk__");
 return retval;
 }
 //-- 1nd Input. A column array containing the objective function
 //--            coefficients.
 int mrowsc = args(0).rows();
 Matrix C (args(0).matrix_value ());
+if (error_state)
+{
+error ("__glpk__: invalid value of C");
+return retval;
+}
 double *c = C.fortran_vec ();
 //-- 2nd Input. A matrix containing the constraints coefficients.
 // If matrix A is NOT a sparse matrix
 // if(!mxIsSparse(A_IN)){
-int mrowsA = args(1).rows();
 Matrix A (args(1).matrix_value ()); // get the matrix
+if (error_state)
+{
+error ("__glpk__: invalid value of A");
+return retval;
+}
+int mrowsA = A.rows ();
 Array<int> rn (mrowsA*mrowsc+1);
 Array<int> cn (mrowsA*mrowsc+1);
 ColumnVector a (mrowsA*mrowsc+1, 0.0);
 volatile int nz = 0;
 //	  }
 //-- 3rd Input. A column array containing the right-hand side value
 //	           for each constraint in the constraint matrix.
 Matrix B (args(2).matrix_value ());
+if (error_state)
+{
+error ("__glpk__: invalid value of b");
+return retval;
+}
 double *b = B.fortran_vec ();
 //-- 4th Input. An array of length mrowsc containing the lower
 //--            bound on each of the variables.
 Matrix LB (args(3).matrix_value ());
+if (error_state)
+{
+error ("__glpk__: invalid value of lb");
+return retval;
+}
 double *lb = LB.fortran_vec ();
 //-- LB argument, default: Free
 Array<int> freeLB (mrowsc);
 for (int i = 0; i < mrowsc; i++)
 //-- 5th Input. An array of at least length numcols containing the upper
 //--            bound on each of the variables.
 Matrix UB (args(4).matrix_value ());
+if (error_state)
+{
+error ("__glpk__: invalid value of ub");
+return retval;
+}
 double *ub = UB.fortran_vec ();
 Array<int> freeUB (mrowsc);
 for (int i = 0; i < mrowsc; i++)
 {
 }
 //-- 6th Input. A column array containing the sense of each constraint
 //--            in the constraint matrix.
 charMatrix CTYPE (args(5).char_matrix_value ());
+if (error_state)
+{
+error ("__glpk__: invalid value of ctype");
+return retval;
+}
 char *ctype = CTYPE.fortran_vec ();
 //-- 7th Input. A column array containing the types of the variables.
 charMatrix VTYPE (args(6).char_matrix_value ());
+if (error_state)
+{
+error ("__glpk__: invalid value of vtype");
+return retval;
+}
 Array<int> vartype (mrowsc);
 volatile int isMIP = 0;
 for (int i = 0; i < mrowsc ; i++)
 {
 }
 //-- 8th Input. Sense of optimization.
 volatile int sense;
 double SENSE = args(7).scalar_value ();
+if (error_state)
+{
+error ("__glpk__: invalid value of sense");
+return retval;
+}
 if (SENSE >= 0)
 sense = 1;
 else
 sense = -1;
 //-- 9th Input. A structure containing the control parameters.
 Octave_map PARAM = args(8).map_value ();
+if (error_state)
+{
+error ("__glpk__: invalid value of param");
+return retval;
+}
 //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 //-- Integer parameters
 //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 //-- Level of messages output by the solver
-if (PARAM.contains ("msglev"))
+OCTAVE_GLPK_GET_INT_PARAM ("msglev", lpxIntParam[0]);
-{
+if (lpxIntParam[0] < 0 || lpxIntParam[0] > 3)
-octave_value tmp = PARAM.contents (PARAM.seek ("msglev"))(0);
+{
+error ("__glpk__: param.msglev must be 0 (no output [default]) or 1 (error messages only) or 2 (normal output) or 3 (full output)");
-double numtmp = tmp.scalar_value ();
+return retval;
-if (numtmp != 0 && numtmp != 1 && numtmp != 2 && numtmp != 3)
+}
-	{
-	  OCTOUT << "'msglev' parameter must be only:\n\t0 - no output,\n\t1 - error messages only),\n\t2 - normal output,\n\t3 - full output [default]\n";
-	  return retval;
-	}
-lpxIntParam[0] = static_cast<int> (numtmp);
-}
 //-- scaling option
-if (PARAM.contains ("scale"))
+OCTAVE_GLPK_GET_INT_PARAM ("scale", lpxIntParam[1]);
-{
+if (lpxIntParam[1] < 0 || lpxIntParam[1] > 2)
-octave_value tmp = PARAM.contents (PARAM.seek ("scale"))(0);
+{
-double numtmp = tmp.scalar_value ();
+error ("__glpk__: param.scale must be 0 (no scaling) or 1 (equilibration scaling [default]) or 2 (geometric mean scaling)");
-if (numtmp != 0 && numtmp != 1 && numtmp != 2)
+return retval;
-	{
-	  OCTOUT << "'scale' parameter must be only:\n\t0 - no scaling,\n\t1 - equilibration scaling,\n\t2 - geometric mean scaling\n";
-	  return retval;
-	}
-lpxIntParam[1] = static_cast<int> (numtmp);
 }
 //-- Dual dimplex option
-if (PARAM.contains ("dual"))
+OCTAVE_GLPK_GET_INT_PARAM ("dual", lpxIntParam[2]);
-{
+if (lpxIntParam[2] < 0 || lpxIntParam[2] > 1)
-octave_value tmp = PARAM.contents (PARAM.seek ("dual"))(0);
+{
-double numtmp = tmp.scalar_value ();
+error ("__glpk__: param.dual must be 0 (do NOT use dual simplex [default]) or 1 (use dual simplex)");
-if (numtmp != 0 && numtmp != 1)
+return retval;
-	{
-	  OCTOUT<<"'dual' parameter must be only:\n\t0 - do not use the dual simplex [default],\n\t1 - use dual simplex\n";
-	  return retval;
-	}
-lpxIntParam[2] = static_cast<int> (numtmp);
 }
 //-- Pricing option
-if (PARAM.contains ("price"))
+OCTAVE_GLPK_GET_INT_PARAM ("price", lpxIntParam[3]);
-{
+if (lpxIntParam[3] < 0 || lpxIntParam[3] > 1)
-octave_value tmp = PARAM.contents (PARAM.seek ("price"))(0);
+{
-double numtmp = tmp.scalar_value();
+error ("__glpk__: param.price must be 0 (textbook pricing) or 1 (steepest edge pricing [default])");
-if (numtmp != 0 && numtmp != 1)
+return retval;
-	{
+}
-	  OCTOUT << "'price' parameter must be only:\n\t0 - textbook pricing,\n\t1 - steepest edge pricing [default]\n";
-	  return retval;
-	}
-lpxIntParam[3] = static_cast<int> (numtmp);
-}
 //-- Solution rounding option
-if (PARAM.contains ("round"))
+OCTAVE_GLPK_GET_INT_PARAM ("round", lpxIntParam[4]);
-{
+if (lpxIntParam[4] < 0 || lpxIntParam[4] > 1)
-octave_value tmp = PARAM.contents (PARAM.seek ("round"))(0);
+{
-double numtmp = tmp.scalar_value ();
+error ("__glpk__: param.round must be 0 (report all primal and dual values [default]) or 1 (replace tiny primal and dual values by exact zero)");
-if (numtmp != 0 && numtmp != 1)
+return retval;
-	{
-	  OCTOUT << "'round' parameter must be only:\n\t0 - report all primal and dual values [default],\n\t1 - replace tiny primal and dual values by exact zero\n";
-	  return retval;
-	}
-lpxIntParam[4] = static_cast<int> (numtmp);
 }
 //-- Simplex iterations limit
-if (PARAM.contains ("itlim"))
+OCTAVE_GLPK_GET_INT_PARAM ("itlim", lpxIntParam[5]);
-{
-octave_value tmp = PARAM.contents (PARAM.seek ("itlim"))(0);
-lpxIntParam[5] = static_cast<int> (tmp.scalar_value ());
-}
 //-- Simplex iterations count
-if (PARAM.contains ("itcnt"))
+OCTAVE_GLPK_GET_INT_PARAM ("itcnt", lpxIntParam[6]);
-{
-octave_value tmp = PARAM.contents (PARAM.seek ("itcnt"))(0);
-lpxIntParam[6] = static_cast<int> (tmp.scalar_value ());
-}
 //-- Output frequency, in iterations
-if (PARAM.contains ("outfrq"))
+OCTAVE_GLPK_GET_INT_PARAM ("outfrq", lpxIntParam[7]);
-{
-octave_value tmp = PARAM.contents (PARAM.seek ("outfrq"))(0);
-lpxIntParam[7] = static_cast<int> (tmp.scalar_value ());
-}
 //-- Branching heuristic option
-if (PARAM.contains("branch"))
+OCTAVE_GLPK_GET_INT_PARAM ("branch", lpxIntParam[14]);
-{
+if (lpxIntParam[14] < 0 || lpxIntParam[14] > 2)
-octave_value tmp = PARAM.contents (PARAM.seek ("branch"))(0);
+{
-double numtmp = tmp.scalar_value ();
+error ("__glpk__: param.branch must be (MIP only) 0 (branch on first variable) or 1 (branch on last variable) or 2 (branch using a heuristic by Driebeck and Tomlin [default]");
-if (numtmp != 0 && numtmp != 1 && numtmp != 2)
+return retval;
-	{
+}
-	  OCTOUT << "'branch' parameter must be only (for MIP only):\n\t0 - branch on the first variable,\n\t1 - branch on the last variable,\n\t2 - branch using a heuristic by Driebeck and Tomlin [default]\n";
-	  return retval;
-	}
-lpxIntParam[14] = static_cast<int> (numtmp);
-}
 //-- Backtracking heuristic option
-if (PARAM.contains ("btrack"))
+OCTAVE_GLPK_GET_INT_PARAM ("btrack", lpxIntParam[15]);
-{
+if (lpxIntParam[15] < 0 || lpxIntParam[15] > 2)
-octave_value tmp = PARAM.contents (PARAM.seek ("btrack"))(0);
+{
-double numtmp = tmp.scalar_value ();
+error ("__glpk__: param.btrack must be (MIP only) 0 (depth first search) or 1 (breadth first search) or 2 (backtrack using the best projection heuristic [default]");
-if (numtmp != 0 && numtmp != 1 && numtmp != 2)
+return retval;
-	{
+}
-	  OCTOUT << "'btrack' parameter must be only (for MIP only):\n\t0 - depth first search,\n\t1 - breadth first search,\n\t2 - backtrack using the best projection heuristic\n";
-	  return retval;
-	}
-lpxIntParam[15] = static_cast<int> (numtmp);
-}
 //-- Presolver option
-if (PARAM.contains ("presol"))
+OCTAVE_GLPK_GET_INT_PARAM ("presol", lpxIntParam[16]);
-{
+if (lpxIntParam[16] < 0 || lpxIntParam[16] > 1)
-octave_value tmp = PARAM.contents (PARAM.seek ("presol"))(0);
+{
-double numtmp = tmp.scalar_value ();
+error ("__glpk__: param.presol must be 0 (do NOT use LP presolver) or 1 (use LP presolver [default])");
-if (numtmp != 0 && numtmp != 1)
+return retval;
-	{
-	  OCTOUT << "'presol' parameter must be only:\n\t0 - LP presolver is ***NOT*** used,\n\t1 - LP presol is used\n";
-	  return retval;
-	}
-lpxIntParam[16] = static_cast<int> (numtmp);
 }
 //-- LPsolver option
 volatile int lpsolver = 1;
-if (PARAM.contains ("lpsolver"))
+OCTAVE_GLPK_GET_INT_PARAM ("lpsolver", lpsolver);
-{
+if (lpsolver < 1 || lpsolver > 2)
-octave_value tmp = PARAM.contents (PARAM.seek ("lpsolver"))(0);
+{
-double numtmp = tmp.scalar_value ();
+error ("__glpk__: param.lpsolver must be 1 (simplex method) or 2 (interior point method)");
-if (numtmp != 1 && numtmp != 2)
+return retval;
-	{
-	  OCTOUT << "'lpsolver' parameter must be only:\n\t1 - simplex method,\n\t2 - interior point method\n";
-	  return retval;
-	}
-lpsolver = static_cast<int> (numtmp);
 }
 //-- Save option
 volatile int save_pb = 0;
-if (PARAM.contains ("save"))
+OCTAVE_GLPK_GET_INT_PARAM ("save", save_pb);
-{
+save_pb = save_pb != 0;
-octave_value tmp = PARAM.contents (PARAM.seek ("save"))(0);
-save_pb = (tmp.scalar_value () != 0);
-}
 //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 //-- Real parameters
 //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 //-- Ratio test option
-if (PARAM.contains ("relax"))
+OCTAVE_GLPK_GET_REAL_PARAM ("relax", 0);
-{
-octave_value tmp = PARAM.contents (PARAM.seek ("relax"))(0);
-lpxRealParam[0] = tmp.scalar_value ();
-}
 //-- Relative tolerance used to check if the current basic solution
 //-- is primal feasible
-if (PARAM.contains ("tolbnd"))
+OCTAVE_GLPK_GET_REAL_PARAM ("tolbnd", 1);
-{
-octave_value tmp = PARAM.contents (PARAM.seek ("tolbn"))(0);
-lpxRealParam[1] = tmp.scalar_value ();
-}
 //-- Absolute tolerance used to check if the current basic solution
 //-- is dual feasible
-if (PARAM.contains ("toldj"))
+OCTAVE_GLPK_GET_REAL_PARAM ("toldj", 2);
-{
-octave_value tmp = PARAM.contents (PARAM.seek ("toldj"))(0);
-lpxRealParam[2] = tmp.scalar_value();
-}
 //-- Relative tolerance used to choose eligible pivotal elements of
 //--	the simplex table in the ratio test
-if (PARAM.contains ("tolpiv"))
+OCTAVE_GLPK_GET_REAL_PARAM ("tolpiv", 3);
-{
-octave_value tmp = PARAM.contents (PARAM.seek ("tolpiv"))(0);
+OCTAVE_GLPK_GET_REAL_PARAM ("objll", 4);
-lpxRealParam[3] = tmp.scalar_value ();
-}
+OCTAVE_GLPK_GET_REAL_PARAM ("objul", 5);
-if (PARAM.contains ("objll"))
+OCTAVE_GLPK_GET_REAL_PARAM ("tmlim", 6);
-{
-octave_value tmp = PARAM.contents (PARAM.seek ("objll"))(0);
+OCTAVE_GLPK_GET_REAL_PARAM ("outdly", 7);
-lpxRealParam[4] = tmp.scalar_value ();
-}
+OCTAVE_GLPK_GET_REAL_PARAM ("tolint", 8);
-if (PARAM.contains ("objul"))
+OCTAVE_GLPK_GET_REAL_PARAM ("tolobj", 9);
-{
-octave_value tmp = PARAM.contents (PARAM.seek ("objul"))(0);
-lpxRealParam[5] = tmp.scalar_value ();
-}
-if (PARAM.contains ("tmlim"))
-{
-octave_value tmp = PARAM.contents (PARAM.seek ("tmlim"))(0);
-lpxRealParam[6] = tmp.scalar_value ();
-}
-if (PARAM.contains ("outdly"))
-{
-octave_value tmp = PARAM.contents (PARAM.seek ("outdly"))(0);
-lpxRealParam[7] = tmp.scalar_value ();
-}
-if (PARAM.contains ("tolint"))
-{
-octave_value tmp = PARAM.contents (PARAM.seek ("tolint"))(0);
-lpxRealParam[8] = tmp.scalar_value ();
-}
-if (PARAM.contains ("tolobj"))
-{
-octave_value tmp = PARAM.contents (PARAM.seek ("tolobj"))(0);
-lpxRealParam[9] = tmp.scalar_value ();
-}
 //-- Assign pointers to the output parameters
 ColumnVector xmin (mrowsc);
 ColumnVector fmin (1);
 ColumnVector status (1);

Mercurial > octave-nkf

comparison src/DLD-FUNCTIONS/__glpk__.cc @ 5240:2bb31f40c339

comparison src/DLD-FUNCTIONS/glpk.cc @ 5240:2bb31f40c339