--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/DLD-FUNCTIONS/__glpk__.cc	Tue Mar 22 16:18:26 2005 +0000
@@ -0,0 +1,653 @@
+/*- --------------------------------------------------------------------
+ --
+ -- Copyright (C) 2005, Nicolo' Giorgetti, All rights reserved.
+ -- E-mail: <giorgetti@dii.unisi.it>.
+ --
+ -- This file is part of GLPKOCT an Octave interface to GLPK.
+ --
+ -- GLPK is free software; you can redistribute it and/or modify it
+ -- under the terms of the GNU General Public License as published by
+ -- the Free Software Foundation; either version 2, or (at your option)
+ -- any later version.
+ --
+ -- GLPK is distributed in the hope that it will be useful, but WITHOUT
+ -- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+ -- or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
+ -- License for more details.
+ --
+ -- You should have received a copy of the GNU General Public License
+ -- along with GLPK; see the file COPYING. If not, write to the Free
+ -- Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
+ -- 02111-1307, USA.
+ --
+ -- ---------------------------------------------------------------------*/
+
+#include <cfloat>
+#include <csetjmp>
+#include <ctime>
+
+//-- Octave headers
+#include <oct.h>
+#include <octave/ov-struct.h>
+#include <octave/config.h>
+#include <octave/error.h>
+
+//-- GLPK C header
+extern "C"{
+#include "glpk.h"
+}
+
+#define OCTOUT octave_stdout
+#define OCTERR octave_stdout
+#define NIntP 17
+#define NRealP 10
+
+int lpxIntParam[NIntP]= {
+   1,
+   1,
+   0,
+   1,
+   0,
+   -1,
+   0,
+   200,
+   1,
+   2,
+   0,
+   1,
+   0,
+   0,
+   2, 
+   2,
+   1
+};
+
+int IParam[NIntP]={
+   LPX_K_MSGLEV,
+   LPX_K_SCALE,
+   LPX_K_DUAL,
+   LPX_K_PRICE,
+   LPX_K_ROUND,
+   LPX_K_ITLIM,
+   LPX_K_ITCNT,
+   LPX_K_OUTFRQ,
+   LPX_K_MPSINFO,
+   LPX_K_MPSOBJ,
+   LPX_K_MPSORIG,
+   LPX_K_MPSWIDE,
+   LPX_K_MPSFREE,
+   LPX_K_MPSSKIP,
+   LPX_K_BRANCH,
+   LPX_K_BTRACK,
+   LPX_K_PRESOL
+};
+
+
+double lpxRealParam[NRealP]={
+   0.07,
+   1e-7,
+   1e-7,
+   1e-9,
+   -DBL_MAX,
+   DBL_MAX,
+   -1.0,
+   0.0,
+   1e-6,
+   1e-7
+};
+
+int RParam[NRealP]={
+   LPX_K_RELAX,
+   LPX_K_TOLBND,
+   LPX_K_TOLDJ,
+   LPX_K_TOLPIV,
+   LPX_K_OBJLL,
+   LPX_K_OBJUL,
+   LPX_K_TMLIM,
+   LPX_K_OUTDLY,
+   LPX_K_TOLINT,
+   LPX_K_TOLOBJ
+};
+
+jmp_buf mark;  //-- Address for long jump to jump to
+int     fperr; //-- Global error number
+
+
+int glpkmex_fault_hook(void *info,  char *msg)
+{    
+    OCTERR<<"*** SEVERE CRITICAL ERROR *** from GLPK !\n\n"<<msg<<" %s\n";
+    longjmp( mark, -1 );
+}
+
+int glpkmex_print_hook(void *info,  char *msg)
+{
+    OCTERR<<msg<<"\n";
+    return 1;
+}
+
+
+int glpk(int sense,int n, int m, double *c,int nz,int *rn,int *cn, 
+         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)
+{
+   
+   LPX *lp;
+   int i,j;
+   int error;
+   clock_t  t_start;
+   int typx=0;
+   int method;
+
+   t_start = clock();
+   
+   lib_set_fault_hook(NULL,glpkmex_fault_hook);
+   
+
+   if (lpxIntParam[0] > 1){
+         lib_set_print_hook(NULL,glpkmex_print_hook);
+   }
+      
+   lp=lpx_create_prob();
+     
+   
+   //-- Set the sense of optimization
+   if (sense==1) lpx_set_obj_dir(lp,LPX_MIN);
+   else lpx_set_obj_dir(lp,LPX_MAX);
+
+   //-- If the problem has integer structural variables switch to MIP
+   if(isMIP) lpx_set_class(lp,LPX_MIP);
+
+   lpx_add_cols(lp,n);   
+   for(i=0;i<n;i++){
+	
+      //-- Define type of the structural variables
+      if (!freeLB[i] && !freeUB[i]){
+        lpx_set_col_bnds(lp,i+1,LPX_DB,lb[i],ub[i]);
+      }else{
+         if (!freeLB[i] && freeUB[i]){
+            lpx_set_col_bnds(lp,i+1,LPX_LO,lb[i],ub[i]);
+         }else{
+            if (freeLB[i] && !freeUB[i]){
+               lpx_set_col_bnds(lp,i+1,LPX_UP,lb[i],ub[i]);
+            }else{
+               lpx_set_col_bnds(lp,i+1,LPX_FR,lb[i],ub[i]);
+            }
+         }
+      }
+      // -- Set the objective coefficient of the corresponding 
+      // -- structural variable. No constant term is assumed.
+	  lpx_set_obj_coef(lp,i+1,c[i]);
+
+      if(isMIP){
+        lpx_set_col_kind(lp,i+1,vartype[i]);
+      }
+   }
+   
+   lpx_add_rows(lp,m);   
+   for(i=0;i<m;i++){
+   
+      /* If the i-th row has no lower bound (types F,U), the
+         corrispondent parameter will be ignored.
+         If the i-th row has no upper bound (types F,L), the corrispondent
+         parameter will be ignored.
+         If the i-th row is of S type, the i-th LB is used, but
+         the i-th UB is ignored.
+      */
+      switch(ctype[i]){
+         case 'F': typx=LPX_FR; break;
+         case 'U': typx=LPX_UP; break;
+         case 'L': typx=LPX_LO; break;
+         case 'S': typx=LPX_FX; break;
+         case 'D': typx=LPX_DB; 
+      }
+      lpx_set_row_bnds(lp,i+1,typx,b[i],b[i]);
+
+   }
+   lpx_load_matrix(lp,nz,rn,cn,a);
+
+   if (save_pb){
+      if(lpx_write_cpxlp(lp, "outpb.lp") != 0){
+        OCTERR<<"Unable to write problem\n";
+        longjmp( mark, -1 );
+      }
+   }
+
+   //-- scale the problem data (if required)
+   //-- if (scale && (!presol || method == 1)) lpx_scale_prob(lp);
+   //-- LPX_K_SCALE=IParam[1]  LPX_K_PRESOL=IParam[16]
+   if (lpxIntParam[1] && (!lpxIntParam[16] || lpsolver!=1)){
+      lpx_scale_prob(lp);
+   }
+   //-- build advanced initial basis (if required)
+   if (lpsolver == 1 && !lpxIntParam[16]){
+      lpx_adv_basis(lp);
+   }
+ 
+   for(i=0;i<NIntP;i++){
+     lpx_set_int_parm(lp,IParam[i],lpxIntParam[i]);
+   }
+   for(i=0;i<NRealP;i++){
+     lpx_set_real_parm(lp,RParam[i],lpxRealParam[i]);
+   } 
+
+   if(lpsolver==1) method='S';
+   else method='T';
+
+   switch(method){
+   case 'S':
+     if(isMIP){
+       method='I';
+       error=lpx_simplex(lp);
+       error=lpx_integer(lp);
+     }else{
+       error=lpx_simplex(lp);
+     }
+     break;
+   case 'T':
+     error=lpx_interior(lp);
+     break;
+   default:
+     insist(method != method);
+   }
+
+   /*
+       error assumes the following results:
+       error=0 <=> No errors
+       error=1 <=> Iteration limit exceeded.
+       error=2 <=> Numerical problems with basis matrix.
+   */
+   if(error==LPX_E_OK){
+     if(isMIP){
+       *status=(double)lpx_mip_status(lp);
+       *fmin=lpx_mip_obj_val(lp);
+     }else{
+       if(lpsolver==1){
+         *status=(double)lpx_get_status(lp);
+         *fmin=lpx_get_obj_val(lp);
+       }else{
+         *status=(double)lpx_ipt_status(lp);
+         *fmin=lpx_ipt_obj_val(lp);
+       }
+     }
+     if(isMIP){
+         for(i=0;i<n;i++)  xmin[i]=lpx_mip_col_val(lp,i+1);
+     }else{
+      /* Primal values */
+      for(i=0;i<n;i++){
+         if(lpsolver==1) xmin[i]=lpx_get_col_prim(lp,i+1);
+         else xmin[i]=lpx_ipt_col_prim(lp,i+1);
+      }
+      /* Dual values */
+      for(i=0; i<m; i++){
+         if(lpsolver==1) lambda[i]=lpx_get_row_dual(lp,i+1);
+         else lambda[i]=lpx_ipt_row_dual(lp,i+1);
+      }
+      /* Reduced costs */
+      for(i=0; i<lpx_get_num_cols(lp); i++){
+         if(lpsolver==1) redcosts[i]=lpx_get_col_dual(lp,i+1);
+         else redcosts[i]=lpx_ipt_col_dual(lp,i+1);
+      }
+     }
+     *time=((double)(clock() - t_start))/CLOCKS_PER_SEC;
+     *mem=(double)lib_env_ptr()->mem_tpeak;      
+
+     lpx_delete_prob(lp);
+     return(0);
+   }
+   lpx_delete_prob(lp);
+   *status=(double)error;
+   return(error);
+}
+
+
+
+
+
+DEFUN_DLD(glpkoct, args, nlhs, "glpkoct: OCT interface for the GLPK library. Don't use glpkoct, use glpk.m instead")
+{
+   // The list of values to return.  See the declaration in oct-obj.h
+   octave_value_list retval;
+
+   int nrhs=args.length();
+	
+   if(nrhs<1){
+      OCTERR<<"Use the script glpk for the optimization\n";
+      return retval;
+   }	
+   
+   //-- 1st Input. Sense of optimization.
+   int sense;
+   double SENSE = args(0).scalar_value ();
+   if (SENSE>=0) sense=1;
+   else sense =-1;
+   
+   //-- 2nd Input. A column array containing the objective function
+   //--            coefficients.
+   int mrowsc=args(1).rows();
+   
+   Matrix C(args(1).matrix_value());
+   double *c=C.fortran_vec();
+   	
+   //-- 3rd Input. A matrix containing the constraints coefficients.
+   // If matrix A is NOT a sparse matrix
+   // if(!mxIsSparse(A_IN)){
+   int mrowsA=args(2).rows();	
+   Matrix A(args(2).matrix_value()); // get the matrix
+   Array<int> rn (mrowsA*mrowsc+1);
+   Array<int> cn (mrowsA*mrowsc+1);
+   ColumnVector a (mrowsA*mrowsc+1, 0.0);
+
+   int nz=0;
+   for(int i=0;i<mrowsA;i++){
+      for(int j=0;j<mrowsc;j++){
+        if(A(i,j) != 0){
+           nz++;
+	   rn(nz)=i+1;
+	   cn(nz)=j+1;
+	   a(nz)=A(i,j);
+        }
+      }
+   }
+// DON'T DELETE THIS PART... REPRESENTS THE SPARSE MATRICES MANIPULATION
+//	  }else{
+//	    int i,j;
+//	    int *jc,*ir;
+//	    double *pr;
+//	    int nelc,count,row;
+//
+//	    /* NOTE: nnz is the actual number of nonzeros and is stored as the
+//       last element of the jc array where the size of the jc array is the
+//       number of columns + 1 */
+//	    nz = *(mxGetJc(A_IN) + mrowsc);
+//	    jc = mxGetJc(A_IN);
+//	    ir = mxGetIr(A_IN);
+//	    pr = mxGetPr(A_IN);
+//
+//       rn=(int *)calloc(nz+1,sizeof(int));
+//	    cn=(int *)calloc(nz+1,sizeof(int));
+//	    a=(double *)calloc(nz+1,sizeof(double));
+//
+//       count=0; row=0;
+//	    for(i=1;i<=mrowsc;i++){
+//	      nelc=jc[i]-jc[i-1];
+//	      for(j=0;j<nelc;j++){
+//		      count++;
+//		      rn[count]=ir[row]+1;
+//		      cn[count]=i;
+//		      a[count]=pr[row];
+//		      row++;
+//	      }
+//	    }
+//	  }
+
+   //-- 4th Input. A column array containing the right-hand side value
+   //	           for each constraint in the constraint matrix.
+   Matrix B(args(3).matrix_value());
+   double *b=B.fortran_vec();
+   	
+   for(int i=0; i< mrowsA; i++){
+   	if (isinf(b[i]) && b[i]<0) b[i]=-octave_Inf;
+   	if (isinf(b[i]) && b[i]>0) b[i]=octave_Inf;
+   }	
+
+      	
+   //-- 5th Input. A column array containing the sense of each constraint
+   //--            in the constraint matrix.
+   charMatrix CTYPE(args(4).char_matrix_value());	
+	
+   char *ctype=CTYPE.fortran_vec();
+	
+   //-- 6th Input. An array of length mrowsc containing the lower
+   //--            bound on each of the variables.
+   Matrix LB(args(5).matrix_value());
+   		 		
+   double *lb=LB.fortran_vec();
+   
+   //-- LB argument, default: Free
+   Array<int> freeLB (mrowsc);
+   for(int i=0;i<mrowsc;i++){
+      if(isinf(lb[i])){
+         freeLB(i)=1;
+         lb[i]=-octave_Inf;
+      }else freeLB(i)=0;
+   }
+   
+   //-- 7th Input. An array of at least length numcols containing the upper
+   //--            bound on each of the variables.
+   Matrix UB(args(6).matrix_value());
+   		
+   double *ub=UB.fortran_vec();
+ 			
+   Array<int> freeUB (mrowsc);
+   for(int i=0;i<mrowsc;i++){
+      if(isinf(ub[i])){
+         freeUB(i)=1;
+	 ub[i]=octave_Inf;
+      }else freeUB(i)=0;
+   }
+		
+   //-- 8th Input. A column array containing the types of the variables.
+   charMatrix VTYPE(args(7).char_matrix_value());
+
+   Array<int> vartype (mrowsc);
+   int isMIP;
+   for (int i = 0; i < mrowsc ; i++){
+      if(VTYPE(i,0)=='I'){
+         isMIP=1;
+	 vartype(i)=LPX_IV;
+      }else{
+         vartype(i)=LPX_CV;
+      }
+   }
+
+   //-- 9th Input. A structure containing the control parameters.
+   Octave_map PARAM = args(8).map_value();
+   
+   //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+   //-- Integer parameters                                             
+   //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+   //-- Level of messages output by the solver
+   if(PARAM.contains("msglev")){
+      octave_value tmp=PARAM.contents(PARAM.seek("msglev"))(0);
+
+      double numtmp=tmp.scalar_value();
+      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]=(int) numtmp;
+   }	
+
+   //-- scaling option
+   if(PARAM.contains("scale")){
+      octave_value tmp=PARAM.contents(PARAM.seek("scale"))(0);
+      double numtmp=tmp.scalar_value(); 
+      if((numtmp != 0) && (numtmp != 1) && (numtmp != 2)){
+         OCTOUT<<"'scale' parameter must be only:\n\t0 - no scaling,\n\t1 - equilibration scaling,\n\t2 - geometric mean scaling\n";
+         return retval;
+   }
+      lpxIntParam[1]=(int) numtmp;
+   }
+
+   //-- Dual dimplex option
+   if(PARAM.contains("dual")){
+      octave_value tmp=PARAM.contents(PARAM.seek("dual"))(0);
+      double numtmp=tmp.scalar_value();
+      if((numtmp != 0) && (numtmp != 1)){
+         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]=(int) numtmp;
+   }
+
+   //-- Pricing option
+   if(PARAM.contains("price")){
+      octave_value tmp=PARAM.contents(PARAM.seek("price"))(0);
+      double numtmp=tmp.scalar_value();
+      if((numtmp != 0) && (numtmp != 1)){
+         OCTOUT<<"'price' parameter must be only:\n\t0 - textbook pricing,\n\t1 - steepest edge pricing [default]\n";
+	 return retval;
+      }
+      lpxIntParam[3]=(int) numtmp;
+   }
+
+   //-- Solution rounding option 
+   if(PARAM.contains("round")){
+      octave_value tmp=PARAM.contents(PARAM.seek("round"))(0);
+      double numtmp=tmp.scalar_value();
+      if((numtmp != 0) && (numtmp != 1)){
+         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]=(int) numtmp;
+   }
+
+   //-- Simplex iterations limit
+   if(PARAM.contains("itlim")){
+      octave_value tmp=PARAM.contents(PARAM.seek("itlim"))(0);
+      lpxIntParam[5]=(int) tmp.scalar_value(); }
+
+   //-- Simplex iterations count
+   if(PARAM.contains("itcnt")){
+      octave_value tmp=PARAM.contents(PARAM.seek("itcnt"))(0);
+      lpxIntParam[6]=(int) tmp.scalar_value(); }
+ 
+   //-- Output frequency, in iterations
+   if(PARAM.contains("outfrq")){
+      octave_value tmp=PARAM.contents(PARAM.seek("outfrq"))(0);
+      lpxIntParam[7]=(int) tmp.scalar_value(); }
+ 
+   //-- Branching heuristic option
+   if(PARAM.contains("branch")){
+      octave_value tmp=PARAM.contents(PARAM.seek("branch"))(0);
+      double numtmp=tmp.scalar_value();			
+      if((numtmp != 0) && (numtmp != 1) && (numtmp != 2)){
+         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]=(int) numtmp;
+   }
+
+   //-- Backtracking heuristic option
+   if(PARAM.contains("btrack")){
+      octave_value tmp=PARAM.contents(PARAM.seek("btrack"))(0);
+      double numtmp=tmp.scalar_value();	
+      if((numtmp != 0) && (numtmp != 1) && (numtmp != 2)){
+         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]=(int) numtmp;
+   }
+
+   //-- Presolver option
+   if(PARAM.contains("presol")){
+      octave_value tmp=PARAM.contents(PARAM.seek("presol"))(0);
+      double numtmp=tmp.scalar_value();
+      if((numtmp != 0) && (numtmp != 1)){
+         OCTOUT<<"'presol' parameter must be only:\n\t0 - LP presolver is ***NOT*** used,\n\t1 - LP presol is used\n";
+	 return retval;
+      }
+      lpxIntParam[16]=(int) numtmp;
+   }
+
+   //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+   //-- Real parameters                                                
+   //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+   //-- Ratio test option
+   if(PARAM.contains("relax")){
+      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_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_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_value tmp=PARAM.contents(PARAM.seek("tolpiv"))(0);
+      lpxRealParam[3]=tmp.scalar_value(); }
+		
+   if(PARAM.contains("objll")){
+      octave_value tmp=PARAM.contents(PARAM.seek("objll"))(0);
+      lpxRealParam[4]=tmp.scalar_value(); }
+		
+   if(PARAM.contains("objul")){
+      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(); }
+	
+
+   //-- 10th Input. If the problem is a LP problem you may select which solver
+   //--      use: RSM (Revised Simplex Method) or IPM (Interior Point Method).
+   //--      If the problem is a MIP problem this field will be ignored.
+   octave_value tmp=args(9).scalar_value();
+   int lpsolver = (int) tmp.scalar_value();
+
+   //-- 11th Input. Saves a copy of the problem if SAVE<>0.
+   tmp=args(10).scalar_value();
+   int save_pb = (tmp.scalar_value() != 0);
+
+   //-- Assign pointers to the output parameters
+   ColumnVector xmin (mrowsc);
+   ColumnVector fmin (1);
+   ColumnVector status (1);
+   ColumnVector lambda (mrowsA);
+   ColumnVector redcosts (mrowsc);
+   ColumnVector time (1);
+   ColumnVector mem (1);
+	
+   int jmpret = setjmp( mark );
+   if (jmpret==0){
+      int error=glpk(sense,mrowsc,mrowsA,c,nz,
+		     rn.fortran_vec(),cn.fortran_vec(),a.fortran_vec(),b,ctype,
+		     freeLB.fortran_vec(),lb,freeUB.fortran_vec(),ub,
+		     vartype.fortran_vec(),isMIP,lpsolver,save_pb,
+		     xmin.fortran_vec(),fmin.fortran_vec(),status.fortran_vec(),
+		     lambda.fortran_vec(),redcosts.fortran_vec(), 
+		     time.fortran_vec(),mem.fortran_vec());
+   }
+
+   // Set the output parameters
+   retval(0)=octave_value(xmin);  // returns xmin
+   retval(1)=octave_value(fmin);   // returns fmin
+   retval(2)=octave_value(status); // returns status
+	
+   // Extra informations
+   Octave_map extra("lambda",octave_value(lambda)); //returns lambda
+   extra.assign("redcosts",octave_value(redcosts)); //returns the reduced costs
+   extra.assign("time",octave_value(time)); //time to solve the optimization pb
+   extra.assign("mem",octave_value(mem));   //memory used
+
+   retval(3)=extra;
+	
+   return retval;
+}
+