/home/coin/SVN-release/OS-2.1.1/Couenne/src/disjunctive/OsiCuts2MatrVec.cpp

Go to the documentation of this file.
00001 /* $Id: OsiCuts2MatrVec.cpp 141 2009-06-03 04:19:19Z pbelotti $ */
00002 /*
00003  * Name:    OsiCuts2MatrVec.cpp
00004  * Author:  Pietro Belotti
00005  * Purpose: turn OsiCuts objects into coefficient matrix and rhs vector
00006  *
00007  * (C) Carnegie-Mellon University, 2008. 
00008  * This file is licensed under the Common Public License (CPL)
00009  */
00010 
00011 #include "CouenneDisjCuts.hpp"
00012 #include "OsiSolverInterface.hpp"
00013 #include "CoinPackedVector.hpp"
00014 #include "CoinPackedMatrix.hpp"
00015 #include "CoinHelperFunctions.hpp"
00016 #include "CouennePrecisions.hpp"
00017 #include "CouenneCutGenerator.hpp"
00018 #include "CouenneProblem.hpp"
00019 
00020 
00021 // add CGLP columns to solver interface; return number of columns
00022 // added (for later removal)
00023 int CouenneDisjCuts::OsiCuts2MatrVec (OsiSolverInterface *cglp,
00024                                       OsiCuts *cuts,
00025                                       int displRow,
00026                                       int displRhs) const {
00027   int
00028     ncols  = 0,
00029     nrcuts = cuts -> sizeRowCuts (),
00030     nccuts = cuts -> sizeColCuts (),
00031     ncgrow = cglp -> getNumRows () - 1,
00032     nnzR   = 0,
00033     ncC    = 0;
00034 
00035   if (!(nrcuts || nccuts))
00036     return 0;
00037 
00038   // count nonzero in row cuts
00039   for (int i=nrcuts; i--;) {
00040 
00041     OsiRowCut *cut = cuts -> rowCutPtr (i);
00042 
00043     if ((cut -> sense () == 'E') ||
00044         (cut -> sense () == 'R')) {
00045 
00046       nnzR += 2 * cut -> row (). getNumElements ();
00047       nrcuts++; // can increase as not used
00048 
00049     } else nnzR += cut -> row (). getNumElements ();
00050   }
00051 
00052   // count bound constraints in column cuts
00053   for (int i=nccuts; i--;) {
00054     OsiColCut *cut = cuts -> colCutPtr (i);
00055     ncC += 
00056       cut -> lbs ().getNumElements () + 
00057       cut -> ubs ().getNumElements ();
00058   }
00059 
00060   int 
00061      nnz     = 2 * (nnzR + 2*nrcuts + 3*ncC),
00062     *indices = new int [nnz],
00063     *start   = new int [nrcuts + ncC + 1],
00064      curel   = 0;
00065 
00066   double 
00067     *elements = new double [nnz],               // for row cuts + col cuts
00068     *collb    = new double [2*(nrcuts + ncC)],  // lower bounds for new columns
00069     *colub    = new double [2*(nrcuts + ncC)],  // upper 
00070     *obj      = new double [2*(nrcuts + ncC)];  // objective coefficient (zero)
00071 
00072   // trivial, lower/upper bounds and objective coefficients
00073   CoinFillN (collb, 2*(nrcuts + ncC), 0.);
00074   CoinFillN (colub, 2*(nrcuts + ncC), 1.);
00075   CoinFillN (obj,   2*(nrcuts + ncC), 0.);
00076 
00077   // scan OsiColCuts ////////////////////////////////////////
00078 
00079   double *saveEl  = elements;
00080   int    *saveInd = indices;
00081 
00082   for (int i = nccuts; i--;) {
00083 
00084     OsiColCut *cut = cuts -> colCutPtr (i);
00085 
00086     const CoinPackedVector
00087       &lbs = cut -> lbs (),
00088       &ubs = cut -> ubs ();
00089 
00090     // lower bounds
00091     const int
00092       *lind = lbs. getIndices (), 
00093        nlcc = lbs. getNumElements ();
00094     const double *lele = lbs. getElements ();
00095 
00096     for (int j = nlcc; j--; lind++, lele++)
00097 
00098       if (couenneCG_ -> Problem () -> Var (*lind) -> Multiplicity () > 0) {
00099         *start++ = curel;
00100         *elements++ = -1.;        *indices++ = displRow + *lind;
00101         if (fabs (*lele) > COUENNE_EPS) 
00102           {*elements++ = -*lele;     *indices++ = displRhs; curel++;}
00103         *elements++ =  1.;        *indices++ = ncgrow;
00104         curel += 2;
00105       }
00106 
00107     // upper bounds
00108     const int
00109       *uind = ubs. getIndices (), 
00110        nucc = ubs. getNumElements ();
00111     const double *uele = ubs. getElements ();
00112 
00113     for (int j = nucc; j--; uind++, uele++)
00114       if (couenneCG_ -> Problem () -> Var (*uind) -> Multiplicity () > 0) {
00115         *start++ = curel;
00116         *elements++ =  1.;        *indices++ = displRow + *uind;
00117         if (fabs (*uele) > COUENNE_EPS) 
00118           {*elements++ = *uele;     *indices++ = displRhs; curel++;}
00119         *elements++ =  1.;        *indices++ = ncgrow;
00120         curel += 2;
00121       }
00122 
00123     ncols += nlcc + nucc;
00124   }
00125 
00126   elements = saveEl;
00127   indices  = saveInd;
00128     
00129   //elements -= (3 * ncols);
00130   //indices  -= (3 * ncols);
00131   start    -= ncols;
00132 
00133   start [ncols] = curel; // may go
00134 
00135 
00136   if (jnlst_ -> ProduceOutput (J_MATRIX, J_DISJCUTS)) {
00137     printf ("%d cuts, have %d cols and cur el is %d. Now for the %d row cuts\n",
00138             nccuts, ncols, curel, nrcuts);
00139 
00140     printf ("matrix (.. %d) %d elements:\n", ncols, curel);
00141 
00142     printf ("start: "); for (int i=0; i<=ncols; i++) printf ("%d ", start [i]);
00143 
00144     printf ("\nElements:\n"); 
00145     for (int i=0, j=0; i<ncols; i++) {
00146       for (int k=0; k<start[i+1] - start[i]; k++, j++) 
00147         printf ("(%d %g) ", indices [j], elements [j]);
00148       printf ("\n");
00149     }
00150   }
00151 
00152   // scan OsiRowCuts /////////////////////////////////////////////
00153 
00154   for (int i = cuts -> sizeRowCuts (); i--;) {
00155 
00156     OsiRowCut *cut = cuts -> rowCutPtr (i);
00157 
00158     const CoinPackedVector &row = cut -> row ();
00159 
00160     const double 
00161       rhs    = cut -> rhs   (),
00162       rng    = cut -> range (),
00163       *rowEl = row. getElements ();
00164 
00165     const int 
00166       *rowIn = row. getIndices (),
00167       rowNE  = row. getNumElements ();
00168 
00169     switch (cut -> sense ()) {
00170 
00171     case 'L': 
00172 
00173       start [ncols++] = curel;
00174       CoinCopyDisp (rowIn, rowNE, indices  + curel, displRow);
00175       CoinCopyN    (rowEl, rowNE, elements + curel);
00176       curel += rowNE;
00177       if (fabs (rhs) > COUENNE_EPS) {indices  [curel] = displRhs;  elements [curel++] = rhs;}
00178       indices  [curel] = ncgrow;    elements [curel++] = 1.;
00179 
00180       break;
00181 
00182     case 'E':
00183     case 'R':
00184 
00185       start [ncols++] = curel;
00186       CoinCopyDisp (rowIn, rowNE, indices  + curel, displRow);
00187       CoinCopyN    (rowEl, rowNE, elements + curel);
00188       curel += rowNE;
00189       if (fabs (rhs+rng) > COUENNE_EPS) {indices  [curel] = displRhs; elements [curel++] = rhs + rng;}
00190       // rng only used here, zero for 'E'
00191       indices  [curel] = ncgrow;   elements [curel++] = 1.;
00192 
00193       start [ncols++] = curel;
00194       CoinCopyDisp (rowIn, rowNE, indices  + curel, displRow);
00195       CoinInvN     (rowEl, rowNE, elements + curel);
00196       curel += rowNE;
00197       if (fabs (rhs) > COUENNE_EPS) {indices  [curel] = displRhs;  elements [curel++] = -rhs;}
00198       indices  [curel] = ncgrow;    elements [curel++] = 1.;
00199 
00200       break;
00201 
00202     case 'G':
00203       start [ncols++] = curel;
00204       CoinCopyDisp (rowIn, rowNE, indices  + curel, displRow);
00205       CoinInvN     (rowEl, rowNE, elements + curel);
00206       curel += rowNE;
00207       if (fabs (rhs) > COUENNE_EPS) {indices  [curel] = displRhs;  elements [curel++] = -rhs;}
00208       indices  [curel] = ncgrow;    elements [curel++] = 1.;
00209 
00210       break;
00211 
00212     default: printf ("unknown type of cut\n");
00213       exit (-1);
00214     }
00215   }
00216 
00217   start [ncols] = curel;
00218 
00219   //  ncols=1;start[1]=3;
00220   /*  *start=0; start[1]=1;
00221   *indices=0;
00222   *elements=4.55;
00223   *collb=-3;
00224   *colub=+3;
00225   *obj=4;*/
00226 
00227 
00228   if (jnlst_ -> ProduceOutput (J_MATRIX, J_DISJCUTS)) {
00229     printf ("===================\nmatrix (.. %d) %d elements:\n", ncols, curel);
00230 
00231     printf ("start: "); for (int i=0; i<=ncols; i++) printf ("%d ", start [i]);
00232 
00233     printf ("\nElements:\n"); 
00234     for (int i=0, j=0; i<ncols; i++) {
00235       for (int k=0; k<start[i+1] - start[i]; k++, j++) 
00236         printf ("(%d %g) ", indices [j], elements [j]);
00237       printf ("\n");
00238     }
00239   }
00240 
00241   /*{
00242     const CoinPackedMatrix *m = cglp->getMatrixByCol();
00243     printf ("before: size_ = %d, start [%d] = %d\n", 
00244             m -> getNumElements (), 
00245             m -> getSizeVectorLengths(),
00246             m -> getVectorStarts () [m -> getSizeVectorLengths()]);
00247             }*/
00248 
00249   cglp -> addCols (ncols,    start,
00250                    indices,  elements,
00251                    collb,    colub,   
00252                    obj);
00253 
00254   /*{
00255     const CoinPackedMatrix *m = cglp->getMatrixByCol();
00256     printf ("after: size_ = %d, start [%d] = %d\n", 
00257             m -> getNumElements (), 
00258             m -> getSizeVectorLengths(),
00259             m -> getVectorStarts () [m -> getSizeVectorLengths()]);
00260             }*/
00261 
00262   delete [] elements;
00263   delete [] collb;
00264   delete [] colub;
00265   delete [] obj;
00266 
00267   delete [] indices;
00268   delete [] start;
00269 
00270   return ncols;
00271 }

Generated on Mon May 3 03:05:19 2010 by  doxygen 1.4.7