Changeset 1034


Ignore:
Timestamp:
Jun 26, 2007 3:00:48 AM (12 years ago)
Author:
forrest
Message:

moving branches/devel to trunk

Location:
trunk
Files:
48 edited

Legend:

Unmodified
Added
Removed
  • trunk

    • Property svn:externals
      •  

        old new  
        11MSVisualStudio   https://projects.coin-or.org/svn/MSVisualStudio/trunk/ExternalsDirs/Clp
        22BuildTools    https://projects.coin-or.org/svn/BuildTools/trunk
        3 Data/Netlib   https://projects.coin-or.org/svn/Data/trunk/Netlib
        4 Data/Sample   https://projects.coin-or.org/svn/Data/trunk/Sample
        5 CoinUtils     https://projects.coin-or.org/svn/CoinUtils/stable/2.0/CoinUtils
         3ThirdParty/Blas https://projects.coin-or.org/svn/BuildTools/ThirdParty/Blas/stable/1.0
         4ThirdParty/Lapack https://projects.coin-or.org/svn/BuildTools/ThirdParty/Lapack/stable/1.0
         5Data/Netlib   https://projects.coin-or.org/svn/Data/stable/1.0/Netlib
         6Data/Sample   https://projects.coin-or.org/svn/Data/stable/1.0/Sample
         7CoinUtils     https://projects.coin-or.org/svn/CoinUtils/trunk/CoinUtils
  • trunk/Clp/src/CbcOrClpParam.cpp

    r883 r1034  
    546546    case SPECIALOPTIONS:
    547547      model->setSpecialOptions(value);
     548#ifdef COIN_HAS_CBC
    548549    case THREADS:
    549550      model->setNumberThreads(value);
    550551      break;
     552#endif
    551553    default:
    552554      break;
     
    578580    value=model->specialOptions();
    579581    break;
     582#ifndef COIN_HAS_CBC
    580583  case THREADS:
    581584    value = model->numberThreads();
     585#endif
    582586  default:
    583587    value=intValue_;
     
    754758    case LOGLEVEL:
    755759      oldValue = model.messageHandler()->logLevel();
    756       model.messageHandler()->setLogLevel(value);
     760      model.messageHandler()->setLogLevel(CoinAbs(value));
    757761      break;
    758762    case SOLVERLOGLEVEL:
     
    764768      model.setIntParam(CbcModel::CbcMaxNumNode,value);
    765769      break;
     770    case MAXSOLS:
     771      oldValue=model.getIntParam(CbcModel::CbcMaxNumSol);
     772      model.setIntParam(CbcModel::CbcMaxNumSol,value);
     773      break;
    766774    case STRONGBRANCHING:
    767775      oldValue=model.numberStrong();
     
    780788      model.setSizeMiniTree(value);
    781789      break;
     790    case CUTPASSINTREE:
     791      oldValue=model.getMaximumCutPasses();
     792      model.setMaximumCutPasses(value);
     793      break;
     794    case CUTPASS:
     795      oldValue=model.getMaximumCutPassesAtRoot();
     796      model.setMaximumCutPassesAtRoot(value);
     797      break;
     798#ifdef COIN_HAS_CBC
     799#ifdef CBC_THREAD
     800    case THREADS:
     801      oldValue=model.getNumberThreads();
     802      model.setNumberThreads(value);
     803      break;
     804#endif
     805#endif
    782806    default:
    783807      break;
     
    803827    value = model.getIntParam(CbcModel::CbcMaxNumNode);
    804828    break;
     829  case MAXSOLS:
     830    value = model.getIntParam(CbcModel::CbcMaxNumSol);
     831    break;
    805832  case STRONGBRANCHING:
    806833    value=model.numberStrong();
     
    815842    value=model.sizeMiniTree();
    816843    break;
     844  case CUTPASSINTREE:
     845    value=model.getMaximumCutPasses();
     846    break;
     847  case CUTPASS:
     848    value=model.getMaximumCutPassesAtRoot();
     849    break;
     850#ifdef COIN_HAS_CBC
     851#ifdef CBC_THREAD
     852  case THREADS:
     853    value = model.getNumberThreads();
     854#endif
     855#endif
    817856  default:
    818857    value=intValue_;
     
    11601199First just try with default settings and look carefully at the log file.  Did cuts help?  Did they take too long?  \
    11611200Look at output to see which cuts were effective and then do some tuning.  You will see that the \
    1162 options for cuts are off, on, root and ifmove.  Off is \
     1201options for cuts are off, on, root and ifmove, forceon.  Off is \
    11631202obvious, on means that this cut generator will be tried in the branch and cut tree (you can fine tune using \
    11641203'depth').  Root means just at the root node while 'ifmove' means that cuts will be used in the tree if they \
    1165 look as if they are doing some good and moving the objective value.  If pre-processing reduced the size of the \
     1204look as if they are doing some good and moving the objective value.  Forceon is same as on but forces code to use \
     1205cut generator at every node.  For probing forceonbut just does fixing probing in tree - not strengthening etc.  \
     1206If pre-processing reduced the size of the \
    11661207problem or strengthened many coefficients then it is probably wise to leave it on.  Switch off heuristics \
    11671208which did not provide solutions.  The other major area to look at is the search.  Hopefully good solutions \
     
    11921233#define REAL_BARRIER
    11931234#else
    1194   parameters[numberParameters-1].append("Uni!versityOfFlorida_dummy");
     1235  parameters[numberParameters-1].append("Uni!versityOfFlorida_dummy");   
    11951236#endif
    11961237#ifdef TAUCS_BARRIER
     
    12141255  parameters[numberParameters-1].append("root");
    12151256  parameters[numberParameters-1].append("ifmove");
     1257  parameters[numberParameters-1].append("forceOn");
    12161258  parameters[numberParameters-1].setLonghelp
    12171259    (
     
    12271269     "This switches on a heuristic which does branch and cut on the problem given by just \
    12281270using variables which have appeared in one or more solutions. \
    1229 It is obviously only tries after two or more solutions."
     1271It obviously only tries after two or more solutions."
    12301272     );
    12311273  parameters[numberParameters++]=
     
    12341276  parameters[numberParameters-1].append("pri!orities");
    12351277  parameters[numberParameters-1].append("column!Order?");
     1278  parameters[numberParameters-1].append("01f!irst?");
     1279  parameters[numberParameters-1].append("01l!ast?");
     1280  parameters[numberParameters-1].append("length!?");
    12361281  parameters[numberParameters-1].setLonghelp
    12371282    (
     
    12411286     );
    12421287#endif
    1243 #if 1
    12441288  parameters[numberParameters++]=
    12451289    CbcOrClpParam("cpp!Generate","Generates C++ code",
     
    125312974 bit in cbc generates size dependent code rather than computed values."
    12541298     );
    1255 #endif
    12561299#ifdef COIN_HAS_CLP
    12571300  parameters[numberParameters++]=
     
    13131356  parameters[numberParameters-1].append("root");
    13141357  parameters[numberParameters-1].append("ifmove");
     1358  parameters[numberParameters-1].append("forceOn");
    13151359  parameters[numberParameters-1].setLonghelp
    13161360    (
     
    13721416  parameters[numberParameters++]=
    13731417    CbcOrClpParam("dualize","Solves dual reformulation",
    1374                   0,1,DUALIZE,false);
     1418                  0,2,DUALIZE,false);
    13751419  parameters[numberParameters-1].setLonghelp
    13761420    (
     
    14571501 is initialized to 'default.mps'."
    14581502     );
     1503#ifdef COIN_HAS_CBC
     1504  parameters[numberParameters++]=
     1505    CbcOrClpParam("extra1","Extra integer parameter 1",
     1506                  -1,INT_MAX,EXTRA1,false);
     1507  parameters[numberParameters-1].setIntValue(-1);
     1508  parameters[numberParameters++]=
     1509    CbcOrClpParam("extra2","Extra integer parameter 2",
     1510                  -1,INT_MAX,EXTRA2,false);
     1511  parameters[numberParameters-1].setIntValue(-1);
     1512  parameters[numberParameters++]=
     1513    CbcOrClpParam("extra3","Extra integer parameter 3",
     1514                  -1,INT_MAX,EXTRA3,false);
     1515  parameters[numberParameters-1].setIntValue(-1);
     1516  parameters[numberParameters++]=
     1517    CbcOrClpParam("extra4","Extra integer parameter 4",
     1518                  -1,INT_MAX,EXTRA4,false);
     1519  parameters[numberParameters-1].setIntValue(-1);
     1520#endif
    14591521#ifdef COIN_HAS_CLP
    14601522  parameters[numberParameters++]=
     
    14871549    parameters[numberParameters-1].append("root");
    14881550    parameters[numberParameters-1].append("ifmove");
    1489   parameters[numberParameters-1].setLonghelp
     1551    parameters[numberParameters-1].append("forceOn");
     1552    parameters[numberParameters-1].setLonghelp
    14901553    (
    14911554     "This switches on flow cover cuts (either at root or in entire tree) \
     
    15181581  parameters[numberParameters-1].append("root");
    15191582  parameters[numberParameters-1].append("ifmove");
     1583  parameters[numberParameters-1].append("forceOn");
    15201584  parameters[numberParameters-1].setLonghelp
    15211585    (
     
    16411705  parameters[numberParameters-1].append("root");
    16421706  parameters[numberParameters-1].append("ifmove");
     1707  parameters[numberParameters-1].append("forceOn");
    16431708  parameters[numberParameters-1].setLonghelp
    16441709    (
    16451710     "This switches on knapsack cuts (either at root or in entire tree) \
     1711See branchAndCut for information on options."
     1712     );
     1713  parameters[numberParameters++]=
     1714    CbcOrClpParam("lift!AndProjectCuts","Whether to use Lift and Project cuts",
     1715                  "off",LANDPCUTS);
     1716  parameters[numberParameters-1].append("on");
     1717  parameters[numberParameters-1].append("root");
     1718  parameters[numberParameters-1].append("ifmove");
     1719  parameters[numberParameters-1].append("forceOn");
     1720  parameters[numberParameters-1].setLonghelp
     1721    (
     1722     "Lift and project cuts - may be expensive to compute. \
    16461723See branchAndCut for information on options."
    16471724     );
     
    16651742  parameters[numberParameters++]=
    16661743    CbcOrClpParam("log!Level","Level of detail in Coin branch and Cut output",
    1667                   -1,63,LOGLEVEL);
     1744                  -63,63,LOGLEVEL);
    16681745  parameters[numberParameters-1].setIntValue(1);
    16691746#endif
     
    17061783  parameters[numberParameters++]=
    17071784    CbcOrClpParam("maxN!odes","Maximum number of nodes to do",
    1708                   1,2147483647,MAXNODES);
     1785                  0,2147483647,MAXNODES);
    17091786  parameters[numberParameters-1].setLonghelp
    17101787    (
    17111788     "This is a repeatable way to limit search.  Normally using time is easier \
    17121789but then the results may not be repeatable."
     1790     );
     1791  parameters[numberParameters++]=
     1792    CbcOrClpParam("maxS!olutions","Maximum number of solutions to get",
     1793                  1,2147483647,MAXSOLS);
     1794  parameters[numberParameters-1].setLonghelp
     1795    (
     1796     "You may want to stop after (say) two solutions or an hour."
    17131797     );
    17141798#endif
     
    17351819  parameters[numberParameters-1].append("root");
    17361820  parameters[numberParameters-1].append("ifmove");
     1821  parameters[numberParameters-1].append("forceOn");
    17371822  parameters[numberParameters-1].setLonghelp
    17381823    (
     
    18211906     );
    18221907#ifdef COIN_HAS_CBC
     1908  parameters[numberParameters++]=
     1909    CbcOrClpParam("node!Strategy","What strategy to use to select nodes",
     1910                  "fewest",NODESTRATEGY);
     1911  parameters[numberParameters-1].append("depth");
     1912  parameters[numberParameters-1].append("upfewest");
     1913  parameters[numberParameters-1].append("downfewest");
     1914  parameters[numberParameters-1].append("updepth");
     1915  parameters[numberParameters-1].append("downdepth");
     1916  parameters[numberParameters-1].setLonghelp
     1917    (
     1918     "Normally before a solution the code will choose node with fewest infeasibilities. \
     1919You can choose depth as the criterion.  You can also say if up or down branch must \
     1920be done first (the up down choice will carry on after solution)."
     1921     );
    18231922  parameters[numberParameters++]=
    18241923    CbcOrClpParam("numberA!nalyze","Number of analysis iterations",
     
    18841983 more lightweight or do more if improvements are still being made.  As Presolve will return\
    18851984 if nothing is being taken out, you should not normally need to use this fine tuning."
    1886      );
     1985     );
     1986#endif
     1987#ifdef COIN_HAS_CBC
     1988  parameters[numberParameters++]=
     1989    CbcOrClpParam("passT!reeCuts","Number of cut passes in tree",
     1990                  -999999,999999,CUTPASSINTREE);
     1991  parameters[numberParameters-1].setLonghelp
     1992    (
     1993     "The default is one pass"
     1994     );
     1995#endif
     1996#ifdef COIN_HAS_CLP
    18871997  parameters[numberParameters++]=
    18881998    CbcOrClpParam("pertV!alue","Method of perturbation",
     
    19512061  parameters[numberParameters-1].append("sos");
    19522062  parameters[numberParameters-1].append("trysos");
     2063  parameters[numberParameters-1].append("equalall");
    19532064  parameters[numberParameters-1].append("strategy");
    19542065  parameters[numberParameters-1].setLonghelp
     
    19582069 Save option saves on file presolved.mps.  equal will turn <= cliques into \
    19592070==.  sos will create sos sets if all 0-1 in sets (well one extra is allowed) \
    1960 and no overlaps.  trysos is same but allows any number extra.  strategy is as \
     2071and no overlaps.  trysos is same but allows any number extra.  equalall will turn all \
     2072valid inequalities into equalities with integer slacks.  strategy is as \
    19612073on but uses CbcStrategy."
    19622074     );
     
    20702182  parameters[numberParameters-1].append("root");
    20712183  parameters[numberParameters-1].append("ifmove");
     2184  parameters[numberParameters-1].append("forceOn");
     2185  parameters[numberParameters-1].append("forceOnBut");
     2186  parameters[numberParameters-1].append("forceOnStrong");
     2187  parameters[numberParameters-1].append("forceOnButStrong");
    20722188  parameters[numberParameters-1].setLonghelp
    20732189    (
    20742190     "This switches on probing cuts (either at root or in entire tree) \
    2075 See branchAndCut for information on options."
    2076      );
     2191See branchAndCut for information on options. \
     2192but strong options do more probing"
     2193     );
     2194  parameters[numberParameters++]=
     2195    CbcOrClpParam("pumpT!une","Dubious ideas for feasibility pump",
     2196                  0,100000000,FPUMPTUNE);
     2197  parameters[numberParameters-1].setLonghelp
     2198    (
     2199     "This fine tunes Feasibility Pump \n\
     2200\t>=1000000 use as accumulate switch\n\
     2201\t>=1000 use index+1 as number of large loops\n\
     2202\t>=100 use 0.05 objvalue as increment\n\
     2203\t>=10 use +0.1 objvalue for cutoff (add)\n\
     2204\t1 == fix ints at bounds, 2 fix all integral ints, 3 and continuous at bounds"
     2205     );
     2206  parameters[numberParameters-1].setIntValue(0);
    20772207#endif
    20782208  parameters[numberParameters++]=
     
    21162246    parameters[numberParameters-1].append("root");
    21172247    parameters[numberParameters-1].append("ifmove");
    2118   parameters[numberParameters-1].setLonghelp
     2248    parameters[numberParameters-1].append("forceOn");
     2249    parameters[numberParameters-1].setLonghelp
    21192250    (
    21202251     "This switches on reduce and split  cuts (either at root or in entire tree) \
     2252See branchAndCut for information on options."
     2253     );
     2254  parameters[numberParameters++]=
     2255    CbcOrClpParam("residual!CapacityCuts","Whether to use Residual Capacity cuts",
     2256                  "off",RESIDCUTS);
     2257  parameters[numberParameters-1].append("on");
     2258  parameters[numberParameters-1].append("root");
     2259  parameters[numberParameters-1].append("ifmove");
     2260  parameters[numberParameters-1].append("forceOn");
     2261  parameters[numberParameters-1].setLonghelp
     2262    (
     2263     "Residual capacity cuts. \
    21212264See branchAndCut for information on options."
    21222265     );
     
    21502293#endif
    21512294#ifdef COIN_HAS_CBC
     2295  parameters[numberParameters++]=
     2296      CbcOrClpParam("Rins","Whether to try Relaxed Induced Neighborhood Search",
     2297                    "off",RINS);
     2298    parameters[numberParameters-1].append("on");
     2299  parameters[numberParameters-1].setLonghelp
     2300    (
     2301     "This switches on Relaxed induced neighborhood Search."
     2302     );
    21522303  parameters[numberParameters++]=
    21532304    CbcOrClpParam("round!ingHeuristic","Whether to use Rounding heuristic",
     
    21782329 is initialized to 'solution.file'.  To read the file use fread(int) twice to pick up number of rows \
    21792330and columns, then fread(double) to pick up objective value, then pick up row activities, row duals, column \
    2180 activities and reduced costs - see bottom of CbcOrClpParam.cpp for code that writes file."
     2331activities and reduced costs - see bottom of CbcOrClpParam.cpp for code that reads or writes file."
    21812332     );
    21822333  parameters[numberParameters++]=
     
    22892440  parameters[numberParameters++]=
    22902441    CbcOrClpParam("sprint!Crash","Whether to try sprint crash",
    2291                   -1,500,SPRINT);
     2442                  -1,5000000,SPRINT);
    22922443  parameters[numberParameters-1].setLonghelp
    22932444    (
     
    23442495 elements may increase."
    23452496     );
     2497#endif
     2498#ifdef COIN_HAS_CBC
     2499  parameters[numberParameters++]=
     2500    CbcOrClpParam("testO!si","Test OsiObject stuff",
     2501                  -1,INT_MAX,TESTOSI,false);
     2502#endif
     2503#ifdef CBC_THREAD
    23462504  parameters[numberParameters++]=
    23472505    CbcOrClpParam("thread!s","Number of threads to try and use",
    2348                   -2,64,THREADS,false);
     2506                  -100,10000,THREADS,false);
     2507  parameters[numberParameters-1].setLonghelp
     2508    (
     2509     "To use multiple threads, set threads to number wanted.  It may be better \
     2510to use one or two more than number of cpus available.  If 100+n then n threads and \
     2511threads used in sub-trees, if 200+n use threads for root cuts, 300+n - both."
     2512     );
    23492513#endif
    23502514#ifdef COIN_HAS_CBC
     
    23522516    CbcOrClpParam("tighten!Factor","Tighten bounds using this times largest \
    23532517activity at continuous solution",
    2354                   1.0,1.0e20,TIGHTENFACTOR,false);
     2518                  1.0e-3,1.0e20,TIGHTENFACTOR,false);
    23552519  parameters[numberParameters-1].setLonghelp
    23562520    (
     
    23742538#endif
    23752539#ifdef COIN_HAS_CBC
     2540  parameters[numberParameters++]=
     2541    CbcOrClpParam("tune!PreProcess","Dubious tuning parameters",
     2542                  0,20000000,PROCESSTUNE,false);
     2543  parameters[numberParameters-1].setLonghelp
     2544    (
     2545     "For making equality cliques this is minimumsize.  Also for adding \
     2546integer slacks.  May be used for more later \
     2547If <1000 that is what it does.  If <1000000 - numberPasses is (value/1000)-1 and tune is tune %1000. \
     2548If >= 1000000! - numberPasses is (value/1000000)-1 and tune is tune %1000000.  In this case if tune is now still >=10000 \
     2549numberPassesPerInnerLoop is changed from 10 to (tune-10000)-1 and tune becomes tune % 10000!!!!! - happy? - \
     2550so to keep normal limit on cliques of 5, do 3 major passes (include presolves) but only doing one tightening pass per major pass - \
     2551you would use 3010005 (I think)"
     2552     );
    23762553  parameters[numberParameters++]=
    23772554    CbcOrClpParam("two!MirCuts","Whether to use Two phase Mixed Integer Rounding cuts",
     
    23802557  parameters[numberParameters-1].append("root");
    23812558  parameters[numberParameters-1].append("ifmove");
     2559  parameters[numberParameters-1].append("forceOn");
    23822560  parameters[numberParameters-1].setLonghelp
    23832561    (
     
    24112589     );
    24122590#endif
     2591  parameters[numberParameters++]=
     2592    CbcOrClpParam("vector","Whether to use vector? Form of matrix in simplex",
     2593                  "off",VECTOR,7,false);
     2594  parameters[numberParameters-1].append("on");
     2595  parameters[numberParameters-1].setLonghelp
     2596    (
     2597     "If this is on and ClpPackedMatrix uses extra column copy in odd format."
     2598     );
    24132599  parameters[numberParameters++]=
    24142600    CbcOrClpParam("verbose","Switches on longer help on single ?",
     
    24582644  }
    24592645}
    2460 #endif
     2646/* Restore a solution from file.
     2647   mode 0 normal, 1 swap rows and columns and primal and dual
     2648   if 2 set then also change signs
     2649*/
     2650void restoreSolution(ClpSimplex * lpSolver,std::string fileName,int mode)
     2651{
     2652  FILE * fp=fopen(fileName.c_str(),"rb");
     2653  if (fp) {
     2654    int numberRows=lpSolver->numberRows();
     2655    int numberColumns=lpSolver->numberColumns();
     2656    int numberRowsFile;
     2657    int numberColumnsFile;
     2658    double objectiveValue;
     2659    fread(&numberRowsFile,sizeof(int),1,fp);
     2660    fread(&numberColumnsFile,sizeof(int),1,fp);
     2661    fread(&objectiveValue,sizeof(double),1,fp);
     2662    double * dualRowSolution = lpSolver->dualRowSolution();
     2663    double * primalRowSolution = lpSolver->primalRowSolution();
     2664    double * dualColumnSolution = lpSolver->dualColumnSolution();
     2665    double * primalColumnSolution = lpSolver->primalColumnSolution();
     2666    if (mode) {
     2667      // swap
     2668      int k=numberRows;
     2669      numberRows=numberColumns;
     2670      numberColumns=k;
     2671      double * temp;
     2672      temp = dualRowSolution;
     2673      dualRowSolution = primalColumnSolution;
     2674      primalColumnSolution=temp;
     2675      temp = dualColumnSolution;
     2676      dualColumnSolution = primalRowSolution;
     2677      primalRowSolution=temp;
     2678    }
     2679    if (numberRows!=numberRowsFile||numberColumns!=numberColumnsFile) {
     2680      std::cout<<"Mismatch on rows and/or columns"<<std::endl;
     2681    } else {
     2682      lpSolver->setObjectiveValue(objectiveValue);
     2683      fread(primalRowSolution,sizeof(double),numberRows,fp);
     2684      fread(dualRowSolution,sizeof(double),numberRows,fp);
     2685      fread(primalColumnSolution,sizeof(double),numberColumns,fp);
     2686      fread(dualColumnSolution,sizeof(double),numberColumns,fp);
     2687      if (mode==3) {
     2688        int i;
     2689        for (i=0;i<numberRows;i++) {
     2690          primalRowSolution[i] = -primalRowSolution[i];
     2691          dualRowSolution[i] = -dualRowSolution[i];
     2692        }
     2693        for (i=0;i<numberColumns;i++) {
     2694          primalColumnSolution[i] = -primalColumnSolution[i];
     2695          dualColumnSolution[i] = -dualColumnSolution[i];
     2696        }
     2697      }
     2698    }
     2699    fclose(fp);
     2700  } else {
     2701    std::cout<<"Unable to open file "<<fileName<<std::endl;
     2702  }
     2703}
     2704#endif
  • trunk/Clp/src/CbcOrClpParam.hpp

    r800 r1034  
    6060    MAXFACTOR,PERTVALUE,MAXITERATION,PRESOLVEPASS,IDIOT,SPRINT,
    6161    OUTPUTFORMAT,SLPVALUE,PRESOLVEOPTIONS,PRINTOPTIONS,SPECIALOPTIONS,
    62     SUBSTITUTION,DUALIZE,VERBOSE,THREADS,CPP,
     62    SUBSTITUTION,DUALIZE,VERBOSE,CPP,PROCESSTUNE,
    6363
    6464    STRONGBRANCHING=151,CUTDEPTH, MAXNODES,NUMBERBEFORE,NUMBERANALYZE,
    65     NUMBERMINI,MIPOPTIONS,MOREMIPOPTIONS,MAXHOTITS,FPUMPITS,CUTPASS,
     65    NUMBERMINI,MIPOPTIONS,MOREMIPOPTIONS,MAXHOTITS,FPUMPITS,MAXSOLS,
     66    FPUMPTUNE,TESTOSI,EXTRA1,EXTRA2,EXTRA3,EXTRA4,CUTPASSINTREE,
     67    THREADS,CUTPASS,
    6668#ifdef COIN_HAS_CBC
    6769    LOGLEVEL ,
     
    7072    DIRECTION=201,DUALPIVOT,SCALING,ERRORSALLOWED,KEEPNAMES,SPARSEFACTOR,
    7173    PRIMALPIVOT,PRESOLVE,CRASH,BIASLU,PERTURBATION,MESSAGES,AUTOSCALE,
    72     CHOLESKY,KKT,BARRIERSCALE,GAMMA,CROSSOVER,PFI,INTPRINT,
     74    CHOLESKY,KKT,BARRIERSCALE,GAMMA,CROSSOVER,PFI,INTPRINT,VECTOR,
    7375   
    7476    NODESTRATEGY = 251,BRANCHSTRATEGY,CUTSSTRATEGY,HEURISTICSTRATEGY,
     
    7678    ROUNDING,SOLVER,CLIQUECUTS,COSTSTRATEGY,FLOWCUTS,MIXEDCUTS,
    7779    TWOMIRCUTS,PREPROCESS,FPUMP,GREEDY,COMBINE,LOCALTREE,USESOLUTION,SOS,
     80    LANDPCUTS,RINS,RESIDCUTS,
    7881   
    7982    DIRECTORY=301,IMPORT,EXPORT,RESTORE,SAVE,DUALSIMPLEX,PRIMALSIMPLEX,EITHERSIMPLEX,
     
    8184    TIGHTEN,FAKEBOUND,HELP,PLUSMINUS,NETWORK,ALLSLACK,REVERSE,BARRIER,NETLIB_BARRIER,NETLIB_TUNE,
    8285    REALLY_SCALE,BASISIN,BASISOUT,SOLVECONTINUOUS,CLEARCUTS,VERSION,STATISTICS,DEBUG,DUMMY,PRINTMASK,
    83     OUTDUPROWS,USERCLP,
     86    OUTDUPROWS,USERCLP,MODELIN,
    8487
    8588    BAB=351,MIPLIB,STRENGTHEN,PRIORITYIN,USERCBC,
  • trunk/Clp/src/ClpCholeskyDense.cpp

    r894 r1034  
    724724  }
    725725}
     726#ifdef CHOL_USING_BLAS
     727// using simple blas interface
     728extern "C"
     729{
     730  /** BLAS Fortran subroutine DGEMM. */
     731  void F77_FUNC(dgemm,DGEMM)(char* transa, char* transb,
     732                             ipfint *m, ipfint *n, ipfint *k,
     733                             const double *alpha, const double *a, ipfint *lda,
     734                             const double *b, ipfint *ldb, const double *beta,
     735                             double *c, ipfint *ldc,
     736                             int transa_len, int transb_len);
     737}
     738#endif
    726739/* Non leaf recursive rectangle rectangle update,
    727740   nUnder is number of rows in iBlock,
     
    734747            int numberBlocks)
    735748{
     749#ifdef CHOL_USING_BLAS
     750    recRecLeaf(above , aUnder ,  aOther, diagonal,work, nUnderK);
     751void
     752ClpCholeskyDense::recRecLeaf(longDouble * above,
     753                             longDouble * aUnder, longDouble *aOther, longDouble * diagonal, longDouble * work,
     754                             int nUnder)
     755  aa = aOther-BLOCK;
     756  for (j = 0; j < BLOCK; j ++) {
     757    aa+=BLOCK;
     758    for (i=0;i< nUnderK;i++) {
     759      longDouble t00=aa[i+0*BLOCK];
     760      for (k=0;k<BLOCK;k++) {
     761        longDouble a00=aUnder[i+k*BLOCK]*work[k];
     762        t00 -= a00 * above[j + k * BLOCK];
     763      }
     764      aa[i]=t00;
     765    }
     766  return;
     767#endif
    736768  if (nDo<=BLOCK&&nUnder<=BLOCK&&nUnderK<=BLOCK) {
     769    assert (nDo==BLOCK&&nUnder==BLOCK);
    737770    recRecLeaf(above , aUnder ,  aOther, diagonal,work, nUnderK);
    738771  } else if (nDo<=nUnderK&&nUnder<=nUnderK) {
  • trunk/Clp/src/ClpCholeskyWssmp.cpp

    r994 r1034  
     1#ifdef WSSMP_BARRIER
    12// Copyright (C) 2003, International Business Machines
    23// Corporation and others.  All Rights Reserved.
     
    1314#include "ClpMessage.hpp"
    1415
    15 #ifdef WSSMP_BARRIER
    1616//#############################################################################
    1717// Constructors / Destructor / Assignment
  • trunk/Clp/src/ClpCholeskyWssmpKKT.cpp

    r994 r1034  
     1#ifdef WSSMP_BARRIER
    12// Copyright (C) 2004, International Business Machines
    23// Corporation and others.  All Rights Reserved.
     
    1314#include "ClpMessage.hpp"
    1415
    15 #ifdef WSSMP_BARRIER
    1616//#############################################################################
    1717// Constructors / Destructor / Assignment
  • trunk/Clp/src/ClpDualRowSteepest.cpp

    r800 r1034  
    364364#endif
    365365  assert (input->packedMode());
    366   assert (updatedColumn->packedMode());
     366  if (!updatedColumn->packedMode()) {
     367    // I think this means empty
     368#ifdef COIN_DEVELOP
     369    printf("updatedColumn not packed mode ClpDualRowSteepest::updateWeights\n");
     370#endif
     371    return 0.0;
     372  }
    367373  double alpha=0.0;
    368374  if (!model_->factorization()->networkBasis()) {
  • trunk/Clp/src/ClpDynamicMatrix.cpp

    r800 r1034  
    230230  numberElements_ = CoinMin(numberElements_,numberElements)+originalMatrix->getNumElements();
    231231  matrix_ = originalMatrix;
    232   zeroElements_ = false;
     232  flags_ &= ~1;
    233233  // resize model (matrix stays same)
    234234  int newRowSize = numberRows+CoinMin(numberSets_,CoinMax(frequency,numberRows))+1;
  • trunk/Clp/src/ClpEventHandler.cpp

    r754 r1034  
    6767  model_= model;
    6868}
     69//#############################################################################
     70// Constructors / Destructor / Assignment
     71//#############################################################################
     72
     73//-------------------------------------------------------------------
     74// Default Constructor
     75//-------------------------------------------------------------------
     76ClpDisasterHandler::ClpDisasterHandler (ClpSimplex * model) :
     77  model_(model)
     78{
     79
     80}
     81
     82//-------------------------------------------------------------------
     83// Copy constructor
     84//-------------------------------------------------------------------
     85ClpDisasterHandler::ClpDisasterHandler (const ClpDisasterHandler & rhs)
     86  : model_(rhs.model_)
     87
     88}
     89
     90//-------------------------------------------------------------------
     91// Destructor
     92//-------------------------------------------------------------------
     93ClpDisasterHandler::~ClpDisasterHandler ()
     94{
     95}
     96
     97//----------------------------------------------------------------
     98// Assignment operator
     99//-------------------------------------------------------------------
     100ClpDisasterHandler &
     101ClpDisasterHandler::operator=(const ClpDisasterHandler& rhs)
     102{
     103  if (this != &rhs) {
     104    model_ = rhs.model_;
     105  }
     106  return *this;
     107}
     108/* set model. */
     109void
     110ClpDisasterHandler::setSimplex(ClpSimplex * model)
     111{
     112  model_= model;
     113}
    69114
    70115
  • trunk/Clp/src/ClpEventHandler.hpp

    r754 r1034  
    8686  //@}
    8787};
     88/** Base class for Clp disaster handling
     89   
     90This is here to allow for disaster handling.  By disaster I mean that Clp
     91would otherwise give up
     92
     93*/
     94
     95class ClpDisasterHandler  {
     96 
     97public:
     98  /**@name Virtual methods that the derived classe should provide.
     99  */
     100  //@{
     101  /// Into simplex
     102  virtual void intoSimplex()=0;
     103  /// Checks if disaster
     104  virtual bool check() const = 0;
     105  /// saves information for next attempt
     106  virtual void saveInfo() =0;
     107  //@}
     108 
     109 
     110  /**@name Constructors, destructor */
     111
     112  //@{
     113  /** Default constructor. */
     114  ClpDisasterHandler(ClpSimplex * model = NULL);
     115  /** Destructor */
     116  virtual ~ClpDisasterHandler();
     117  // Copy
     118  ClpDisasterHandler(const ClpDisasterHandler&);
     119  // Assignment
     120  ClpDisasterHandler& operator=(const ClpDisasterHandler&);
     121  /// Clone
     122  virtual ClpDisasterHandler * clone() const =0;
     123
     124  //@}
     125 
     126  /**@name Sets/gets */
     127
     128  //@{
     129  /** set model. */
     130  void setSimplex(ClpSimplex * model);
     131  /// Get model
     132  inline ClpSimplex * simplex() const
     133  { return model_;};
     134  //@}
     135 
     136 
     137protected:
     138  /**@name Data members
     139     The data members are protected to allow access for derived classes. */
     140  //@{
     141  /// Pointer to simplex
     142  ClpSimplex * model_;
     143  //@}
     144};
    88145#endif
  • trunk/Clp/src/ClpGubDynamicMatrix.cpp

    r800 r1034  
    168168  numberElements_ = CoinMin(numberElements_,numberElements)+originalMatrix->getNumElements();
    169169  matrix_ = originalMatrix;
    170   zeroElements_ = false;
     170  flags_ &= ~1;
    171171  // resize model (matrix stays same)
    172172  model->resize(numberRows,numberNeeded);
  • trunk/Clp/src/ClpGubMatrix.cpp

    r800 r1034  
    40734073// Correct sequence in and out to give true value
    40744074void
    4075 ClpGubMatrix::correctSequence(int & sequenceIn, int & sequenceOut) const
     4075ClpGubMatrix::correctSequence(const ClpSimplex * model,int & sequenceIn, int & sequenceOut)
    40764076{
    4077   sequenceIn = trueSequenceIn_;
    4078   sequenceOut = trueSequenceOut_;
     4077  if (sequenceIn!=-999) {
     4078    sequenceIn = trueSequenceIn_;
     4079    sequenceOut = trueSequenceOut_;
     4080  }
    40794081}
  • trunk/Clp/src/ClpGubMatrix.hpp

    r754 r1034  
    153153  virtual int synchronize(ClpSimplex * model,int mode);
    154154  /// Correct sequence in and out to give true value
    155   virtual void correctSequence(int & sequenceIn, int & sequenceOut) const;
     155  virtual void correctSequence(const ClpSimplex * model,int & sequenceIn, int & sequenceOut) ;
    156156  //@}
    157157
  • trunk/Clp/src/ClpLinearObjective.cpp

    r754 r1034  
    198198    return 0.0;
    199199  }
     200}
     201// Return objective value (without any ClpModel offset) (model may be NULL)
     202double
     203ClpLinearObjective::objectiveValue(const ClpSimplex * model, const double * solution) const
     204{
     205  const double * cost = objective_;
     206  if (model&&model->costRegion())
     207    cost = model->costRegion();
     208  double currentObj=0.0;
     209  for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
     210    currentObj += cost[iColumn]*solution[iColumn];
     211  }
     212  return currentObj;
    200213}
    201214//-------------------------------------------------------------------
  • trunk/Clp/src/ClpLinearObjective.hpp

    r754 r1034  
    4545                            double & predictedObj,
    4646                            double & thetaObj);
     47  /// Return objective value (without any ClpModel offset) (model may be NULL)
     48  virtual double objectiveValue(const ClpSimplex * model, const double * solution) const ;
    4749  /// Resize objective
    4850  virtual void resize(int newNumberColumns) ;
  • trunk/Clp/src/ClpMain.cpp

    r828 r1034  
    6060
    6161int mainTest (int argc, const char *argv[],int algorithm,
    62               ClpSimplex empty, bool doPresolve,int switchOff);
     62              ClpSimplex empty, bool doPresolve,int switchOff,bool doVector);
    6363static void statistics(ClpSimplex * originalModel, ClpSimplex * model);
    6464static void generateCode(const char * fileName,int type);
     
    8686    int presolveOptions=0;
    8787    int doCrash=0;
     88    int doVector=0;
    8889    int doSprint=-1;
    8990    // set reasonable defaults
     
    475476              doCrash=action;
    476477              break;
     478            case VECTOR:
     479              doVector=action;
     480              break;
    477481            case MESSAGES:
    478482              models[iModel].messageHandler()->setPrefix(action!=0);
     
    590594              solveOptions.setSpecialOption(4,presolveOptions);
    591595              solveOptions.setSpecialOption(5,printOptions&1);
     596              if (doVector) {
     597                ClpMatrixBase * matrix = models[iModel].clpMatrix();
     598                if (dynamic_cast< ClpPackedMatrix*>(matrix)) {
     599                  ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
     600                  clpMatrix->makeSpecialColumnCopy();
     601                }
     602              }
    592603              if (method==ClpSolve::useDual) {
    593604                // dual
     
    663674              }
    664675              if (dualize) {
    665                 ((ClpSimplexOther *) models+iModel)->restoreFromDual(model2);
     676                int returnCode=((ClpSimplexOther *) models+iModel)->restoreFromDual(model2);
    666677                delete model2;
     678                if (returnCode&&dualize!=2)
     679                  models[iModel].primal(1);
    667680              }
    668681              if (status>=0)
     
    781794                fileName = "-";
    782795              } else {
     796                // See if .lp
     797                {
     798                  const char * c_name = field.c_str();
     799                  int length = strlen(c_name);
     800                  if (length>3&&!strncmp(c_name+length-3,".lp",3))
     801                    gmpl=-1; // .lp
     802                }
    783803                bool absolutePath;
    784804                if (dirsep=='/') {
     
    805825                } else {
    806826                  fileName = directory+field;
    807                   // See if gmpl (model & data)
     827                  // See if gmpl (model & data) - or even lp file
    808828                  int length = field.size();
    809829                  int percent = field.find('%');
     
    816836                    printf("GMPL model file %s and data file %s\n",
    817837                           fileName.c_str(),gmplData.c_str());
    818                   }
     838                  }
    819839                }
    820840                std::string name=fileName;
     
    842862                                                 keepImportNames!=0,
    843863                                                 allowImportErrors!=0);
    844                 else
     864                else if (gmpl>0)
    845865                  status= models[iModel].readGMPL(fileName.c_str(),
    846866                                                  (gmpl==2) ? gmplData.c_str() : NULL,
    847867                                                  keepImportNames!=0);
     868                else
     869                  status= models[iModel].readLp(fileName.c_str(),1.0e-12);
    848870                if (!status||(status>0&&allowImportErrors)) {
    849871                  goodModels[iModel]=true;
     
    937959                bool deleteModel2=false;
    938960                ClpSimplex * model2 = models+iModel;
     961                if (dualize) {
     962                  model2 = ((ClpSimplexOther *) model2)->dualOfModel();
     963                  printf("Dual of model has %d rows and %d columns\n",
     964                         model2->numberRows(),model2->numberColumns());
     965                  model2->setOptimizationDirection(1.0);
     966                  preSolve=0; // as picks up from model
     967                }
    939968                if (preSolve) {
    940969                  ClpPresolve pinfo;
     
    10311060                }
    10321061#else
    1033                 if (dualize) {
    1034                   ClpSimplex * model3 = ((ClpSimplexOther *) model2)->dualOfModel();
    1035                   printf("Dual of model has %d rows and %d columns\n",
    1036                          model3->numberRows(),model3->numberColumns());
    1037                   model3->writeMps(fileName.c_str(),(outputFormat-1)/2,1+((outputFormat-1)&1));
    1038                   delete model3;
    1039                 } else {
    1040                   model2->writeMps(fileName.c_str(),(outputFormat-1)/2,1+((outputFormat-1)&1));
    1041                 }
     1062                model2->writeMps(fileName.c_str(),(outputFormat-1)/2,1+((outputFormat-1)&1));
    10421063#endif
    10431064                if (deleteModel2)
     
    13791400              models[iModel].setSpecialOptions(0);
    13801401              mainTest(nFields,fields,algorithm,models[iModel],
    1381                        (preSolve!=0),specialOptions);
     1402                       (preSolve!=0),specialOptions,doVector!=0);
    13821403            }
    13831404            break;
     
    13981419              if (models[iModel].numberRows())
    13991420                algorithm=7;
    1400               mainTest(nFields,fields,algorithm,models[iModel],(preSolve!=0),specialOptions);
     1421              mainTest(nFields,fields,algorithm,models[iModel],(preSolve!=0),specialOptions,doVector!=0);
    14011422            }
    14021423            break;
  • trunk/Clp/src/ClpMatrixBase.cpp

    r800 r1034  
    548548// Correct sequence in and out to give true value
    549549void
    550 ClpMatrixBase::correctSequence(int & sequenceIn, int & sequenceOut) const
     550ClpMatrixBase::correctSequence(const ClpSimplex * model,int & sequenceIn, int & sequenceOut)
    551551{
    552552}
  • trunk/Clp/src/ClpMatrixBase.hpp

    r800 r1034  
    227227      mode=12  - after factorize but before permute stuff
    228228      mode=13  - at end of simplex to delete stuff
     229
    229230  */
    230231  virtual int generalExpanded(ClpSimplex * model,int mode,int & number);
     
    242243  /// Returns reduced cost of a variable
    243244  double reducedCost(ClpSimplex * model,int sequence) const;
    244   /// Correct sequence in and out to give true value
    245   virtual void correctSequence(int & sequenceIn, int & sequenceOut) const;
     245  /// Correct sequence in and out to give true value (if both -1 maybe do whole matrix)
     246  virtual void correctSequence(const ClpSimplex * model,int & sequenceIn, int & sequenceOut) ;
    246247  //@}
    247248 
  • trunk/Clp/src/ClpMessage.cpp

    r754 r1034  
    100100  {CLP_BAD_STRING_VALUES,3005,1,"%d string elements had no values associated with them"},
    101101  {CLP_CRUNCH_STATS,61,2,"Crunch %d (%d) rows, %d (%d) columns and %d (%d) elements"},
     102  {CLP_GENERAL,1000,1,"%s"},
    102103  {CLP_DUMMY_END,999999,0,""}
    103104};
     
    122123     addMessage(message->internalNumber,oneMessage);
    123124     message ++;
    124 }
     125  }
     126  // Put into compact form
     127  toCompact();
    125128
    126129  // now override any language ones
  • trunk/Clp/src/ClpMessage.hpp

    r754 r1034  
    9999  CLP_BAD_STRING_VALUES,
    100100  CLP_CRUNCH_STATS,
     101  CLP_GENERAL,
    101102  CLP_DUMMY_END
    102103};
  • trunk/Clp/src/ClpModel.cpp

    r1015 r1034  
    4141
    4242//#############################################################################
    43 
    44 ClpModel::ClpModel () :
     43ClpModel::ClpModel (bool emptyMessages) :
    4544
    4645  optimizationDirection_(1),
     
    7776  lengthNames_(0),
    7877  numberThreads_(0),
     78  specialOptions_(0),
    7979#ifndef CLP_NO_STD
    8080  defaultHandler_(true),
     
    103103  handler_->setLogLevel(1);
    104104  eventHandler_ = new ClpEventHandler();
    105   messages_ = ClpMessage();
    106   coinMessages_ = CoinMessage();
     105  if (!emptyMessages) {
     106    messages_ = ClpMessage();
     107    coinMessages_ = CoinMessage();
     108  }
    107109  CoinSeedRandom(1234567);
    108110}
     
    157159  eventHandler_=NULL;
    158160  whatsChanged_=0;
     161  specialOptions_ = 0;
    159162}
    160163//#############################################################################
     
    288291  gutsOfLoadModel(numrows, numcols,
    289292                  collb, colub, obj, rowlb, rowub, rowObjective);
    290   CoinPackedMatrix matrix(true,numrows,numcols,start[numcols],
     293  int numberElements = start ? start[numcols] : 0;
     294  CoinPackedMatrix matrix(true,numrows,numrows ? numcols : 0,numberElements,
    291295                              value,index,start,NULL);
    292296  matrix_ = new ClpPackedMatrix(matrix);
     
    319323ClpModel::loadProblem (  CoinModel & modelObject,bool tryPlusMinusOne)
    320324{
    321   if (modelObject.numberElements()==0)
     325  if (modelObject.numberColumns()==0&&modelObject.numberRows()==0)
    322326    return 0;
    323327  int numberErrors = 0;
     
    675679  userPointer_ = rhs.userPointer_;
    676680  scalingFlag_ = rhs.scalingFlag_;
     681  specialOptions_ = rhs.specialOptions_;
    677682  if (trueCopy) {
    678683#ifndef CLP_NO_STD
     
    707712    rowUpper_ = ClpCopyOfArray ( rhs.rowUpper_, numberRows_ );
    708713    columnLower_ = ClpCopyOfArray ( rhs.columnLower_, numberColumns_ );
     714    int scaleLength = ((specialOptions_&131072)==0) ? 1 : 2;
    709715    columnUpper_ = ClpCopyOfArray ( rhs.columnUpper_, numberColumns_ );
    710     rowScale_ = ClpCopyOfArray(rhs.rowScale_,numberRows_);
    711     columnScale_ = ClpCopyOfArray(rhs.columnScale_,numberColumns_);
     716    rowScale_ = ClpCopyOfArray(rhs.rowScale_,numberRows_*scaleLength);
     717    columnScale_ = ClpCopyOfArray(rhs.columnScale_,numberColumns_*scaleLength);
    712718    if (rhs.objective_)
    713719      objective_  = rhs.objective_->clone();
     
    830836      return false;
    831837    break;
    832   case ClpLastIntParam:
     838  default:
    833839    return false;
    834840  }
     
    875881    break;
    876882   
    877   case ClpLastDblParam:
     883  default:
    878884    return false;
    879885  }
     
    893899    break;
    894900
    895   case ClpLastStrParam:
     901  default:
    896902    return false;
    897903  }
     
    986992  CoinPackedMatrix matrix2;
    987993  matrix_=new ClpPackedMatrix(matrix2);
     994}
     995/* Really clean up matrix.
     996   a) eliminate all duplicate AND small elements in matrix
     997   b) remove all gaps and set extraGap_ and extraMajor_ to 0.0
     998   c) reallocate arrays and make max lengths equal to lengths
     999   d) orders elements
     1000   returns number of elements eliminated or -1 if not ClpMatrix
     1001*/
     1002int
     1003ClpModel::cleanMatrix(double threshold)
     1004{
     1005  ClpPackedMatrix * matrix = (dynamic_cast< ClpPackedMatrix*>(matrix_));
     1006  if (matrix) {
     1007    return matrix->getPackedMatrix()->cleanMatrix(threshold);
     1008  } else {
     1009    return -1;
     1010  }
    9881011}
    9891012// Resizes
     
    13051328    }
    13061329#endif
    1307     if (elements)
    1308       matrix_->appendMatrix(number,0,rowStarts,columns,elements);
     1330    //if (elements)
     1331    matrix_->appendMatrix(number,0,rowStarts,columns,elements);
    13091332  }
    13101333}
     
    17861809    }
    17871810#endif
    1788     if (elements)
    1789       matrix_->appendMatrix(number,1,columnStarts,rows,elements);
     1811    //if (elements)
     1812    matrix_->appendMatrix(number,1,columnStarts,rows,elements);
    17901813  }
    17911814}
     
    22582281{
    22592282  double * array = NULL;
    2260   if (problemStatus_==1&&!secondaryStatus_)
     2283  if (problemStatus_==1&&!secondaryStatus_) {
    22612284    array = ClpCopyOfArray(ray_,numberRows_);
     2285#if 0
     2286    // clean up
     2287    double largest=1.0e-30;
     2288    double smallest=COIN_DBL_MAX;
     2289    int i;
     2290    for (i=0;i<numberRows_;i++) {
     2291      double value = fabs(array[i]);
     2292      smallest = CoinMin(smallest,value);
     2293      largest = CoinMax(largest,value);
     2294    }
     2295#endif
     2296  }
    22622297  return array;
    22632298}
     
    23622397  }
    23632398  m.messageHandler()->setPrefix(savePrefix);
    2364   if (!status||ignoreErrors) {
     2399  if (!status||(ignoreErrors&&(status>0&&status<100000))) {
    23652400    loadProblem(*m.getMatrixByCol(),
    23662401                m.getColLower(),m.getColUpper(),
     
    27652800  strParam_[ClpProbName] = rhs->strParam_[ClpProbName];
    27662801#endif
    2767 
     2802  specialOptions_ = rhs->specialOptions_;
    27682803  optimizationDirection_ = rhs->optimizationDirection_;
    27692804  objectiveValue_=rhs->objectiveValue_;
     
    30083043  int iRow;
    30093044  for (iRow=first; iRow<last;iRow++) {
    3010     rowNames_[iRow]= rowNames[iRow-first];
    3011     maxLength = CoinMax(maxLength,(unsigned int) strlen(rowNames[iRow-first]));
     3045    if (rowNames[iRow-first]&&strlen(rowNames[iRow-first])) {
     3046      rowNames_[iRow]= rowNames[iRow-first];
     3047      maxLength = CoinMax(maxLength,(unsigned int) strlen(rowNames[iRow-first]));
     3048    } else {
     3049      maxLength = CoinMax(maxLength,(unsigned int) 8);
     3050      char name[9];
     3051      sprintf(name,"R%7.7d",iRow);
     3052      rowNames_[iRow]=name;
     3053    }
    30123054  }
    30133055  // May be too big - but we would have to check both rows and columns to be exact
     
    30243066  int iColumn;
    30253067  for (iColumn=first; iColumn<last;iColumn++) {
    3026     columnNames_[iColumn]= columnNames[iColumn-first];
    3027     maxLength = CoinMax(maxLength,(unsigned int) strlen(columnNames[iColumn-first]));
     3068    if (columnNames[iColumn-first]&&strlen(columnNames[iColumn-first])) {
     3069      columnNames_[iColumn]= columnNames[iColumn-first];
     3070      maxLength = CoinMax(maxLength,(unsigned int) strlen(columnNames[iColumn-first]));
     3071    } else {
     3072      maxLength = CoinMax(maxLength,(unsigned int) 8);
     3073      char name[9];
     3074      sprintf(name,"C%7.7d",iColumn);
     3075      columnNames_[iColumn]=name;
     3076    }
    30283077  }
    30293078  // May be too big - but we would have to check both rows and columns to be exact
     
    34163465  columnScale_ = NULL;
    34173466}
     3467void
     3468ClpModel::setSpecialOptions(unsigned int value)
     3469{
     3470  specialOptions_=value;
     3471}
     3472/* This creates a coinModel object
     3473 */
     3474CoinModel *
     3475ClpModel::createCoinModel() const
     3476{
     3477  CoinModel * coinModel = new CoinModel();
     3478  CoinPackedMatrix matrixByRow;
     3479  matrixByRow.reverseOrderedCopyOf(*matrix());
     3480  coinModel->setObjectiveOffset(objectiveOffset());
     3481  coinModel->setProblemName(problemName().c_str());
     3482
     3483  // Build by row from scratch
     3484  const double * element = matrixByRow.getElements();
     3485  const int * column = matrixByRow.getIndices();
     3486  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
     3487  const int * rowLength = matrixByRow.getVectorLengths();
     3488  int i;
     3489  for (i=0;i<numberRows_;i++) {
     3490    coinModel->addRow(rowLength[i],column+rowStart[i],
     3491                      element+rowStart[i],rowLower_[i],rowUpper_[i]);
     3492  }
     3493  // Now do column part
     3494  const double * objective = this->objective();
     3495  for (i=0;i<numberColumns_;i++) {
     3496    coinModel->setColumnBounds(i,columnLower_[i],columnUpper_[i]);
     3497    coinModel->setColumnObjective(i,objective[i]);
     3498  }
     3499  for ( i=0;i<numberColumns_;i++) {
     3500    if (isInteger(i))
     3501      coinModel->setColumnIsInteger(i,true);
     3502  }
     3503  // do names
     3504  for (i=0;i<numberRows_;i++) {
     3505    char temp[30];
     3506    strcpy(temp,rowName(i).c_str());
     3507    int length = strlen(temp);
     3508    for (int j=0;j<length;j++) {
     3509      if (temp[j]=='-')
     3510        temp[j]='_';
     3511    }
     3512    coinModel->setRowName(i,temp);
     3513  }
     3514  for (i=0;i<numberColumns_;i++) {
     3515    char temp[30];
     3516    strcpy(temp,columnName(i).c_str());
     3517    int length = strlen(temp);
     3518    for (int j=0;j<length;j++) {
     3519      if (temp[j]=='-')
     3520        temp[j]='_';
     3521    }
     3522    coinModel->setColumnName(i,temp);
     3523  }
     3524  ClpQuadraticObjective * obj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
     3525  if (obj) {
     3526    const CoinPackedMatrix * quadObj = obj->quadraticObjective();
     3527    // add in quadratic
     3528    const double * element = quadObj->getElements();
     3529    const int * row = quadObj->getIndices();
     3530    const CoinBigIndex * columnStart = quadObj->getVectorStarts();
     3531    const int * columnLength = quadObj->getVectorLengths();
     3532    for (i=0;i<numberColumns_;i++) {
     3533      int nels = columnLength[i];
     3534      if (nels) {
     3535        CoinBigIndex start = columnStart[i];
     3536        double constant = coinModel->getColumnObjective(i);
     3537        char temp[100000];
     3538        char temp2[30];
     3539        sprintf(temp,"%g",constant);
     3540        for (CoinBigIndex k=start;k<start+nels;k++) {
     3541          int kColumn = row[k];
     3542          double value = element[k];
     3543#if 1
     3544          // ampl gives twice with assumed 0.5
     3545          if (kColumn<i)
     3546            continue;
     3547          else if (kColumn==i)
     3548            value *= 0.5;
     3549#endif
     3550          if (value==1.0)
     3551            sprintf(temp2,"+%s",coinModel->getColumnName(kColumn));
     3552          else if (value==-1.0)
     3553            sprintf(temp2,"-%s",coinModel->getColumnName(kColumn));
     3554          else if (value>0.0)
     3555            sprintf(temp2,"+%g*%s",value,coinModel->getColumnName(kColumn));
     3556          else
     3557            sprintf(temp2,"%g*%s",value,coinModel->getColumnName(kColumn));
     3558          strcat(temp,temp2);
     3559          assert (strlen(temp)<100000);
     3560        }
     3561        coinModel->setObjective(i,temp);
     3562        if (logLevel()>2)
     3563          printf("el for objective column %s is %s\n",coinModel->getColumnName(i),temp);
     3564      }
     3565    }
     3566  }
     3567  return coinModel;
     3568}
    34183569//#############################################################################
    34193570// Constructors / Destructor / Assignment
  • trunk/Clp/src/ClpModel.hpp

    r800 r1034  
    4646  //@{
    4747    /// Default constructor
    48     ClpModel (  );
     48    ClpModel (bool emptyMessages=false  );
    4949
    5050  /** Copy constructor. May scale depending on mode
     
    249249  /// Create empty ClpPackedMatrix
    250250  void createEmptyMatrix();
     251  /** Really clean up matrix (if ClpPackedMatrix).
     252      a) eliminate all duplicate AND small elements in matrix
     253      b) remove all gaps and set extraGap_ and extraMajor_ to 0.0
     254      c) reallocate arrays and make max lengths equal to lengths
     255      d) orders elements
     256      returns number of elements eliminated or -1 if not ClpPackedMatrix
     257  */
     258  int cleanMatrix(double threshold=1.0e-20);
    251259#ifndef CLP_NO_STD
    252260  /// Drops names - makes lengthnames 0 and names empty
     
    268276  void setColumnName(int colIndex, std::string & name) ;
    269277#endif
     278  /** This creates a coinModel object
     279  */
     280  CoinModel * createCoinModel() const;
    270281
    271282    /** Write the problem in MPS format to the specified file.
     
    361372       4 - scaled problem optimal - unscaled problem has primal and dual infeasibilities
    362373       5 - giving up in primal with flagged variables
    363        6 - failed due to empty problem check
     374       6 - failed due to empty problem check 
    364375       7 - postSolve says not optimal
     376       8 - failed due to bad element check
    365377       100 up - translation of enum from ClpEventHandler
    366378   */
     
    745757    // Get an integer parameter
    746758    inline bool getIntParam(ClpIntParam key, int& value) const {
    747       if (key!=ClpLastIntParam) {
     759      if (key<ClpLastIntParam) {
    748760        value = intParam_[key];
    749761        return true;
     
    754766    // Get an double parameter
    755767    inline bool getDblParam(ClpDblParam key, double& value) const {
    756       if (key!=ClpLastDblParam) {
     768      if (key<ClpLastDblParam) {
    757769        value = dblParam_[key];
    758770        return true;
     
    764776    // Get a string parameter
    765777    inline bool getStrParam(ClpStrParam key, std::string& value) const {
    766       if (key!=ClpLastStrParam) {
     778      if (key<ClpLastStrParam) {
    767779        value = strParam_[key];
    768780        return true;
     
    774786    /// Create C++ lines to get to current state
    775787    void generateCpp( FILE * fp);
     788  /** For advanced options
     789      1 - Don't keep changing infeasibility weight
     790      2 - Keep nonLinearCost round solves
     791      4 - Force outgoing variables to exact bound (primal)
     792      8 - Safe to use dense initial factorization
     793      16 -Just use basic variables for operation if column generation
     794      32 -Clean up with primal before strong branching
     795      64 -Treat problem as feasible until last minute (i.e. minimize infeasibilities)
     796      128 - Switch off all matrix sanity checks
     797      256 - No row copy
     798      512 - If not in values pass, solution guaranteed, skip as much as possible
     799      1024 - In branch and bound
     800      2048 - Don't bother to re-factorize if < 20 iterations
     801      4096 - Skip some optimality checks
     802      8192 - Do Primal when cleaning up primal
     803      16384 - In fast dual (so we can switch off things)
     804      32678 - called from Osi
     805      65356 - keep arrays around as much as possible
     806      131072 - scale factor arrays have inverse values at end
     807      NOTE - many applications can call Clp but there may be some short cuts
     808             which are taken which are not guaranteed safe from all applications.
     809             Vetted applications will have a bit set and the code may test this
     810             At present I expect a few such applications - if too many I will
     811             have to re-think.  It is up to application owner to change the code
     812             if she/he needs these short cuts.  I will not debug unless in Coin
     813             repository.  See COIN_CLP_VETTED comments.
     814      0x01000000 is Cbc (and in branch and bound)
     815      0x02000000 is in a different branch and bound
     816  */
     817#define COIN_CBC_USING_CLP 0x01000000
     818  inline unsigned int specialOptions() const
     819  { return specialOptions_;};
     820  void setSpecialOptions(unsigned int value);
    776821  //@}
    777822
     
    902947  /// Number of threads (not very operational)
    903948  int numberThreads_;
     949  /** For advanced options
     950      See get and set for meaning
     951  */
     952  unsigned int specialOptions_;
    904953  /// Message handler
    905954  CoinMessageHandler * handler_;
  • trunk/Clp/src/ClpNetworkMatrix.cpp

    r930 r1034  
    567567      int iRowP = indices_[j+1];
    568568      indexRowU[numberElements]=iRowM;
    569       rowCount[iRowM++];
     569      rowCount[iRowM]++;
    570570      elementU[numberElements]=-1.0;
    571571      indexRowU[numberElements+1]=iRowP;
    572       rowCount[iRowP++];
     572      rowCount[iRowP]++;
    573573      elementU[numberElements+1]=1.0;
    574574      numberElements+=2;
     
    584584      if (iRowM>=0) {
    585585        indexRowU[numberElements]=iRowM;
    586         rowCount[iRowM++];
     586        rowCount[iRowM]++;
    587587        elementU[numberElements++]=-1.0;
    588588      }
    589589      if (iRowP>=0) {
    590590        indexRowU[numberElements]=iRowP;
    591         rowCount[iRowP++];
     591        rowCount[iRowP]++;
    592592        elementU[numberElements++]=1.0;
    593593      }
  • trunk/Clp/src/ClpNonLinearCost.cpp

    r892 r1034  
    7979  feasibleCost_=0.0;
    8080  infeasibilityWeight_ = -1.0;
     81  double * cost = model_->costRegion();
     82  // check if all 0
     83  int iSequence;
     84  bool allZero=true;
     85  for (iSequence=0;iSequence<numberTotal1;iSequence++) {
     86    if (cost[iSequence]) {
     87      allZero=false;
     88      break;
     89    }
     90  }
     91  if (allZero)
     92    model_->setInfeasibilityCost(1.0);
    8193  double infeasibilityCost = model_->infeasibilityCost();
    8294  sumInfeasibilities_=0.0;
     
    94106  infeasible_=NULL;
    95107
    96   int iSequence;
    97108  double * upper = model_->upperRegion();
    98109  double * lower = model_->lowerRegion();
    99   double * cost = model_->costRegion();
    100110
    101111  // See how we are storing things
     
    609619            // possibly below
    610620            lowerValue = lower_[iRange+1];
    611             if (value<lowerValue-primalTolerance) {
    612               value = lowerValue-value;
     621            if (value-lowerValue<-primalTolerance) {
     622              value = lowerValue-value-primalTolerance;
    613623#ifndef NDEBUG
    614624              if(value>1.0e15)
     
    626636            // possibly above
    627637            upperValue = lower_[iRange];
    628             if (value>upperValue+primalTolerance) {
    629               value = value-upperValue;
     638            if (value-upperValue>primalTolerance) {
     639              value = value-upperValue-primalTolerance;
    630640#ifndef NDEBUG
    631641              if(value>1.0e15)
     
    805815      case ClpSimplex::basic:
    806816      case ClpSimplex::superBasic:
    807         if (value<upperValue+primalTolerance) {
    808           if (value>lowerValue-primalTolerance) {
     817        if (value-upperValue<=primalTolerance) {
     818          if (value-lowerValue>=-primalTolerance) {
    809819            // feasible
    810820            //newWhere=CLP_FEASIBLE;
     
    812822            // below
    813823            newWhere=CLP_BELOW_LOWER;
    814             double infeasibility = lowerValue-value;
     824            double infeasibility = lowerValue-value-primalTolerance;
    815825            sumInfeasibilities_ += infeasibility;
    816826            largestInfeasibility_ = CoinMax(largestInfeasibility_,infeasibility);
     
    822832          // above
    823833          newWhere = CLP_ABOVE_UPPER;
    824           double infeasibility = value-upperValue;
     834          double infeasibility = value-upperValue-primalTolerance;
    825835          sumInfeasibilities_ += infeasibility;
    826836          largestInfeasibility_ = CoinMax(largestInfeasibility_,infeasibility);
     
    12651275      // get correct place
    12661276      int newWhere=CLP_FEASIBLE;
    1267       if (value<upperValue+primalTolerance) {
    1268         if (value>lowerValue-primalTolerance) {
     1277      if (value-upperValue<=primalTolerance) {
     1278        if (value-lowerValue>=-primalTolerance) {
    12691279          // feasible
    12701280          //newWhere=CLP_FEASIBLE;
     
    13781388      // get correct place
    13791389      int newWhere=CLP_FEASIBLE;
    1380       if (value<upperValue+primalTolerance) {
    1381         if (value>lowerValue-primalTolerance) {
     1390      if (value-upperValue<=primalTolerance) {
     1391        if (value-lowerValue>=-primalTolerance) {
    13821392          // feasible
    13831393          //newWhere=CLP_FEASIBLE;
     
    15211531    // get correct place
    15221532    int newWhere=CLP_FEASIBLE;
    1523     if (value<upperValue+primalTolerance) {
    1524       if (value>lowerValue-primalTolerance) {
     1533    if (value-upperValue<=primalTolerance) {
     1534      if (value-lowerValue>=-primalTolerance) {
    15251535        // feasible
    15261536        //newWhere=CLP_FEASIBLE;
     
    16001610    cost_[start+2] = costValue+infeasibilityCost;
    16011611    double primalTolerance = model_->currentPrimalTolerance();
    1602     if (solutionValue>=lowerValue-primalTolerance) {
    1603       if (solutionValue<=upperValue+primalTolerance) {
     1612    if (solutionValue-lowerValue>=-primalTolerance) {
     1613      if (solutionValue-upperValue<=primalTolerance) {
    16041614        iRange=start+1;
    16051615      } else {
     
    17341744      value=lowerValue;
    17351745    int newWhere=CLP_FEASIBLE;
    1736     if (value<upperValue+primalTolerance) {
    1737       if (value>lowerValue-primalTolerance) {
     1746    if (value-upperValue<=primalTolerance) {
     1747      if (value-lowerValue>=-primalTolerance) {
    17381748        // feasible
    17391749        //newWhere=CLP_FEASIBLE;
  • trunk/Clp/src/ClpObjective.hpp

    r754 r1034  
    4848                            double & predictedObj,
    4949                            double & thetaObj)=0;
     50  /// Return objective value (without any ClpModel offset) (model may be NULL)
     51  virtual double objectiveValue(const ClpSimplex * model, const double * solution) const = 0;
    5052  /// Resize objective
    5153  virtual void resize(int newNumberColumns) = 0;
     
    5860   */
    5961  virtual int markNonlinear(char * which);
     62  /// Say we have new primal solution - so may need to recompute
     63  virtual void newXValues() {};
    6064  //@}
    6165 
  • trunk/Clp/src/ClpPackedMatrix.cpp

    r1002 r1034  
    88#include "CoinIndexedVector.hpp"
    99#include "CoinHelperFunctions.hpp"
     10//#define THREAD
    1011
    1112#include "ClpSimplex.hpp"
     
    1516#include "ClpQuadraticObjective.hpp"
    1617#endif
    17 //#define THREAD
    1818// at end to get min/max!
    1919#include "ClpPackedMatrix.hpp"
     
    4444    matrix_(NULL),
    4545    numberActiveColumns_(0),
    46     zeroElements_(false),
    47     hasGaps_(true),
    48     rowCopy_(NULL)
     46    flags_(2),
     47    rowCopy_(NULL),
     48    columnCopy_(NULL)
    4949{
    5050  setType(1);
     
    5959  matrix_ = new CoinPackedMatrix(*(rhs.matrix_));
    6060  numberActiveColumns_ = rhs.numberActiveColumns_;
    61   zeroElements_ = rhs.zeroElements_;
    62   hasGaps_ = rhs.hasGaps_;
     61  flags_ = rhs.flags_;
    6362  int numberRows = getNumRows();
    6463  if (rhs.rhsOffset_&&numberRows) {
     
    6867  }
    6968  if (rhs.rowCopy_) {
     69    assert ((flags_&4)!=0);
    7070    rowCopy_ = new ClpPackedMatrix2(*rhs.rowCopy_);
    7171  } else {
    7272    rowCopy_ = NULL;
     73  }
     74  if (rhs.columnCopy_) {
     75    assert ((flags_&8)!=0);
     76    columnCopy_ = new ClpPackedMatrix3(*rhs.columnCopy_);
     77  } else {
     78    columnCopy_ = NULL;
    7379  }
    7480}
     
    8187
    8288  matrix_ = rhs;
    83   zeroElements_ = false;
    84   hasGaps_ = true;
     89  flags_ = 2;
    8590  numberActiveColumns_ = matrix_->getNumCols();
    8691  rowCopy_ = NULL;
     92  columnCopy_=NULL;
    8793  setType(1);
    8894 
     
    95101  numberActiveColumns_ = matrix_->getNumCols();
    96102  rowCopy_ = NULL;
    97   zeroElements_ = false;
    98   hasGaps_ = true;
     103  flags_ = 2;
     104  columnCopy_=NULL;
    99105  setType(1);
    100106 
     
    108114  delete matrix_;
    109115  delete rowCopy_;
     116  delete columnCopy_;
    110117}
    111118
     
    121128    matrix_ = new CoinPackedMatrix(*(rhs.matrix_));
    122129    numberActiveColumns_ = rhs.numberActiveColumns_;
    123     zeroElements_ = rhs.zeroElements_;
    124     hasGaps_ = rhs.hasGaps_;
     130    flags_ = rhs.flags_;
    125131    delete rowCopy_;
     132    delete columnCopy_;
    126133    if (rhs.rowCopy_) {
     134      assert ((flags_&4)!=0);
    127135      rowCopy_ = new ClpPackedMatrix2(*rhs.rowCopy_);
    128136    } else {
    129137      rowCopy_ = NULL;
     138    }
     139    if (rhs.columnCopy_) {
     140      assert ((flags_&8)!=0);
     141      columnCopy_ = new ClpPackedMatrix3(*rhs.columnCopy_);
     142    } else {
     143      columnCopy_ = NULL;
    130144    }
    131145  }
     
    160174                                 numberColumns,whichColumns);
    161175  numberActiveColumns_ = matrix_->getNumCols();
    162   zeroElements_ = rhs.zeroElements_;
    163   hasGaps_ = rhs.hasGaps_;
    164176  rowCopy_ = NULL;
     177  flags_ = rhs.flags_;
     178  columnCopy_=NULL;
    165179}
    166180ClpPackedMatrix::ClpPackedMatrix (
     
    173187                                 numberColumns,whichColumns);
    174188  numberActiveColumns_ = matrix_->getNumCols();
    175   zeroElements_ = false;
    176   hasGaps_=true;
    177189  rowCopy_ = NULL;
     190  flags_ = 2;
     191  columnCopy_=NULL;
    178192  setType(1);
    179193}
     
    190204  //copy->matrix_->removeGaps();
    191205  copy->numberActiveColumns_ = copy->matrix_->getNumCols();
    192   copy->hasGaps_=false;
     206  copy->flags_ = flags_&(~2); // no gaps
    193207  return copy;
    194208}
     
    227241  const int * columnLength = matrix_->getVectorLengths();
    228242  const double * elementByColumn = matrix_->getElements();
    229   if (!hasGaps_) {
     243  if (!(flags_&2)) {
    230244    if (scalar==1.0) {
    231245      CoinBigIndex start=columnStart[0];
     
    326340    const double * elementByColumn = matrix_->getElements();
    327341    if (!spare) {
    328       if (!hasGaps_) {
     342      if (!(flags_&2)) {
    329343        CoinBigIndex start=columnStart[0];
    330344        for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
     
    359373      for (iRow=0;iRow<numberRows;iRow++)
    360374        spare[iRow] = x[iRow]*rowScale[iRow];
    361       if (!hasGaps_) {
     375      if (!(flags_&2)) {
    362376        CoinBigIndex start=columnStart[0];
    363377        for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
     
    419433#endif
    420434  bool packed = rowArray->packedMode();
    421   double factor = 0.27;
     435  double factor = 0.30;
    422436  // We may not want to do by row if there may be cache problems
    423437  // It would be nice to find L2 cache size - for moment 512K
     
    440454  double factor2 = factor*multiplierX;
    441455  if (packed&&rowCopy_&&numberInRowArray>2&&numberInRowArray>factor2*numberRows&&
    442       numberInRowArray<0.9*numberRows) {
     456      numberInRowArray<0.9*numberRows&&0) {
    443457    rowCopy_->transposeTimes(model,rowCopy->getPackedMatrix(),rowArray,y,columnArray);
    444458    return;
     
    447461    // do by column
    448462    // If no gaps - can do a bit faster
    449     if (!hasGaps_) {
     463    if (!(flags_&2)||columnCopy_) {
    450464      transposeTimesByColumn( model,  scalar,
    451465                              rowArray, y, columnArray);
     
    668682        }
    669683      }
     684#if 0
    670685      double value = 0.0;
    671686      CoinBigIndex j;
     
    692707        index[numberNonZero++]=iColumn;
    693708      }
     709#else
     710      if (!columnCopy_)
     711        gutsOfTransposeTimesUnscaled(pi,columnArray,zeroTolerance);
     712      else
     713        columnCopy_->transposeTimes(model,pi,columnArray);
     714      numberNonZero = columnArray->getNumElements();
     715#endif
    694716    } else {
    695717      // scaled
     
    707729      }
    708730      const double * columnScale = model->columnScale();
     731#if 0
    709732      double value = 0.0;
    710733      double scale=columnScale[0];
     
    725748        }
    726749        value = 0.0;
    727         for (j=start; j<end;j++) {
    728           int iRow = row[j];
    729           value += pi[iRow]*elementByColumn[j];
     750        if (model->getColumnStatus(iColumn+1)!=ClpSimplex::basic) {
     751          for (j=start; j<end;j++) {
     752            int iRow = row[j];
     753            value += pi[iRow]*elementByColumn[j];
     754          }
    730755        }
    731756      }
     
    735760        index[numberNonZero++]=iColumn;
    736761      }
     762#else
     763      if (!columnCopy_)
     764        gutsOfTransposeTimesScaled(pi,columnScale,columnArray,zeroTolerance);
     765      else
     766        columnCopy_->transposeTimes(model,pi,columnArray);
     767      numberNonZero = columnArray->getNumElements();
     768#endif
    737769    }
    738770    // zero out
     
    899931    int numberOriginal = 0;
    900932    if (packed) {
     933#if 0
    901934      numberNonZero=0;
    902935      // and set up mark as char array
     
    937970        }
    938971      }
     972#else
     973      gutsOfTransposeTimesByRowGE3(rowArray,columnArray,y,zeroTolerance,scalar);
     974      numberNonZero = columnArray->getNumElements();
     975#endif
    939976    } else {
    940977      double * markVector = y->denseVector();
     
    9821019    double value;
    9831020    if (packed) {
     1021#if 0
    9841022      int iRow0 = whichRow[0];
    9851023      int iRow1 = whichRow[1];
     
    10431081        }
    10441082      }
     1083#else
     1084      gutsOfTransposeTimesByRowEQ2(rowArray,columnArray,y,zeroTolerance,scalar);
     1085      numberNonZero = columnArray->getNumElements();
     1086#endif
    10451087    } else {
    10461088      int iRow = whichRow[0];
     
    10821124    CoinBigIndex j;
    10831125    if (packed) {
     1126#if 0
    10841127      double value = pi[0]*scalar;
    10851128      for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
     
    10911134        }
    10921135      }
     1136#else
     1137      gutsOfTransposeTimesByRowEQ1(rowArray,columnArray,zeroTolerance,scalar);
     1138      numberNonZero = columnArray->getNumElements();
     1139#endif
    10931140    } else {
    10941141      double value = pi[iRow]*scalar;
     
    11051152  columnArray->setNumElements(numberNonZero);
    11061153  y->setNumElements(0);
     1154}
     1155// Meat of transposeTimes by column when not scaled
     1156void
     1157ClpPackedMatrix::gutsOfTransposeTimesUnscaled(const double * pi,CoinIndexedVector * output, const double zeroTolerance) const
     1158{
     1159  int numberNonZero=0;
     1160  int * index = output->getIndices();
     1161  double * array = output->denseVector();
     1162  // get matrix data pointers
     1163  const int * row = matrix_->getIndices();
     1164  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
     1165  const double * elementByColumn = matrix_->getElements();
     1166  double value = 0.0;
     1167  CoinBigIndex j;
     1168  CoinBigIndex end = columnStart[1];
     1169  for (j=columnStart[0]; j<end;j++) {
     1170    int iRow = row[j];
     1171    value += pi[iRow]*elementByColumn[j];
     1172  }
     1173  int iColumn;
     1174  for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
     1175    CoinBigIndex start = end;
     1176    end = columnStart[iColumn+2];
     1177    if (fabs(value)>zeroTolerance) {
     1178      array[numberNonZero]=value;
     1179      index[numberNonZero++]=iColumn;
     1180    }
     1181    value = 0.0;
     1182    for (j=start; j<end;j++) {
     1183      int iRow = row[j];
     1184      value += pi[iRow]*elementByColumn[j];
     1185    }
     1186  }
     1187  if (fabs(value)>zeroTolerance) {
     1188    array[numberNonZero]=value;
     1189    index[numberNonZero++]=iColumn;
     1190  }
     1191  output->setNumElements(numberNonZero);
     1192}
     1193// Meat of transposeTimes by column when scaled
     1194void
     1195ClpPackedMatrix::gutsOfTransposeTimesScaled(const double * pi,const double * columnScale,
     1196                                           CoinIndexedVector * output, const double zeroTolerance) const
     1197{
     1198  int numberNonZero=0;
     1199  int * index = output->getIndices();
     1200  double * array = output->denseVector();
     1201  // get matrix data pointers
     1202  const int * row = matrix_->getIndices();
     1203  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
     1204  const double * elementByColumn = matrix_->getElements();
     1205  double value = 0.0;
     1206  double scale=columnScale[0];
     1207  CoinBigIndex j;
     1208  CoinBigIndex end = columnStart[1];
     1209  for (j=columnStart[0]; j<end;j++) {
     1210    int iRow = row[j];
     1211    value += pi[iRow]*elementByColumn[j];
     1212  }
     1213  int iColumn;
     1214  for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
     1215    value *= scale;
     1216    CoinBigIndex start = end;
     1217    scale = columnScale[iColumn+1];
     1218    end = columnStart[iColumn+2];
     1219    if (fabs(value)>zeroTolerance) {
     1220      array[numberNonZero]=value;
     1221      index[numberNonZero++]=iColumn;
     1222    }
     1223    value = 0.0;
     1224    for (j=start; j<end;j++) {
     1225      int iRow = row[j];
     1226      value += pi[iRow]*elementByColumn[j];
     1227    }
     1228  }
     1229  value *= scale;
     1230  if (fabs(value)>zeroTolerance) {
     1231    array[numberNonZero]=value;
     1232    index[numberNonZero++]=iColumn;
     1233  }
     1234  output->setNumElements(numberNonZero);
     1235}
     1236// Meat of transposeTimes by row n > 2 if packed
     1237void
     1238ClpPackedMatrix::gutsOfTransposeTimesByRowGE3(const CoinIndexedVector * piVector, CoinIndexedVector * output,
     1239                                             CoinIndexedVector * spareVector, const double tolerance, const double scalar) const
     1240{
     1241  double * pi = piVector->denseVector();
     1242  int numberNonZero=0;
     1243  int * index = output->getIndices();
     1244  double * array = output->denseVector();
     1245  int numberInRowArray = piVector->getNumElements();
     1246  const int * column = getIndices();
     1247  const CoinBigIndex * rowStart = getVectorStarts();
     1248  const double * element = getElements();
     1249  const int * whichRow = piVector->getIndices();
     1250  // ** Row copy is already scaled
     1251  int iRow;
     1252  int i;
     1253  // and set up mark as char array
     1254  char * marked = (char *) (index+output->capacity());
     1255  int * lookup = spareVector->getIndices();
     1256#if 1
     1257  for (i=0;i<numberInRowArray;i++) {
     1258    iRow = whichRow[i];
     1259    double value = pi[i]*scalar;
     1260    CoinBigIndex j;
     1261    for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
     1262      int iColumn = column[j];
     1263      double elValue = element[j];
     1264      if (!marked[iColumn]) {
     1265        marked[iColumn]=1;
     1266        lookup[iColumn]=numberNonZero;
     1267        array[numberNonZero] = value*elValue;
     1268        index[numberNonZero++]=iColumn;
     1269      } else {
     1270        int k = lookup[iColumn];
     1271        array[k] += value*elValue;
     1272      }
     1273    }
     1274  }
     1275#else
     1276  int nextRow = whichRow[0];
     1277  CoinBigIndex nextStart = rowStart[nextRow];
     1278  CoinBigIndex nextEnd = rowStart[nextRow+1];
     1279  whichRow[numberInRowArray]=0; // for electricfence etc
     1280  for (i=0;i<numberInRowArray;i++) {
     1281    iRow = nextRow;
     1282    nextRow = whichRow[i+1];
     1283    CoinBigIndex start = nextStart;
     1284    CoinBigIndex end = nextEnd;
     1285    double value = pi[i]*scalar;
     1286    CoinBigIndex j;
     1287    nextRow = whichRow[i+1];
     1288    nextStart = rowStart[nextRow];
     1289    nextEnd = rowStart[nextRow+1];
     1290    for (j=start;j<end;j++) {
     1291      int iColumn = column[j];
     1292      double elValue = element[j];
     1293      if (!marked[iColumn]) {
     1294        marked[iColumn]=1;
     1295        lookup[iColumn]=numberNonZero;
     1296        array[numberNonZero] = value*elValue;
     1297        index[numberNonZero++]=iColumn;
     1298      } else {
     1299        int k = lookup[iColumn];
     1300        array[k] += value*elValue;
     1301      }
     1302    }
     1303  }
     1304#endif
     1305  // get rid of tiny values and zero out marked
     1306  for (i=0;i<numberNonZero;i++) {
     1307    int iColumn = index[i];
     1308    marked[iColumn]=0;
     1309    double value = array[i];
     1310    while (fabs(value)<=tolerance) {
     1311      numberNonZero--;
     1312      value = array[numberNonZero];
     1313      iColumn = index[numberNonZero];
     1314      marked[iColumn]=0;
     1315      if (i<numberNonZero) {
     1316        array[numberNonZero]=0.0;
     1317        array[i] = value;
     1318        index[i] = iColumn;
     1319      } else {
     1320        array[i]=0.0;
     1321        value =1.0; // to force end of while
     1322      }
     1323    }
     1324  }
     1325  output->setNumElements(numberNonZero);
     1326  spareVector->setNumElements(0);
     1327}
     1328// Meat of transposeTimes by row n == 2 if packed
     1329void
     1330ClpPackedMatrix::gutsOfTransposeTimesByRowEQ2(const CoinIndexedVector * piVector, CoinIndexedVector * output,
     1331                                   CoinIndexedVector * spareVector, const double tolerance, const double scalar) const
     1332{
     1333  double * pi = piVector->denseVector();
     1334  int numberNonZero=0;
     1335  int * index = output->getIndices();
     1336  double * array = output->denseVector();
     1337  const int * column = getIndices();
     1338  const CoinBigIndex * rowStart = getVectorStarts();
     1339  const double * element = getElements();
     1340  const int * whichRow = piVector->getIndices();
     1341  int iRow0 = whichRow[0];
     1342  int iRow1 = whichRow[1];
     1343  double pi0 = pi[0];
     1344  double pi1 = pi[1];
     1345  if (rowStart[iRow0+1]-rowStart[iRow0]>
     1346      rowStart[iRow1+1]-rowStart[iRow1]) {
     1347    // do one with fewer first
     1348    iRow0=iRow1;
     1349    iRow1=whichRow[0];
     1350    pi0=pi1;
     1351    pi1=pi[0];
     1352  }
     1353  // and set up mark as char array
     1354  char * marked = (char *) (index+output->capacity());
     1355  int * lookup = spareVector->getIndices();
     1356  double value = pi0*scalar;
     1357  CoinBigIndex j;
     1358  for (j=rowStart[iRow0];j<rowStart[iRow0+1];j++) {
     1359    int iColumn = column[j];
     1360    double elValue = element[j];
     1361    double value2 = value*elValue;
     1362    array[numberNonZero] = value2;
     1363    marked[iColumn]=1;
     1364    lookup[iColumn]=numberNonZero;
     1365    index[numberNonZero++]=iColumn;
     1366  }
     1367  int numberOriginal = numberNonZero;
     1368  value = pi1*scalar;
     1369  for (j=rowStart[iRow1];j<rowStart[iRow1+1];j++) {
     1370    int iColumn = column[j];
     1371    double elValue = element[j];
     1372    double value2 = value*elValue;
     1373    // I am assuming no zeros in matrix
     1374    if (marked[iColumn]) {
     1375      int iLookup = lookup[iColumn];
     1376      array[iLookup] += value2;
     1377    } else {
     1378      if (fabs(value2)>tolerance) {
     1379        array[numberNonZero] = value2;
     1380        index[numberNonZero++]=iColumn;
     1381      }
     1382    }
     1383  }
     1384  // get rid of tiny values and zero out marked
     1385  int i;
     1386  int iFirst=numberNonZero;
     1387  for (i=0;i<numberOriginal;i++) {
     1388    int iColumn = index[i];
     1389    marked[iColumn]=0;
     1390    if (fabs(array[i])<=tolerance) {
     1391      if (numberNonZero>numberOriginal) {
     1392        numberNonZero--;
     1393        double value = array[numberNonZero];
     1394        array[numberNonZero]=0.0;
     1395        array[i]=value;
     1396        index[i]=index[numberNonZero];
     1397      } else {
     1398        iFirst = i;
     1399      }
     1400    }
     1401  }
     1402 
     1403  if (iFirst<numberNonZero) {
     1404    int n=iFirst;
     1405    for (i=n;i<numberOriginal;i++) {
     1406      int iColumn = index[i];
     1407      double value = array[i];
     1408      array[i]=0.0;
     1409      if (fabs(value)>tolerance) {
     1410        array[n]=value;
     1411        index[n++]=iColumn;
     1412      }
     1413    }
     1414    for (;i<numberNonZero;i++) {
     1415      int iColumn = index[i];
     1416      double value = array[i];
     1417      array[i]=0.0;
     1418      array[n]=value;
     1419      index[n++]=iColumn;
     1420    }
     1421    numberNonZero=n;
     1422  }
     1423  output->setNumElements(numberNonZero);
     1424  spareVector->setNumElements(0);
     1425}
     1426// Meat of transposeTimes by row n == 1 if packed
     1427void
     1428ClpPackedMatrix::gutsOfTransposeTimesByRowEQ1(const CoinIndexedVector * piVector, CoinIndexedVector * output,
     1429                                   const double tolerance, const double scalar) const
     1430{
     1431  double * pi = piVector->denseVector();
     1432  int numberNonZero=0;
     1433  int * index = output->getIndices();
     1434  double * array = output->denseVector();
     1435  const int * column = getIndices();
     1436  const CoinBigIndex * rowStart = getVectorStarts();
     1437  const double * element = getElements();
     1438  int iRow=piVector->getIndices()[0];
     1439  numberNonZero=0;
     1440  CoinBigIndex j;
     1441  double value = pi[0]*scalar;
     1442  for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
     1443    int iColumn = column[j];
     1444    double elValue = element[j];
     1445    double value2 = value*elValue;
     1446    if (fabs(value2)>tolerance) {
     1447      array[numberNonZero] = value2;
     1448      index[numberNonZero++]=iColumn;
     1449    }
     1450  }
     1451  output->setNumElements(numberNonZero);
    11071452}
    11081453/* Return <code>x *A in <code>z</code> but
     
    11291474  assert (!rowArray->packedMode());
    11301475  columnArray->setPacked();
    1131   if (!hasGaps_&&numberToDo>5) {
     1476  if (!(flags_&2)&&numberToDo>5) {
    11321477    // no gaps and a reasonable number
    11331478    if (!rowScale) {
     
    12221567  bool packed = pi->packedMode();
    12231568  // factor should be smaller if doing both with two pi vectors
    1224   double factor = 0.27;
     1569  double factor = 0.30;
    12251570  // We may not want to do by row if there may be cache problems
    12261571  // It would be nice to find L2 cache size - for moment 512K
     
    12391584  if (!packed)
    12401585    factor *= 0.9;
    1241   return ((numberInRowArray>factor*numberRows||!model->rowCopy())&&!hasGaps_);
     1586  return ((numberInRowArray>factor*numberRows||!model->rowCopy())&&!(flags_&2));
    12421587}
    12431588#ifndef CLP_ALL_ONE_FILE
     
    12911636        pi[iRow]=piOld[i];
    12921637      }
    1293       CoinBigIndex j;
    1294       CoinBigIndex end = columnStart[0];
    1295       for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
    1296         CoinBigIndex start = end;
    1297         end = columnStart[iColumn+1];
    1298         ClpSimplex::Status status = model->getStatus(iColumn);
    1299         if (status==ClpSimplex::basic||status==ClpSimplex::isFixed) continue;
    1300         double value = 0.0;
    1301         for (j=start; j<end;j++) {
    1302           int iRow = row[j];
    1303           value -= pi[iRow]*elementByColumn[j];
    1304         }
    1305         if (fabs(value)>zeroTolerance) {
    1306           // and do other array
    1307           double modification = 0.0;
    1308           for (j=start; j<end;j++) {
    1309             int iRow = row[j];
    1310             modification += piWeight[iRow]*elementByColumn[j];
    1311           }
    1312           double thisWeight = weights[iColumn];
    1313           double pivot = value*scaleFactor;
    1314           double pivotSquared = pivot * pivot;
    1315           thisWeight += pivotSquared * devex + pivot * modification;
    1316           if (thisWeight<DEVEX_TRY_NORM) {
    1317             if (referenceIn<0.0) {
    1318               // steepest
    1319               thisWeight = CoinMax(DEVEX_TRY_NORM,DEVEX_ADD_ONE+pivotSquared);
    1320             } else {
    1321               // exact
    1322               thisWeight = referenceIn*pivotSquared;
    1323               if (reference(iColumn))
    1324                 thisWeight += 1.0;
    1325               thisWeight = CoinMax(thisWeight,DEVEX_TRY_NORM);
    1326             }
    1327           }
    1328           weights[iColumn] = thisWeight;
    1329           if (!killDjs) {
    1330             array[numberNonZero]=value;
    1331             index[numberNonZero++]=iColumn;
    1332           }
    1333         }
     1638      if (!columnCopy_) {
     1639        CoinBigIndex j;
     1640        CoinBigIndex end = columnStart[0];
     1641        for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
     1642          CoinBigIndex start = end;
     1643          end = columnStart[iColumn+1];
     1644          ClpSimplex::Status status = model->getStatus(iColumn);
     1645          if (status==ClpSimplex::basic||status==ClpSimplex::isFixed) continue;
     1646          double value = 0.0;
     1647          for (j=start; j<end;j++) {
     1648            int iRow = row[j];
     1649            value -= pi[iRow]*elementByColumn[j];
     1650          }
     1651          if (fabs(value)>zeroTolerance) {
     1652            // and do other array
     1653            double modification = 0.0;
     1654            for (j=start; j<end;j++) {
     1655              int iRow = row[j];
     1656              modification += piWeight[iRow]*elementByColumn[j];
     1657            }
     1658            double thisWeight = weights[iColumn];
     1659            double pivot = value*scaleFactor;
     1660            double pivotSquared = pivot * pivot;
     1661            thisWeight += pivotSquared * devex + pivot * modification;
     1662            if (thisWeight<DEVEX_TRY_NORM) {
     1663              if (referenceIn<0.0) {
     1664                // steepest
     1665                thisWeight = CoinMax(DEVEX_TRY_NORM,DEVEX_ADD_ONE+pivotSquared);
     1666              } else {
     1667                // exact
     1668                thisWeight = referenceIn*pivotSquared;
     1669                if (reference(iColumn))
     1670                  thisWeight += 1.0;
     1671                thisWeight = CoinMax(thisWeight,DEVEX_TRY_NORM);
     1672              }
     1673            }
     1674            weights[iColumn] = thisWeight;
     1675            if (!killDjs) {
     1676              array[numberNonZero]=value;
     1677              index[numberNonZero++]=iColumn;
     1678            }
     1679          }
     1680        }
     1681      } else {
     1682        // use special column copy
     1683        // reset back
     1684        if (killDjs)
     1685          scaleFactor=0.0;
     1686        columnCopy_->transposeTimes2(model,pi,dj1,piWeight,referenceIn,devex,
     1687                                     reference,weights,scaleFactor);
     1688        numberNonZero = dj1->getNumElements();
    13341689      }
    13351690    } else {
     
    13401695        pi[iRow]=piOld[i]*rowScale[iRow];
    13411696      }
    1342       const double * columnScale = model->columnScale();
    1343       CoinBigIndex j;
    1344       CoinBigIndex end = columnStart[0];
    1345       for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
    1346         CoinBigIndex start = end;
    1347         end = columnStart[iColumn+1];
    1348         ClpSimplex::Status status = model->getStatus(iColumn);
    1349         if (status==ClpSimplex::basic||status==ClpSimplex::isFixed) continue;
    1350         double scale=columnScale[iColumn];
    1351         double value = 0.0;
    1352         for (j=start; j<end;j++) {
    1353           int iRow = row[j];
    1354           value -= pi[iRow]*elementByColumn[j];
    1355         }
    1356         value *= scale;
    1357         if (fabs(value)>zeroTolerance) {
    1358           double modification = 0.0;
    1359           for (j=start; j<end;j++) {
    1360             int iRow = row[j];
    1361             modification += piWeight[iRow]*elementByColumn[j]*rowScale[iRow];
    1362           }
    1363           modification *= scale;
    1364           double thisWeight = weights[iColumn];
    1365           double pivot = value*scaleFactor;
    1366           double pivotSquared = pivot * pivot;
    1367           thisWeight += pivotSquared * devex + pivot * modification;
    1368           if (thisWeight<DEVEX_TRY_NORM) {
    1369             if (referenceIn<0.0) {
    1370               // steepest
    1371               thisWeight = CoinMax(DEVEX_TRY_NORM,DEVEX_ADD_ONE+pivotSquared);
    1372             } else {
    1373               // exact
    1374               thisWeight = referenceIn*pivotSquared;
    1375               if (reference(iColumn))
    1376                 thisWeight += 1.0;
    1377               thisWeight = CoinMax(thisWeight,DEVEX_TRY_NORM);
    1378             }
    1379           }
    1380           weights[iColumn] = thisWeight;
    1381           if (!killDjs) {
    1382             array[numberNonZero]=value;
    1383             index[numberNonZero++]=iColumn;
    1384           }
    1385         }
     1697      // can also scale piWeight as not used again
     1698      int numberWeight = pi2->getNumElements();
     1699      const int * indexWeight = pi2->getIndices();
     1700      for (i=0;i<numberWeight;i++) {
     1701        int iRow = indexWeight[i];
     1702        piWeight[iRow] *= rowScale[iRow];
     1703      }
     1704      if (!columnCopy_) {
     1705        const double * columnScale = model->columnScale();
     1706        CoinBigIndex j;
     1707        CoinBigIndex end = columnStart[0];
     1708        for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
     1709          CoinBigIndex start = end;
     1710          end = columnStart[iColumn+1];
     1711          ClpSimplex::Status status = model->getStatus(iColumn);
     1712          if (status==ClpSimplex::basic||status==ClpSimplex::isFixed) continue;
     1713          double scale=columnScale[iColumn];
     1714          double value = 0.0;
     1715          for (j=start; j<end;j++) {
     1716            int iRow = row[j];
     1717            value -= pi[iRow]*elementByColumn[j];
     1718          }
     1719          value *= scale;
     1720          if (fabs(value)>zeroTolerance) {
     1721            double modification = 0.0;
     1722            for (j=start; j<end;j++) {
     1723              int iRow = row[j];
     1724              modification += piWeight[iRow]*elementByColumn[j];
     1725            }
     1726            modification *= scale;
     1727            double thisWeight = weights[iColumn];
     1728            double pivot = value*scaleFactor;
     1729            double pivotSquared = pivot * pivot;
     1730            thisWeight += pivotSquared * devex + pivot * modification;
     1731            if (thisWeight<DEVEX_TRY_NORM) {
     1732              if (referenceIn<0.0) {
     1733                // steepest
     1734                thisWeight = CoinMax(DEVEX_TRY_NORM,DEVEX_ADD_ONE+pivotSquared);
     1735              } else {
     1736                // exact
     1737                thisWeight = referenceIn*pivotSquared;
     1738                if (reference(iColumn))
     1739                  thisWeight += 1.0;
     1740                thisWeight = CoinMax(thisWeight,DEVEX_TRY_NORM);
     1741              }
     1742            }
     1743            weights[iColumn] = thisWeight;
     1744            if (!killDjs) {
     1745              array[numberNonZero]=value;
     1746              index[numberNonZero++]=iColumn;
     1747            }
     1748          }
     1749        }
     1750      } else {
     1751        // use special column copy
     1752        // reset back
     1753        if (killDjs)
     1754          scaleFactor=0.0;
     1755        columnCopy_->transposeTimes2(model,pi,dj1,piWeight,referenceIn,devex,
     1756                                     reference,weights,scaleFactor);
     1757        numberNonZero = dj1->getNumElements();
    13861758      }
    13871759    }
     
    14421814    } else {
    14431815      // scaled
     1816      // can also scale piWeight as not used again
     1817      int numberWeight = pi2->getNumElements();
     1818      const int * indexWeight = pi2->getIndices();
     1819      for (int i=0;i<numberWeight;i++) {
     1820        int iRow = indexWeight[i];
     1821        piWeight[iRow] *= rowScale[iRow];
     1822      }
    14441823      const double * columnScale = model->columnScale();
    14451824      CoinBigIndex j;
     
    14611840          for (j=start; j<end;j++) {
    14621841            int iRow = row[j];
    1463             modification += piWeight[iRow]*elementByColumn[j]*rowScale[iRow];
     1842            modification += piWeight[iRow]*elementByColumn[j];
    14641843          }
    14651844          modification *= scale;
     
    16181997  const int * row = matrix_->getIndices();
    16191998  const double * elementByColumn = matrix_->getElements();
    1620   if (!zeroElements_) {
     1999  if ((flags_&1)==0) {
    16212000    if (!rowScale) {
    16222001      // no scaling
     
    17002079ClpPackedMatrix::scale(ClpModel * model) const
    17012080{
     2081#ifndef NDEBUG
     2082  //checkFlags();
     2083#endif
    17022084  int numberRows = model->numberRows();
    17032085  int numberColumns = matrix_->getNumCols();
     
    17232105  const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
    17242106  const double * element = rowCopy->getElements();
    1725   double * rowScale = new double [numberRows];
    1726   double * columnScale = new double [numberColumns];
     2107  int scaleLength = ((model->specialOptions()&131072)==0) ? 1 : 2;
     2108  double * rowScale = new double [numberRows*scaleLength];
     2109  double * columnScale = new double [numberColumns*scaleLength];
    17272110  // we are going to mark bits we are interested in
    17282111  char * usefulRow = new char [numberRows];
     
    21082491    }
    21092492#endif
     2493    if (scaleLength>1) {
     2494      // make copy
     2495      double * inverseScale = rowScale + numberRows;
     2496      for (iRow=0;iRow<numberRows;iRow++)
     2497        inverseScale[iRow] = 1.0/rowScale[iRow] ;
     2498      inverseScale = columnScale + numberColumns;
     2499      for (iColumn=0;iColumn<numberColumns;iColumn++)
     2500        inverseScale[iColumn] = 1.0/columnScale[iColumn] ;
     2501    }
    21102502    model->setRowScale(rowScale);
    21112503    model->setColumnScale(columnScale);
     
    21772569  const double * elementByColumn = matrix_->getElements();
    21782570  CoinBigIndex i;
     2571  int * index = rowArray->getIndices();
     2572  double * array = rowArray->denseVector();
     2573  int number = 0;
    21792574  if (!rowScale) {
    2180     CoinBigIndex j=columnStart[iColumn];
    2181     rowArray->createPacked(columnLength[iColumn],
    2182                            row+j,elementByColumn+j);
     2575    for (i=columnStart[iColumn];
     2576         i<columnStart[iColumn]+columnLength[iColumn];i++) {
     2577      int iRow = row[i];
     2578      double value = elementByColumn[i];
     2579      if (value) {
     2580        array[number]=value;
     2581        index[number++]=iRow;
     2582      }
     2583    }
     2584    rowArray->setNumElements(number);
     2585    rowArray->setPackedMode(true);
    21832586  } else {
    21842587    // apply scaling
    21852588    double scale = model->columnScale()[iColumn];
    2186     int * index = rowArray->getIndices();
    2187     double * array = rowArray->denseVector();
    2188     int number = 0;
    21892589    for (i=columnStart[iColumn];
    21902590         i<columnStart[iColumn]+columnLength[iColumn];i++) {
    21912591      int iRow = row[i];
    2192       array[number]=elementByColumn[i]*scale*rowScale[iRow];
    2193       index[number++]=iRow;
     2592      double value = elementByColumn[i]*scale*rowScale[iRow];
     2593      if (value) {
     2594        array[number]=value;
     2595        index[number++]=iRow;
     2596      }
    21942597    }
    21952598    rowArray->setNumElements(number);
     
    22842687  int numberColumns = matrix_->getNumCols();
    22852688  // Say no gaps
    2286   hasGaps_=false;
    2287   if (check==14) {
     2689  flags_ &= ~2;
     2690  if (check==14||check==10) {
    22882691    for (iColumn=0;iColumn<numberColumns;iColumn++) {
    22892692      CoinBigIndex start = columnStart[iColumn];
    22902693      CoinBigIndex end = start + columnLength[iColumn];
    22912694      if (end!=columnStart[iColumn+1]) {
    2292         hasGaps_=true;
     2695        flags_ |= 2;
    22932696        break;
    22942697      }
     
    22962699    return true;
    22972700  }
    2298   assert (check==15);
    2299   int * mark = new int [numberRows];
    2300   int i;
    2301   for (i=0;i<numberRows;i++)
    2302     mark[i]=-1;
    2303   for (iColumn=0;iColumn<numberColumns;iColumn++) {
    2304     CoinBigIndex j;
    2305     CoinBigIndex start = columnStart[iColumn];
    2306     CoinBigIndex end = start + columnLength[iColumn];
    2307     if (end!=columnStart[iColumn+1])
    2308       hasGaps_=true;
    2309     for (j=start;j<end;j++) {
    2310       double value = fabs(elementByColumn[j]);
    2311       int iRow = row[j];
    2312       if (iRow<0||iRow>=numberRows) {
     2701  assert (check==15||check==11);
     2702  if (check==15) {
     2703    int * mark = new int [numberRows];
     2704    int i;
     2705    for (i=0;i<numberRows;i++)
     2706      mark[i]=-1;
     2707    for (iColumn=0;iColumn<numberColumns;iColumn++) {
     2708      CoinBigIndex j;
     2709      CoinBigIndex start = columnStart[iColumn];
     2710      CoinBigIndex end = start + columnLength[iColumn];
     2711      if (end!=columnStart[iColumn+1])
     2712        flags_ |= 2;
     2713      for (j=start;j<end;j++) {
     2714        double value = fabs(elementByColumn[j]);
     2715        int iRow = row[j];
     2716        if (iRow<0||iRow>=numberRows) {
    23132717#ifndef COIN_BIG_INDEX
    2314         printf("Out of range %d %d %d %g\n",iColumn,j,row[j],elementByColumn[j]);
     2718          printf("Out of range %d %d %d %g\n",iColumn,j,row[j],elementByColumn[j]);
    23152719#elif COIN_BIG_INDEX==1
    2316         printf("Out of range %d %ld %d %g\n",iColumn,j,row[j],elementByColumn[j]);
     2720          printf("Out of range %d %ld %d %g\n",iColumn,j,row[j],elementByColumn[j]);
    23172721#else
    2318         printf("Out of range %d %lld %d %g\n",iColumn,j,row[j],elementByColumn[j]);
    2319 #endif
    2320         return false;
    2321       }
    2322       if (mark[iRow]==-1) {
    2323         mark[iRow]=j;
    2324       } else {
    2325         // duplicate
    2326         numberDuplicate++;
    2327       }
    2328       //printf("%d %d %d %g\n",iColumn,j,row[j],elementByColumn[j]);
    2329       if (!value)
    2330         zeroElements_ = true; // there are zero elements
    2331       if (value<smallest) {
    2332         numberSmall++;
    2333       } else if (!(value<=largest)) {
    2334         numberLarge++;
    2335         if (firstBadColumn<0) {
    2336           firstBadColumn=iColumn;
    2337           firstBadRow=row[j];
    2338           firstBadElement=elementByColumn[j];
    2339         }
    2340       }
    2341     }
    2342     //clear mark
    2343     for (j=columnStart[iColumn];
    2344          j<columnStart[iColumn]+columnLength[iColumn];j++) {
    2345       int iRow = row[j];
    2346       mark[iRow]=-1;
    2347     }
    2348   }
    2349   delete [] mark;
     2722          printf("Out of range %d %lld %d %g\n",iColumn,j,row[j],elementByColumn[j]);
     2723#endif
     2724          return false;
     2725        }
     2726        if (mark[iRow]==-1) {
     2727          mark[iRow]=j;
     2728        } else {
     2729          // duplicate
     2730          numberDuplicate++;
     2731        }
     2732        //printf("%d %d %d %g\n",iColumn,j,row[j],elementByColumn[j]);
     2733        if (!value)
     2734          flags_ |= 1; // there are zero elements
     2735        if (value<smallest) {
     2736          numberSmall++;
     2737        } else if (!(value<=largest)) {
     2738          numberLarge++;
     2739          if (firstBadColumn<0) {
     2740            firstBadColumn=iColumn;
     2741            firstBadRow=row[j];
     2742            firstBadElement=elementByColumn[j];
     2743          }
     2744        }
     2745      }
     2746      //clear mark
     2747      for (j=columnStart[iColumn];
     2748           j<columnStart[iColumn]+columnLength[iColumn];j++) {
     2749        int iRow = row[j];
     2750        mark[iRow]=-1;
     2751      }
     2752    }
     2753    delete [] mark;
     2754  } else {
     2755    // just check for out of range - not for duplicates
     2756    for (iColumn=0;iColumn<numberColumns;iColumn++) {
     2757      CoinBigIndex j;
     2758      CoinBigIndex start = columnStart[iColumn];
     2759      CoinBigIndex end = start + columnLength[iColumn];
     2760      if (end!=columnStart[iColumn+1])
     2761        flags_ |= 2;
     2762      for (j=start;j<end;j++) {
     2763        double value = fabs(elementByColumn[j]);
     2764        int iRow = row[j];
     2765        if (iRow<0||iRow>=numberRows) {
     2766#ifndef COIN_BIG_INDEX
     2767          printf("Out of range %d %d %d %g\n",iColumn,j,row[j],elementByColumn[j]);
     2768#elif COIN_BIG_INDEX==1
     2769          printf("Out of range %d %ld %d %g\n",iColumn,j,row[j],elementByColumn[j]);
     2770#else
     2771          printf("Out of range %d %lld %d %g\n",iColumn,j,row[j],elementByColumn[j]);
     2772#endif
     2773          return false;
     2774        }
     2775        if (!value)
     2776          flags_ |= 1; // there are zero elements
     2777        if (value<smallest) {
     2778          numberSmall++;
     2779        } else if (!(value<=largest)) {
     2780          numberLarge++;
     2781          if (firstBadColumn<0) {
     2782            firstBadColumn=iColumn;
     2783            firstBadRow=iRow;
     2784            firstBadElement=value;
     2785          }
     2786        }
     2787      }
     2788    }
     2789  }
    23502790  if (numberLarge) {
    23512791    model->messageHandler()->message(CLP_BAD_MATRIX,model->messages())
     
    23692809  // If smallest >0.0 then there can't be zero elements
    23702810  if (smallest>0.0)
    2371     zeroElements_=false;
     2811    flags_ &= ~1;;
    23722812  return true;
    23732813}
     
    26883128  rhsOffset(model,true);
    26893129}
     3130// Gets rid of special copies
     3131void
     3132ClpPackedMatrix::clearCopies()
     3133{
     3134  delete rowCopy_;
     3135  delete columnCopy_;
     3136  rowCopy_=NULL;
     3137  columnCopy_=NULL;
     3138  flags_ &= ~(4+8);
     3139 }
    26903140// makes sure active columns correct
    26913141int
     
    26933143{
    26943144  numberActiveColumns_ = matrix_->getNumCols();
     3145#if 0
     3146  ClpMatrixBase * rowCopyBase = reverseOrderedCopy();
     3147  ClpPackedMatrix* rowCopy =
     3148    dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
     3149  // Make sure it is really a ClpPackedMatrix
     3150  assert (rowCopy!=NULL);
     3151
     3152  const int * column = rowCopy->getIndices();
     3153  const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
     3154  const int * rowLength = rowCopy->getVectorLengths();
     3155  const double * element = rowCopy->getElements();
     3156  for (int i=0;i<rowCopy->getNumRows();i++) {
     3157    if (!rowLength[i])
     3158      printf("zero row %d\n",i);
     3159  }
     3160  delete rowCopy;
     3161#endif
    26953162  return 0;
    26963163}
     
    27773244ClpPackedMatrix::reallyScale(const double * rowScale, const double * columnScale)
    27783245{
     3246  clearCopies();
    27793247  int numberColumns = matrix_->getNumCols();
    27803248  const int * row = matrix_->getIndices();
     
    27983266  if (matrix_->getNumCols())
    27993267    matrix_->deleteCols(numDel,indDel);
     3268  clearCopies();
    28003269  numberActiveColumns_ = matrix_->getNumCols();
    28013270  // may now have gaps
    2802   hasGaps_=true;
     3271  flags_ |= 2;
    28033272  matrix_->setExtraGap(1.0e-50);
    28043273}
     
    28093278  if (matrix_->getNumRows())
    28103279    matrix_->deleteRows(numDel,indDel);
     3280  clearCopies();
    28113281  numberActiveColumns_ = matrix_->getNumCols();
    28123282  // may now have gaps
    2813   hasGaps_=true;
     3283  flags_ |= 2;
    28143284  matrix_->setExtraGap(1.0e-50);
    28153285}
     
    28213291  matrix_->appendCols(number,columns);
    28223292  numberActiveColumns_ = matrix_->getNumCols();
     3293  clearCopies();
    28233294}
    28243295// Append Rows
     
    28293300  numberActiveColumns_ = matrix_->getNumCols();
    28303301  // may now have gaps
    2831   hasGaps_=true;
     3302  flags_ |= 2;
     3303  clearCopies();
    28323304}
    28333305#endif
     
    28643336    numberErrors=matrix_->appendCols(number,starts,index,element,numberOther);
    28653337  }
     3338  clearCopies();
    28663339  numberActiveColumns_ = matrix_->getNumCols();
    28673340  return numberErrors;
     
    28763349    delete rowCopy_;
    28773350    rowCopy_=NULL;
     3351    flags_ &= ~4;
     3352  } else {
     3353    flags_ |= 4;
     3354  }
     3355}
     3356void
     3357ClpPackedMatrix::specialColumnCopy(ClpSimplex * model)
     3358{
     3359  delete columnCopy_;
     3360  if ((flags_&8)!=0)
     3361    columnCopy_ = new ClpPackedMatrix3(model,matrix_);
     3362  else
     3363    columnCopy_=NULL;
     3364}
     3365// Say we don't want special column copy
     3366void
     3367ClpPackedMatrix::releaseSpecialColumnCopy()
     3368{
     3369  flags_ &= ~8;
     3370  delete columnCopy_;
     3371  columnCopy_=NULL;
     3372}
     3373// Correct sequence in and out to give true value
     3374void
     3375ClpPackedMatrix::correctSequence(const ClpSimplex * model,int & sequenceIn, int & sequenceOut)
     3376{
     3377  if (columnCopy_) {
     3378    if (sequenceIn!=-999) {
     3379      if (sequenceIn!=sequenceOut) {
     3380        if (sequenceIn<numberActiveColumns_)
     3381          columnCopy_->swapOne(model,this,sequenceIn);
     3382        if (sequenceOut<numberActiveColumns_)
     3383          columnCopy_->swapOne(model,this,sequenceOut);
     3384      }
     3385    } else {
     3386      // do all
     3387      columnCopy_->sortBlocks(model);
     3388    }
     3389  }
     3390}
     3391// Check validity
     3392void
     3393ClpPackedMatrix::checkFlags() const
     3394{
     3395  int iColumn;
     3396  // get matrix data pointers
     3397  //const int * row = matrix_->getIndices();
     3398  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
     3399  const int * columnLength = matrix_->getVectorLengths();
     3400  const double * elementByColumn = matrix_->getElements();
     3401  if (!zeros()) {
     3402    for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
     3403      CoinBigIndex j;
     3404      for (j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
     3405        if (!elementByColumn[j])
     3406          abort();
     3407      }
     3408    }
     3409  }
     3410  if ((flags_&2)==0) {
     3411    for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
     3412      if (columnStart[iColumn+1]!=columnStart[iColumn]+columnLength[iColumn]) {
     3413        abort();
     3414      }
     3415    }
    28783416  }
    28793417}
     
    29273465  // tune
    29283466#if 0
    2929   int chunkY[5]={4096,8192,16384,32768,65535};
     3467  int chunkY[7]={1024,2048,4096,8192,16384,32768,65535};
    29303468  int its = model->maximumIterations();
    29313469  if (its>=1000000&&its<1000999) {
    29323470    its -= 1000000;
    29333471    its = its/10;
    2934     if (its>=5) abort();
     3472    if (its>=7) abort();
    29353473    chunk=chunkY[its];
    29363474    printf("chunk size %d\n",chunk);
     
    37174255  }
    37184256}
     4257/* Default constructor. */
     4258ClpPackedMatrix3::ClpPackedMatrix3()
     4259  : numberBlocks_(0),
     4260    numberColumns_(0),
     4261    column_(NULL),
     4262    start_(NULL),
     4263    row_(NULL),
     4264    element_(NULL),
     4265    block_(NULL)
     4266{
     4267}
     4268/* Constructor from copy. */
     4269ClpPackedMatrix3::ClpPackedMatrix3(ClpSimplex * model,const CoinPackedMatrix * columnCopy)
     4270  : numberBlocks_(0),
     4271    numberColumns_(0),
     4272    column_(NULL),
     4273    start_(NULL),
     4274    row_(NULL),
     4275    element_(NULL),
     4276    block_(NULL)
     4277{
     4278#define MINBLOCK 6
     4279#define MAXBLOCK 100
     4280#define MAXUNROLL 10
     4281  numberColumns_ = columnCopy->getNumCols();
     4282  int numberRows = columnCopy->getNumRows();
     4283  int * counts = new int[numberRows+1];
     4284  CoinZeroN(counts,numberRows+1);
     4285  CoinBigIndex nels=0;
     4286  CoinBigIndex nZeroEl=0;
     4287  int iColumn;
     4288  // get matrix data pointers
     4289  const int * row = columnCopy->getIndices();
     4290  const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
     4291  const int * columnLength = columnCopy->getVectorLengths();
     4292  const double * elementByColumn = columnCopy->getElements();
     4293  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     4294    CoinBigIndex start = columnStart[iColumn];
     4295    int n = columnLength[iColumn];
     4296    CoinBigIndex end = start + n;
     4297    int kZero=0;
     4298    for (CoinBigIndex j=start;j<end;j++) {
     4299      if(!elementByColumn[j])
     4300        kZero++;
     4301    }
     4302    n -= kZero;
     4303    nZeroEl += kZero;
     4304    nels += n;
     4305    counts[n]++;
     4306  }
     4307  int nZeroColumns = counts[0];
     4308  counts[0]=-1;
     4309  numberColumns_ -= nZeroColumns;
     4310  column_ = new int [2*numberColumns_];
     4311  int * lookup = column_ + numberColumns_;
     4312  row_ = new int[nels];
     4313  element_ = new double[nels];
     4314  int nOdd=0;
     4315  CoinBigIndex nInOdd=0;
     4316  int i;
     4317  for (i=1;i<=numberRows;i++) {
     4318    int n = counts[i];
     4319    if (n) {
     4320      if (n<MINBLOCK||i>MAXBLOCK) {
     4321        nOdd += n;
     4322        counts[i]=-1;
     4323        nInOdd += n*i;
     4324      } else {
     4325        numberBlocks_++;
     4326      }
     4327    } else {
     4328      counts[i]=-1;
     4329    }
     4330  }
     4331  start_ = new CoinBigIndex [nOdd+1];
     4332  // even if no blocks do a dummy one
     4333  numberBlocks_ = CoinMax(numberBlocks_,1);
     4334  block_ = (blockStruct *) new char [numberBlocks_*sizeof(blockStruct)];
     4335  memset(block_,0,numberBlocks_*sizeof(blockStruct));
     4336  // Fill in what we can
     4337  int nTotal=nOdd;
     4338  // in case no blocks
     4339  block_->startIndices_=nTotal;
     4340  nels=nInOdd;
     4341  int nBlock=0;
     4342  for (i=0;i<=CoinMin(MAXBLOCK,numberRows);i++) {
     4343    if (counts[i]>0) {
     4344      blockStruct * block = block_ + nBlock;
     4345      int n=counts[i];
     4346      counts[i]= nBlock; // backward pointer
     4347      nBlock++;
     4348      block->startIndices_ = nTotal;
     4349      block->startElements_ = nels;
     4350      block->numberElements_ = i;
     4351      // up counts
     4352      nTotal += n;
     4353      nels += n*i;
     4354    }
     4355  }
     4356  // fill
     4357  start_[0]=0;
     4358  nOdd=0;
     4359  nInOdd=0;
     4360  const double * columnScale = model->columnScale();
     4361  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     4362    CoinBigIndex start = columnStart[iColumn];
     4363    int n = columnLength[iColumn];
     4364    CoinBigIndex end = start + n;
     4365    int kZero=0;
     4366    for (CoinBigIndex j=start;j<end;j++) {
     4367      if(!elementByColumn[j])
     4368        kZero++;
     4369    }
     4370    n -= kZero;
     4371    if (n) {
     4372      int iBlock = counts[n];
     4373      if (iBlock>=0) {
     4374        blockStruct * block = block_ + iBlock;
     4375        int k = block->numberInBlock_;
     4376        block->numberInBlock_ ++;
     4377        column_[block->startIndices_+k]=iColumn;
     4378        lookup[iColumn]=k;
     4379        CoinBigIndex put = block->startElements_+k*n;
     4380        for (CoinBigIndex j=start;j<end;j++) {
     4381          double value = elementByColumn[j];
     4382          if(value) {
     4383            if (columnScale)
     4384              value *= columnScale[iColumn];
     4385            element_[put] = value;
     4386            row_[put++]=row[j];
     4387          }
     4388        }
     4389      } else {
     4390        // odd ones
     4391        for (CoinBigIndex j=start;j<end;j++) {
     4392          double value = elementByColumn[j];
     4393          if(value) {
     4394            if (columnScale)
     4395              value *= columnScale[iColumn];
     4396            element_[nInOdd] = value;
     4397            row_[nInOdd++]=row[j];
     4398          }
     4399        }
     4400        column_[nOdd]=iColumn;
     4401        lookup[iColumn]=-1;
     4402        nOdd++;
     4403        start_[nOdd]=nInOdd;
     4404      }
     4405    }
     4406  }
     4407  delete [] counts;
     4408}
     4409/* Destructor */
     4410ClpPackedMatrix3::~ClpPackedMatrix3()
     4411{
     4412  delete [] column_;
     4413  delete [] start_;
     4414  delete [] row_;
     4415  delete [] element_;
     4416  delete [] block_;
     4417}
     4418/* The copy constructor. */
     4419ClpPackedMatrix3::ClpPackedMatrix3(const ClpPackedMatrix3 & rhs)
     4420  : numberBlocks_(rhs.numberBlocks_),
     4421    numberColumns_(rhs.numberColumns_),
     4422    column_(NULL),
     4423    start_(NULL),
     4424    row_(NULL),
     4425    element_(NULL),
     4426    block_(NULL)
     4427{
     4428  if (rhs.numberBlocks_) {
     4429    block_ = CoinCopyOfArray(rhs.block_,numberBlocks_);
     4430    column_ = CoinCopyOfArray(rhs.column_,2*numberColumns_);
     4431    int numberOdd = block_->startIndices_;
     4432    start_ = CoinCopyOfArray(rhs.start_,numberOdd+1);
     4433    blockStruct * lastBlock = block_ + (numberBlocks_-1);
     4434    CoinBigIndex numberElements = lastBlock->startElements_+
     4435      lastBlock->numberInBlock_*lastBlock->numberElements_;
     4436    row_ = CoinCopyOfArray(rhs.row_,numberElements);
     4437    element_ = CoinCopyOfArray(rhs.element_,numberElements);
     4438  }
     4439}
     4440ClpPackedMatrix3&
     4441ClpPackedMatrix3::operator=(const ClpPackedMatrix3 & rhs)
     4442{
     4443  if (this != &rhs) {
     4444    delete [] column_;
     4445    delete [] start_;
     4446    delete [] row_;
     4447    delete [] element_;
     4448    delete [] block_;
     4449    numberBlocks_ = rhs.numberBlocks_;
     4450    numberColumns_ = rhs.numberColumns_;
     4451    if (rhs.numberBlocks_) {
     4452      block_ = CoinCopyOfArray(rhs.block_,numberBlocks_);
     4453      column_ = CoinCopyOfArray(rhs.column_,2*numberColumns_);
     4454      int numberOdd = block_->startIndices_;
     4455      start_ = CoinCopyOfArray(rhs.start_,numberOdd+1);
     4456      blockStruct * lastBlock = block_ + (numberBlocks_-1);
     4457      CoinBigIndex numberElements = lastBlock->startElements_+
     4458        lastBlock->numberInBlock_*lastBlock->numberElements_;
     4459      row_ = CoinCopyOfArray(rhs.row_,numberElements);
     4460      element_ = CoinCopyOfArray(rhs.element_,numberElements);
     4461    } else {
     4462      column_ = NULL;
     4463      start_ = NULL;
     4464      row_ = NULL;
     4465      element_ = NULL;
     4466      block_ = NULL;
     4467    }
     4468  }
     4469  return *this;
     4470}
     4471/* Sort blocks */
     4472void
     4473ClpPackedMatrix3::sortBlocks(const ClpSimplex * model)
     4474{
     4475  int * lookup = column_ + numberColumns_;
     4476  for (int iBlock=0;iBlock<numberBlocks_;iBlock++) {
     4477    blockStruct * block = block_ + iBlock;
     4478    int numberInBlock = block->numberInBlock_;
     4479    int nel = block->numberElements_;
     4480    int * row = row_ + block->startElements_;
     4481    double * element = element_ + block->startElements_;
     4482    int * column = column_ + block->startIndices_;
     4483    int lastPrice=0;
     4484    int firstNotPrice=numberInBlock-1;
     4485    while (lastPrice<=firstNotPrice) {
     4486      // find first basic or fixed
     4487      int iColumn=numberInBlock;
     4488      for (;lastPrice<=firstNotPrice;lastPrice++) {
     4489        iColumn = column[lastPrice];
     4490        if (model->getColumnStatus(iColumn)==ClpSimplex::basic||
     4491            model->getColumnStatus(iColumn)==ClpSimplex::isFixed)
     4492          break;
     4493      }
     4494      // find last non basic or fixed
     4495      int jColumn=-1;
     4496      for (;firstNotPrice>lastPrice;firstNotPrice--) {
     4497        jColumn = column[firstNotPrice];
     4498        if (model->getColumnStatus(jColumn)!=ClpSimplex::basic&&
     4499            model->getColumnStatus(jColumn)!=ClpSimplex::isFixed)
     4500          break;
     4501      }
     4502      if (firstNotPrice>lastPrice) {
     4503        assert (column[lastPrice]==iColumn);
     4504        assert (column[firstNotPrice]==jColumn);
     4505        // need to swap
     4506        column[firstNotPrice]=iColumn;
     4507        lookup[iColumn]=firstNotPrice;
     4508        column[lastPrice]=jColumn;
     4509        lookup[jColumn]=lastPrice;
     4510        double * elementA = element + lastPrice*nel;
     4511        int * rowA = row + lastPrice*nel;
     4512        double * elementB = element + firstNotPrice*nel;
     4513        int * rowB = row + firstNotPrice*nel;
     4514        for (int i=0;i<nel;i++) {
     4515          int temp = rowA[i];
     4516          double tempE = elementA[i];
     4517          rowA[i]=rowB[i];
     4518          elementA[i]=elementB[i];
     4519          rowB[i]=temp;
     4520          elementB[i]=tempE;
     4521        }
     4522        firstNotPrice--;
     4523        lastPrice++;
     4524      } else if (lastPrice==firstNotPrice) {
     4525        // make sure correct side
     4526        iColumn = column[lastPrice];
     4527        if (model->getColumnStatus(iColumn)!=ClpSimplex::basic&&
     4528            model->getColumnStatus(iColumn)!=ClpSimplex::isFixed)
     4529          lastPrice++;
     4530        break;
     4531      }
     4532    }
     4533    block->numberPrice_=lastPrice;
     4534#ifndef NDEBUG
     4535    // check
     4536    int i;
     4537    for (i=0;i<lastPrice;i++) {
     4538      int iColumn = column[i];
     4539      assert (model->getColumnStatus(iColumn)!=ClpSimplex::basic&&
     4540              model->getColumnStatus(iColumn)!=ClpSimplex::isFixed);
     4541      assert (lookup[iColumn]==i);
     4542    }
     4543    for (;i<numberInBlock;i++) {
     4544      int iColumn = column[i];
     4545      assert (model->getColumnStatus(iColumn)==ClpSimplex::basic||
     4546              model->getColumnStatus(iColumn)==ClpSimplex::isFixed);
     4547      assert (lookup[iColumn]==i);
     4548    }
     4549#endif
     4550  }
     4551}
     4552// Swap one variable
     4553void
     4554ClpPackedMatrix3::swapOne(const ClpSimplex * model,const ClpPackedMatrix * matrix,
     4555                          int iColumn)
     4556{
     4557  int * lookup = column_ + numberColumns_;
     4558  // position in block
     4559  int kA = lookup[iColumn];
     4560  if (kA<0)
     4561    return; // odd one
     4562  // get matrix data pointers
     4563  const CoinPackedMatrix * columnCopy = matrix->getPackedMatrix();
     4564  //const int * row = columnCopy->getIndices();
     4565  const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
     4566  const int * columnLength = columnCopy->getVectorLengths();
     4567  const double * elementByColumn = columnCopy->getElements();
     4568  CoinBigIndex start = columnStart[iColumn];
     4569  int n = columnLength[iColumn];
     4570  if (matrix->zeros()) {
     4571    CoinBigIndex end = start + n;
     4572    for (CoinBigIndex j=start;j<end;j++) {
     4573      if(!elementByColumn[j])
     4574        n--;
     4575    }
     4576  }
     4577  // find block - could do binary search
     4578  int iBlock=CoinMin(n,numberBlocks_)-1;
     4579  while (block_[iBlock].numberElements_!=n)
     4580    iBlock--;
     4581  blockStruct * block = block_ + iBlock;
     4582  int nel = block->numberElements_;
     4583  int * row = row_ + block->startElements_;
     4584  double * element = element_ + block->startElements_;
     4585  int * column = column_ + block->startIndices_;
     4586  assert (column[kA]==iColumn);
     4587  bool moveUp = (model->getColumnStatus(iColumn)==ClpSimplex::basic||
     4588                 model->getColumnStatus(iColumn)==ClpSimplex::isFixed);
     4589  int lastPrice=block->numberPrice_;
     4590  int kB;
     4591  if (moveUp) {
     4592    kB=lastPrice-1;
     4593    block->numberPrice_--;
     4594  } else {
     4595    kB=lastPrice;
     4596    block->numberPrice_++;
     4597  }
     4598  int jColumn = column[kB];
     4599  column[kA]=jColumn;
     4600  lookup[jColumn]=kA;
     4601  column[kB]=iColumn;
     4602  lookup[iColumn]=kB;
     4603  double * elementA = element + kB*nel;
     4604  int * rowA = row + kB*nel;
     4605  double * elementB = element + kA*nel;
     4606  int * rowB = row + kA*nel;
     4607  for (int i=0;i<nel;i++) {
     4608    int temp = rowA[i];
     4609    double tempE = elementA[i];
     4610    rowA[i]=rowB[i];
     4611    elementA[i]=elementB[i];
     4612    rowB[i]=temp;
     4613    elementB[i]=tempE;
     4614  }
     4615#if 1
     4616#ifndef NDEBUG
     4617  // check
     4618  int i;
     4619  for (i=0;i<block->numberPrice_;i++) {
     4620    int iColumn = column[i];
     4621    if (iColumn!=model->sequenceIn()&&iColumn!=model->sequenceOut())
     4622      assert (model->getColumnStatus(iColumn)!=ClpSimplex::basic&&
     4623              model->getColumnStatus(iColumn)!=ClpSimplex::isFixed);
     4624    assert (lookup[iColumn]==i);
     4625  }
     4626  int numberInBlock = block->numberInBlock_;
     4627  for (;i<numberInBlock;i++) {
     4628    int iColumn = column[i];
     4629    if (iColumn!=model->sequenceIn()&&iColumn!=model->sequenceOut())
     4630      assert (model->getColumnStatus(iColumn)==ClpSimplex::basic||
     4631              model->getColumnStatus(iColumn)==ClpSimplex::isFixed);
     4632    assert (lookup[iColumn]==i);
     4633  }
     4634#endif
     4635#endif
     4636}
     4637/* Return <code>x * -1 * A in <code>z</code>.
     4638   Note - x packed and z will be packed mode
     4639   Squashes small elements and knows about ClpSimplex */
     4640void
     4641ClpPackedMatrix3::transposeTimes(const ClpSimplex * model,
     4642                                 const double * pi,
     4643                                 CoinIndexedVector * output) const
     4644{
     4645  int numberNonZero=0;
     4646  int * index = output->getIndices();
     4647  double * array = output->denseVector();
     4648  double zeroTolerance = model->factorization()->zeroTolerance();
     4649  double value = 0.0;
     4650  CoinBigIndex j;
     4651  int numberOdd = block_->startIndices_;
     4652  if (numberOdd) {
     4653    // A) as probably long may be worth unrolling
     4654    CoinBigIndex end = start_[1];
     4655    for (j=start_[0]; j<end;j++) {
     4656      int iRow = row_[j];
     4657      value += pi[iRow]*element_[j];
     4658    }
     4659    int iColumn;
     4660    // int jColumn=column_[0];
     4661   
     4662    for (iColumn=0;iColumn<numberOdd-1;iColumn++) {
     4663      CoinBigIndex start = end;
     4664      end = start_[iColumn+2];
     4665      if (fabs(value)>zeroTolerance) {
     4666        array[numberNonZero]=value;
     4667        index[numberNonZero++]=column_[iColumn];
     4668        //index[numberNonZero++]=jColumn;
     4669      }
     4670      // jColumn = column_[iColumn+1];
     4671      value = 0.0;
     4672      //if (model->getColumnStatus(jColumn)!=ClpSimplex::basic) {
     4673        for (j=start; j<end;j++) {
     4674          int iRow = row_[j];
     4675          value += pi[iRow]*element_[j];
     4676        }
     4677        //}
     4678    }
     4679    if (fabs(value)>zeroTolerance) {
     4680      array[numberNonZero]=value;
     4681      index[numberNonZero++]=column_[iColumn];
     4682      //index[numberNonZero++]=jColumn;
     4683    }
     4684  }
     4685  for (int iBlock=0;iBlock<numberBlocks_;iBlock++) {
     4686    // B) Can sort so just do nonbasic (and nonfixed)
     4687    // C) Can do two at a time (if so put odd one into start_)
     4688    // D) can use switch
     4689    blockStruct * block = block_ + iBlock;
     4690    //int numberPrice = block->numberInBlock_;
     4691    int numberPrice = block->numberPrice_;
     4692    int nel = block->numberElements_;
     4693    int * row = row_ + block->startElements_;
     4694    double * element = element_ + block->startElements_;
     4695    int * column = column_ + block->startIndices_;
     4696#if 0
     4697    // two at a time
     4698    if ((numberPrice&1)!=0) {
     4699      double value = 0.0;
     4700      int nel2=nel;
     4701      for (; nel2;nel2--) {
     4702        int iRow = *row++;
     4703        value += pi[iRow]*(*element++);
     4704      }
     4705      if (fabs(value)>zeroTolerance) {
     4706        array[numberNonZero]=value;
     4707        index[numberNonZero++]=*column;
     4708      }
     4709      column++;
     4710    }
     4711    numberPrice = numberPrice>>1;
     4712    switch ((nel%2)) {
     4713      // 2 k +0
     4714    case 0:
     4715      for (;numberPrice;numberPrice--) {
     4716        double value0 = 0.0;
     4717        double value1 = 0.0;
     4718        int nel2=nel;
     4719        for (; nel2;nel2--) {
     4720          int iRow0 = *row;
     4721          int iRow1 = *(row+nel);
     4722          row++;
     4723          double element0 = *element;
     4724          double element1 = *(element+nel);
     4725          element++;
     4726          value0 += pi[iRow0]*element0;
     4727          value1 += pi[iRow1]*element1;
     4728        }
     4729        row += nel;
     4730        element += nel;
     4731        if (fabs(value0)>zeroTolerance) {
     4732          array[numberNonZero]=value0;
     4733          index[numberNonZero++]=*column;
     4734        }
     4735        column++;
     4736        if (fabs(value1)>zeroTolerance) {
     4737          array[numberNonZero]=value1;
     4738          index[numberNonZero++]=*column;
     4739        }
     4740        column++;
     4741      } 
     4742      break;
     4743      // 2 k +1
     4744    case 1:
     4745      for (;numberPrice;numberPrice--) {
     4746        double value0;
     4747        double value1;
     4748        int nel2=nel-1;
     4749        {
     4750          int iRow0 = row[0];
     4751          int iRow1 = row[nel];
     4752          double pi0 = pi[iRow0];
     4753          double pi1 = pi[iRow1];
     4754          value0 = pi0*element[0];
     4755          value1 = pi1*element[nel];
     4756          row++;
     4757          element++;
     4758        }
     4759        for (; nel2;nel2--) {
     4760          int iRow0 = *row;
     4761          int iRow1 = *(row+nel);
     4762          row++;
     4763          double element0 = *element;
     4764          double element1 = *(element+nel);
     4765          element++;
     4766          value0 += pi[iRow0]*element0;
     4767          value1 += pi[iRow1]*element1;
     4768        }
     4769        row += nel;
     4770        element += nel;
     4771        if (fabs(value0)>zeroTolerance) {
     4772          array[numberNonZero]=value0;
     4773          index[numberNonZero++]=*column;
     4774        }
     4775        column++;
     4776        if (fabs(value1)>zeroTolerance) {
     4777          array[numberNonZero]=value1;
     4778          index[numberNonZero++]=*column;
     4779        }
     4780        column++;
     4781      }
     4782      break;
     4783    }
     4784#else
     4785    for (;numberPrice;numberPrice--) {
     4786      double value = 0.0;
     4787      int nel2=nel;
     4788      for (; nel2;nel2--) {
     4789        int iRow = *row++;
     4790        value += pi[iRow]*(*element++);
     4791      }
     4792      if (fabs(value)>zeroTolerance) {
     4793        array[numberNonZero]=value;
     4794        index[numberNonZero++]=*column;
     4795      }
     4796      column++;
     4797    }
     4798#endif
     4799  }
     4800  output->setNumElements(numberNonZero);
     4801}
     4802// Updates two arrays for steepest
     4803void
     4804ClpPackedMatrix3::transposeTimes2(const ClpSimplex * model,
     4805                                  const double * pi, CoinIndexedVector * output,
     4806                                  const double * piWeight,
     4807                                  double referenceIn, double devex,
     4808                                  // Array for exact devex to say what is in reference framework
     4809                                  unsigned int * reference,
     4810                                  double * weights, double scaleFactor)
     4811{
     4812  int numberNonZero=0;
     4813  int * index = output->getIndices();
     4814  double * array = output->denseVector();
     4815  double zeroTolerance = model->factorization()->zeroTolerance();
     4816  double value = 0.0;
     4817  bool killDjs = (scaleFactor==0.0);
     4818  if (!scaleFactor)
     4819    scaleFactor=1.0;
     4820  int numberOdd = block_->startIndices_;
     4821  int iColumn;
     4822  CoinBigIndex end = start_[0];
     4823  for (iColumn=0;iColumn<numberOdd;iColumn++) {
     4824    CoinBigIndex start = end;
     4825    CoinBigIndex j;
     4826    int jColumn = column_[iColumn];
     4827    end = start_[iColumn+1];
     4828    value = 0.0;
     4829    if (model->getColumnStatus(jColumn)!=ClpSimplex::basic) {
     4830      for (j=start; j<end;j++) {
     4831        int iRow = row_[j];
     4832        value -= pi[iRow]*element_[j];
     4833      }
     4834      if (fabs(value)>zeroTolerance) {
     4835        // and do other array
     4836        double modification = 0.0;
     4837        for (j=start; j<end;j++) {
     4838          int iRow = row_[j];
     4839          modification += piWeight[iRow]*element_[j];
     4840        }
     4841        double thisWeight = weights[jColumn];
     4842        double pivot = value*scaleFactor;
     4843        double pivotSquared = pivot * pivot;
     4844        thisWeight += pivotSquared * devex + pivot * modification;
     4845        if (thisWeight<DEVEX_TRY_NORM) {
     4846          if (referenceIn<0.0) {
     4847            // steepest
     4848            thisWeight = CoinMax(DEVEX_TRY_NORM,DEVEX_ADD_ONE+pivotSquared);
     4849          } else {
     4850            // exact
     4851            thisWeight = referenceIn*pivotSquared;
     4852            if (reference(jColumn))
     4853              thisWeight += 1.0;
     4854            thisWeight = CoinMax(thisWeight,DEVEX_TRY_NORM);
     4855          }
     4856        }
     4857        weights[jColumn] = thisWeight;
     4858        if (!killDjs) {
     4859          array[numberNonZero]=value;
     4860          index[numberNonZero++]=jColumn;
     4861        }
     4862      }
     4863    }
     4864  }
     4865  for (int iBlock=0;iBlock<numberBlocks_;iBlock++) {
     4866    // B) Can sort so just do nonbasic (and nonfixed)
     4867    // C) Can do two at a time (if so put odd one into start_)
     4868    // D) can use switch
     4869    blockStruct * block = block_ + iBlock;
     4870    //int numberPrice = block->numberInBlock_;
     4871    int numberPrice = block->numberPrice_;
     4872    int nel = block->numberElements_;
     4873    int * row = row_ + block->startElements_;
     4874    double * element = element_ + block->startElements_;
     4875    int * column = column_ + block->startIndices_;
     4876    for (;numberPrice;numberPrice--) {
     4877      double value = 0.0;
     4878      int nel2=nel;
     4879      for (; nel2;nel2--) {
     4880        int iRow = *row++;
     4881        value -= pi[iRow]*(*element++);
     4882      }
     4883      if (fabs(value)>zeroTolerance) {
     4884        int jColumn = *column;
     4885        // back to beginning
     4886        row -= nel;
     4887        element -= nel;
     4888        // and do other array
     4889        double modification = 0.0;
     4890        nel2=nel;
     4891        for (; nel2;nel2--) {
     4892          int iRow = *row++;
     4893          modification += piWeight[iRow]*(*element++);
     4894        }
     4895        double thisWeight = weights[jColumn];
     4896        double pivot = value*scaleFactor;
     4897        double pivotSquared = pivot * pivot;
     4898        thisWeight += pivotSquared * devex + pivot * modification;
     4899        if (thisWeight<DEVEX_TRY_NORM) {
     4900          if (referenceIn<0.0) {
     4901            // steepest
     4902            thisWeight = CoinMax(DEVEX_TRY_NORM,DEVEX_ADD_ONE+pivotSquared);
     4903          } else {
     4904            // exact
     4905            thisWeight = referenceIn*pivotSquared;
     4906            if (reference(jColumn))
     4907              thisWeight += 1.0;
     4908            thisWeight = CoinMax(thisWeight,DEVEX_TRY_NORM);
     4909          }
     4910        }
     4911        weights[jColumn] = thisWeight;
     4912        if (!killDjs) {
     4913          array[numberNonZero]=value;
     4914          index[numberNonZero++]=jColumn;
     4915        }
     4916      }
     4917      column++;
     4918    }
     4919  }
     4920  output->setNumElements(numberNonZero);
     4921  output->setPackedMode(true);
     4922}
    37194923#ifdef CLP_ALL_ONE_FILE
    37204924#undef reference
  • trunk/Clp/src/ClpPackedMatrix.hpp

    r800 r1034  
    44#define ClpPackedMatrix_H
    55
    6 
    76#include "CoinPragma.hpp"
    87
     
    1615
    1716class ClpPackedMatrix2;
     17class ClpPackedMatrix3;
    1818class ClpPackedMatrix : public ClpMatrixBase {
    1919 
     
    260260  inline void setMatrixNull()
    261261  { matrix_=NULL;};
     262  /// Say we want special column copy
     263  inline void makeSpecialColumnCopy()
     264  { flags_ |= 8;};
     265  /// Say we don't want special column copy
     266  void releaseSpecialColumnCopy();
     267  /// Are there zeros?
     268  inline bool zeros() const
     269  { return ((flags_&1)!=0);};
     270  /// Do we want special column copy
     271  inline bool wantsSpecialColumnCopy() const
     272  { return ((flags_&8)!=0);};
    262273   //@}
    263274
     
    299310  /// make special row copy
    300311  void specialRowCopy(ClpSimplex * model,const ClpMatrixBase * rowCopy);
    301    //@}
     312  /// make special column copy
     313  void specialColumnCopy(ClpSimplex * model);
     314  /// Correct sequence in and out to give true value
     315  virtual void correctSequence(const ClpSimplex * model,int & sequenceIn, int & sequenceOut) ;
     316   //@}
     317private:
     318  /// Meat of transposeTimes by column when not scaled
     319  void gutsOfTransposeTimesUnscaled(const double * pi,CoinIndexedVector * output, const double tolerance) const;
     320  /// Meat of transposeTimes by column when scaled
     321  void gutsOfTransposeTimesScaled(const double * pi,const double * columnScale, CoinIndexedVector * output, const double tolerance) const;
     322  /// Meat of transposeTimes by row n > 2 if packed
     323  void gutsOfTransposeTimesByRowGE3(const CoinIndexedVector * piVector, CoinIndexedVector * output,
     324                                   CoinIndexedVector * spareVector, const double tolerance, const double scalar) const;
     325  /// Meat of transposeTimes by row n == 2 if packed
     326  void gutsOfTransposeTimesByRowEQ2(const CoinIndexedVector * piVector, CoinIndexedVector * output,
     327                                   CoinIndexedVector * spareVector, const double tolerance, const double scalar) const;
     328  /// Meat of transposeTimes by row n == 1 if packed
     329  void gutsOfTransposeTimesByRowEQ1(const CoinIndexedVector * piVector, CoinIndexedVector * output,
     330                                   const double tolerance, const double scalar) const;
     331  /// Gets rid of special copies
     332  void clearCopies();
    302333   
    303334   
    304335protected:
     336  /// Check validity
     337  void checkFlags() const;
    305338   /**@name Data members
    306339      The data members are protected to allow access for derived classes. */
     
    310343  /// number of active columns (normally same as number of columns)
    311344  int numberActiveColumns_;
    312   /// Zero element flag - set true if any zero elements
    313   bool zeroElements_;
    314   /// Gaps flag - set true if column start and length don't say contiguous
    315   bool hasGaps_;
     345  /** Flags -
     346      1 - has zero elements
     347      2 - has gaps
     348      4 - has special row copy
     349      8 - has special column copy
     350  */
     351  int flags_;
    316352  /// Special row copy
    317353  ClpPackedMatrix2 * rowCopy_;
     354  /// Special column copy
     355  ClpPackedMatrix3 * columnCopy_;
    318356   //@}
    319357};
     
    404442   //@}
    405443};
     444typedef struct {
     445  CoinBigIndex startElements_; // point to data
     446  int startIndices_; // point to column_
     447  int numberInBlock_;
     448  int numberPrice_; // at beginning
     449  int numberElements_; // number elements per column
     450} blockStruct;
     451class ClpPackedMatrix3 {
     452 
     453public:
     454  /**@name Useful methods */
     455  //@{
     456    /** Return <code>x * -1 * A in <code>z</code>.
     457        Note - x packed and z will be packed mode
     458        Squashes small elements and knows about ClpSimplex */
     459  void transposeTimes(const ClpSimplex * model,
     460                      const double * pi,
     461                      CoinIndexedVector * output) const;
     462  /// Updates two arrays for steepest
     463  void transposeTimes2(const ClpSimplex * model,
     464                       const double * pi, CoinIndexedVector * dj1,
     465                       const double * piWeight,
     466                       double referenceIn, double devex,
     467                       // Array for exact devex to say what is in reference framework
     468                       unsigned int * reference,
     469                       double * weights, double scaleFactor);
     470  //@}
     471
     472
     473  /**@name Constructors, destructor */
     474   //@{
     475   /** Default constructor. */
     476   ClpPackedMatrix3();
     477   /** Constructor from copy. */
     478   ClpPackedMatrix3(ClpSimplex * model,const CoinPackedMatrix * columnCopy);
     479   /** Destructor */
     480   virtual ~ClpPackedMatrix3();
     481   //@}
     482
     483   /**@name Copy method */
     484   //@{
     485   /** The copy constructor. */
     486   ClpPackedMatrix3(const ClpPackedMatrix3&);
     487   ClpPackedMatrix3& operator=(const ClpPackedMatrix3&);
     488   //@}
     489   /**@name Sort methods */
     490   //@{
     491   /** Sort blocks */
     492  void sortBlocks(const ClpSimplex * model);
     493  /// Swap one variable
     494  void swapOne(const ClpSimplex * model,const ClpPackedMatrix * matrix,
     495                int iColumn);
     496   //@}
     497   
     498   
     499protected:
     500   /**@name Data members
     501      The data members are protected to allow access for derived classes. */
     502   //@{
     503  /// Number of blocks
     504  int numberBlocks_;
     505  /// Number of columns
     506  int numberColumns_;
     507  /// Column indices and reverse lookup (within block)
     508  int * column_;
     509  /// Starts for odd/long vectors
     510  CoinBigIndex * start_;
     511  /// Rows
     512  int * row_;
     513  /// Elements
     514  double * element_;
     515  /// Blocks (ordinary start at 0 and go to first block)
     516  blockStruct * block_;
     517   //@}
     518};
    406519
    407520#endif
  • trunk/Clp/src/ClpPresolve.cpp

    r850 r1034  
    1313
    1414#include "CoinPackedMatrix.hpp"
     15#include "ClpPackedMatrix.hpp"
    1516#include "ClpSimplex.hpp"
    1617#ifndef SLIM_CLP
     
    481482    }
    482483
     484    if ((presolveActions_&16384)!=0)
     485      prob->setPresolveOptions(prob->presolveOptions()|16384);
    483486    // Check number rows dropped
    484487    int lastDropped=0;
     
    14321435  // Messages
    14331436  CoinMessages messages = originalModel->coinMessages();
     1437  // Only go round 100 times even if integer preprocessing
     1438  int totalPasses=100;
    14341439  while (result==-1) {
    14351440
     
    14561461    if (!keepIntegers)
    14571462      presolvedModel_->deleteIntegerInformation();
     1463    totalPasses--;
    14581464
    14591465    double ratio=2.0;
     
    16331639        }
    16341640        presolvedModel_->copyNames(rowNames,columnNames);
     1641      } else {
     1642        presolvedModel_->setLengthNames(0);
    16351643      }
    16361644#endif
     
    16961704                                                       <<numberChanges
    16971705                                                       <<CoinMessageEol;
    1698           if (!result) {
     1706          if (!result&&totalPasses>0) {
    16991707            result = -1; // round again
    17001708          }
  • trunk/Clp/src/ClpPresolve.hpp

    r754 r1034  
    124124  { if (doSingleton) presolveActions_  &= ~512; else presolveActions_ |= 512;};
    125125  /// Set whole group
    126   inline bool presolveActions() const
     126  inline int presolveActions() const
    127127  { return presolveActions_&0xffff;};
    128128  inline void setPresolveActions(int action)
  • trunk/Clp/src/ClpPrimalColumnSteepest.cpp

    r886 r1034  
    983983  if (sequenceOut>=0)
    984984    outgoingWeight=weights_[sequenceOut];
    985    
    986   // put row of tableau in rowArray and columnArray (packed)
    987   // get subset which have nonzero tableau elements
    988   transposeTimes2(updates,spareColumn1,alternateWeights_,spareColumn2,spareRow2,
    989                   -scaleFactor);
     985  // update row weights here so we can scale alternateWeights_
    990986  // update weights
    991987  double * weight;
     
    10081004    int iSequence = index[j];
    10091005    double value2 = updateBy[j];
    1010     updateBy[j]=0.0;
    10111006    ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
    10121007    double value;
     
    11221117    }
    11231118  }
     1119  // put row of tableau in rowArray and columnArray (packed)
     1120  // get subset which have nonzero tableau elements
     1121  transposeTimes2(updates,spareColumn1,alternateWeights_,spareColumn2,spareRow2,
     1122                  -scaleFactor);
     1123  // zero updateBy
     1124  CoinZeroN(updateBy,number);
    11241125  alternateWeights_->clear();
    11251126  // columns
     
    15371538    bool needSubset =  (mode_<4||numberSwitched_>1||mode_>=10);
    15381539
    1539     if (needSubset) {
    1540       // now update weight update array
    1541       model_->factorization()->updateColumnTranspose(spareRow2,
    1542                                                      alternateWeights_);
    1543       transposeTimes2(updates,spareColumn1,alternateWeights_,spareColumn2,spareRow2,0.0);
    1544     } else {
    1545       // put row of tableau in rowArray and columnArray
    1546       model_->clpMatrix()->transposeTimes(model_,-1.0,
    1547                                           updates,spareColumn2,spareColumn1);
    1548     }
    15491540    double * weight;
    15501541    double * other = alternateWeights_->denseVector();
     
    15551546    updateBy = updates->denseVector();
    15561547    weight = weights_+numberColumns;
    1557    
    15581548    if (needSubset) {
    1559       // Exact
     1549      // now update weight update array
     1550      model_->factorization()->updateColumnTranspose(spareRow2,
     1551                                                     alternateWeights_);
     1552      // do alternateWeights_ here so can scale
    15601553      for (j=0;j<number;j++) {
    15611554        int iSequence = index[j];
     
    15631556        // row has -1
    15641557        double pivot = - updateBy[j];
    1565         updateBy[j]=0.0;
    15661558        double modification = other[iSequence];
    15671559        double pivotSquared = pivot * pivot;
     
    15821574        weight[iSequence] = thisWeight;
    15831575      }
     1576      transposeTimes2(updates,spareColumn1,alternateWeights_,spareColumn2,spareRow2,0.0);
     1577    } else {
     1578      // put row of tableau in rowArray and columnArray
     1579      model_->clpMatrix()->transposeTimes(model_,-1.0,
     1580                                          updates,spareColumn2,spareColumn1);
     1581    }
     1582   
     1583    if (needSubset) {
     1584      CoinZeroN(updateBy,number);
    15841585    } else if (mode_==4) {
    15851586      // Devex
  • trunk/Clp/src/ClpQuadraticObjective.cpp

    r800 r1034  
    1111// Constructors / Destructor / Assignment
    1212//#############################################################################
    13 
    1413//-------------------------------------------------------------------
    1514// Default Constructor
     
    129128          printf("number above = %d, number below = %d, error\n",
    130129                 numberAbove,numberBelow);
     130          abort();
    131131        }
    132132      } else {
     
    262262    delete quadraticObjective_;
    263263    quadraticObjective_ = NULL;
     264    delete [] objective_;
     265    delete [] gradient_;
    264266    ClpObjective::operator=(rhs);
    265267    numberColumns_=rhs.numberColumns_;
     
    893895  return CoinMin(theta,maximumTheta);
    894896}
     897// Return objective value (without any ClpModel offset) (model may be NULL)
     898double
     899ClpQuadraticObjective::objectiveValue(const ClpSimplex * model, const double * solution) const
     900{
     901  bool scaling=false;
     902  if (model&&(model->rowScale()||
     903              model->objectiveScale()!=1.0))
     904    scaling=true;
     905  const double * cost = NULL;
     906  if (model)
     907    cost = model->costRegion();
     908  if (!cost) {
     909    // not in solve
     910    cost = objective_;
     911    scaling=false;
     912  }
     913  double linearCost =0.0;
     914  int numberColumns = model->numberColumns();
     915  int numberTotal = numberColumns;
     916  double currentObj=0.0;
     917  for (int iColumn=0;iColumn<numberTotal;iColumn++) {
     918    linearCost += cost[iColumn]*solution[iColumn];
     919  }
     920  if (!activated_||!quadraticObjective_) {
     921    currentObj=linearCost;
     922    return currentObj;
     923  }
     924  assert (model);
     925  const int * columnQuadratic = quadraticObjective_->getIndices();
     926  const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
     927  const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
     928  const double * quadraticElement = quadraticObjective_->getElements();
     929  double c=0.0;
     930  if (!scaling) {
     931    if (!fullMatrix_) {
     932      int iColumn;
     933      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     934        double valueI = solution[iColumn];
     935        CoinBigIndex j;
     936        for (j=columnQuadraticStart[iColumn];
     937             j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     938          int jColumn = columnQuadratic[j];
     939          double valueJ = solution[jColumn];
     940          double elementValue = quadraticElement[j];
     941          if (iColumn!=jColumn) {
     942            c += valueI*valueJ*elementValue;
     943          } else {
     944            c += 0.5*valueI*valueI*elementValue;
     945          }
     946        }
     947      }
     948    } else {
     949      // full matrix stored
     950      int iColumn;
     951      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     952        double valueI = solution[iColumn];
     953        CoinBigIndex j;
     954        for (j=columnQuadraticStart[iColumn];
     955             j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     956          int jColumn = columnQuadratic[j];
     957          double valueJ = solution[jColumn];
     958          double elementValue = quadraticElement[j];
     959          valueJ *= elementValue;
     960          c += valueI*valueJ;
     961        }
     962      }
     963      c *= 0.5;
     964    }
     965  } else {
     966    // scaling
     967    // for now only if half
     968    assert (!fullMatrix_);
     969    const double * columnScale = model->columnScale();
     970    double direction = model->objectiveScale();
     971    // direction is actually scale out not scale in
     972    if (direction)
     973      direction = 1.0/direction;
     974    if (!columnScale) {
     975      for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
     976        double valueI = solution[iColumn];
     977        CoinBigIndex j;
     978        for (j=columnQuadraticStart[iColumn];
     979             j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     980          int jColumn = columnQuadratic[j];
     981          double valueJ = solution[jColumn];
     982          double elementValue = quadraticElement[j];
     983          elementValue *= direction;
     984          if (iColumn!=jColumn) {
     985            c += valueI*valueJ*elementValue;
     986          } else {
     987            c += 0.5*valueI*valueI*elementValue;
     988          }
     989        }
     990      }
     991    } else {
     992      // scaling
     993      for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
     994        double valueI = solution[iColumn];
     995        double scaleI = columnScale[iColumn]*direction;
     996        CoinBigIndex j;
     997        for (j=columnQuadraticStart[iColumn];
     998             j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     999          int jColumn = columnQuadratic[j];
     1000          double valueJ = solution[jColumn];
     1001          double elementValue = quadraticElement[j];
     1002          elementValue *= scaleI*columnScale[jColumn];
     1003          if (iColumn!=jColumn) {
     1004            c += valueI*valueJ*elementValue;
     1005          } else {
     1006            c += 0.5*valueI*valueI*elementValue;
     1007          }
     1008        }
     1009      }
     1010    }
     1011  }
     1012  currentObj = c+linearCost;
     1013  return currentObj;
     1014}
    8951015// Scale objective
    8961016void
  • trunk/Clp/src/ClpQuadraticObjective.hpp

    r754 r1034  
    4747                            double & predictedObj,
    4848                            double & thetaObj);
     49  /// Return objective value (without any ClpModel offset) (model may be NULL)
     50  virtual double objectiveValue(const ClpSimplex * model, const double * solution) const ;
    4951  virtual void resize(int newNumberColumns) ;
    5052  /// Delete columns in  objective
  • trunk/Clp/src/ClpSimplex.cpp

    r995 r1034  
    2828#include "ClpHelperFunctions.hpp"
    2929#include "CoinModel.hpp"
     30#include "CoinLpIO.hpp"
    3031#include <cfloat>
    3132
     
    3536//#############################################################################
    3637
    37 ClpSimplex::ClpSimplex () :
     38ClpSimplex::ClpSimplex (bool emptyMessages) :
    3839
    39   ClpModel(),
     40  ClpModel(emptyMessages),
    4041  columnPrimalInfeasibility_(0.0),
    4142  rowPrimalInfeasibility_(0.0),
     
    5152  largestPrimalError_(0.0),
    5253  largestDualError_(0.0),
    53   largestSolutionError_(0.0),
     54  alphaAccuracy_(-1.0),
    5455  dualBound_(1.0e10),
    5556  alpha_(0.0),
     
    101102  savedSolution_(NULL),
    102103  numberTimesOptimal_(0),
     104  disasterArea_(NULL),
    103105  changeMade_(1),
    104106  algorithm_(0),
     
    106108  perturbation_(100),
    107109  nonLinearCost_(NULL),
    108   specialOptions_(0),
    109110  lastBadIteration_(-999999),
    110111  lastFlaggedIteration_(-999999),
     
    131132  saveStatus_=NULL;
    132133  // get an empty factorization so we can set tolerances etc
    133   factorization_ = new ClpFactorization();
     134  getEmptyFactorization();
    134135  // Say sparse
    135136  factorization_->sparseThreshold(1);
     
    162163  largestPrimalError_(0.0),
    163164  largestDualError_(0.0),
    164   largestSolutionError_(0.0),
     165  alphaAccuracy_(-1.0),
    165166  dualBound_(1.0e10),
    166167  alpha_(0.0),
     
    212213  savedSolution_(NULL),
    213214  numberTimesOptimal_(0),
     215  disasterArea_(NULL),
    214216  changeMade_(1),
    215217  algorithm_(0),
     
    217219  perturbation_(100),
    218220  nonLinearCost_(NULL),
    219   specialOptions_(0),
    220221  lastBadIteration_(-999999),
    221222  lastFlaggedIteration_(-999999),
     
    242243  saveStatus_=NULL;
    243244  // get an empty factorization so we can set tolerances etc
    244   factorization_ = new ClpFactorization();
     245  getEmptyFactorization();
    245246  // say Steepest pricing
    246247  dualRowPivot_ = new ClpDualRowSteepest();
     
    327328ClpSimplex::~ClpSimplex ()
    328329{
     330  setPersistenceFlag(0);
    329331  gutsOfDelete(0);
    330332  delete nonLinearCost_;
     
    430432      if (!numberBasic) {
    431433        //printf("no errors on basic - going to all slack - numberOut %d\n",numberOut);
    432         allSlackBasis();
     434        allSlackBasis(true);
    433435      }
    434436      CoinSort_2(save, save + numberOut, sort,
     
    14591461    progressFlag_ |= 1; // making real progress
    14601462  if (sequenceIn_!=sequenceOut_) {
     1463    if (alphaAccuracy_>0.0) {
     1464      double value = fabs(alpha_);
     1465      if (value>1.0)
     1466        alphaAccuracy_ *= value;
     1467      else
     1468        alphaAccuracy_ /= value;
     1469    }
    14611470    //assert( getStatus(sequenceOut_)== basic);
    14621471    setStatus(sequenceIn_,basic);
     
    15151524  int in = sequenceIn_;
    15161525  int out = sequenceOut_;
    1517   matrix_->correctSequence(in,out);
     1526  matrix_->correctSequence(this,in,out);
    15181527  int cycle=progress_->cycle(in,out,
    15191528                            directionIn_,directionOut_);
     
    16031612  largestPrimalError_(0.0),
    16041613  largestDualError_(0.0),
    1605   largestSolutionError_(0.0),
     1614  alphaAccuracy_(-1.0),
    16061615  dualBound_(1.0e10),
    16071616  alpha_(0.0),
     
    16531662  savedSolution_(NULL),
    16541663  numberTimesOptimal_(0),
     1664  disasterArea_(NULL),
    16551665  changeMade_(1),
    16561666  algorithm_(0),
     
    16581668  perturbation_(100),
    16591669  nonLinearCost_(NULL),
    1660   specialOptions_(0),
    16611670  lastBadIteration_(-999999),
    16621671  lastFlaggedIteration_(-999999),
     
    16861695  primalColumnPivot_ = NULL;
    16871696  gutsOfDelete(0);
    1688   specialOptions_ =0;
    16891697  delete nonLinearCost_;
    16901698  nonLinearCost_ = NULL;
     
    17081716  largestPrimalError_(0.0),
    17091717  largestDualError_(0.0),
    1710   largestSolutionError_(0.0),
     1718  alphaAccuracy_(-1.0),
    17111719  dualBound_(1.0e10),
    17121720  alpha_(0.0),
     
    17581766  savedSolution_(NULL),
    17591767  numberTimesOptimal_(0),
     1768  disasterArea_(NULL),
    17601769  changeMade_(1),
    17611770  algorithm_(0),
     
    17631772  perturbation_(100),
    17641773  nonLinearCost_(NULL),
    1765   specialOptions_(0),
    17661774  lastBadIteration_(-999999),
    17671775  lastFlaggedIteration_(-999999),
     
    17881796  saveStatus_=NULL;
    17891797  // get an empty factorization so we can set tolerances etc
    1790   factorization_ = new ClpFactorization();
     1798  getEmptyFactorization();
    17911799  // say Steepest pricing
    17921800  dualRowPivot_ = new ClpDualRowSteepest();
     
    18021810  if (this != &rhs) {
    18031811    gutsOfDelete(0);
    1804     specialOptions_=0;
    18051812    delete nonLinearCost_;
    18061813    nonLinearCost_ = NULL;
     
    18851892  }
    18861893  if (rhs.factorization_) {
     1894    delete factorization_;
    18871895    factorization_ = new ClpFactorization(*rhs.factorization_);
    18881896  } else {
     
    19021910  largestPrimalError_ = rhs.largestPrimalError_;
    19031911  largestDualError_ = rhs.largestDualError_;
    1904   largestSolutionError_ = rhs.largestSolutionError_;
     1912  alphaAccuracy_ = rhs.alphaAccuracy_;
    19051913  dualBound_ = rhs.dualBound_;
    19061914  alpha_ = rhs.alpha_;
     
    19321940  primalColumnPivot_ = rhs.primalColumnPivot_->clone(true);
    19331941  numberTimesOptimal_ = rhs.numberTimesOptimal_;
     1942  disasterArea_ = NULL;
    19341943  changeMade_ = rhs.changeMade_;
    19351944  algorithm_ = rhs.algorithm_;
     
    19371946  perturbation_ = rhs.perturbation_;
    19381947  infeasibilityCost_ = rhs.infeasibilityCost_;
    1939   specialOptions_ = rhs.specialOptions_;
    19401948  lastBadIteration_ = rhs.lastBadIteration_;
    19411949  lastFlaggedIteration_ = rhs.lastFlaggedIteration_;
     
    19952003  }
    19962004  int i;
    1997   for (i=0;i<6;i++) {
    1998     delete rowArray_[i];
    1999     rowArray_[i]=NULL;
    2000     delete columnArray_[i];
    2001     columnArray_[i]=NULL;
     2005  if ((specialOptions_&65536)==0) {
     2006    for (i=0;i<6;i++) {
     2007      delete rowArray_[i];
     2008      rowArray_[i]=NULL;
     2009      delete columnArray_[i];
     2010      columnArray_[i]=NULL;
     2011    }
    20022012  }
    20032013  delete rowCopy_;
     
    20092019    delete auxiliaryModel_;
    20102020    auxiliaryModel_ = NULL;
    2011     delete factorization_;
    2012     factorization_ = NULL;
     2021    setEmptyFactorization();
    20132022    delete [] pivotVariable_;
    20142023    pivotVariable_=NULL;
     
    20402049  objectiveValue_ = 0.0;
    20412050  // now look at primal solution
    2042   columnPrimalInfeasibility_=0.0;
    2043   columnPrimalSequence_=-1;
    2044   rowPrimalInfeasibility_=0.0;
    2045   rowPrimalSequence_=-1;
    2046   largestSolutionError_=0.0;
    20472051  solution = rowActivityWork_;
    20482052  sumPrimalInfeasibilities_=0.0;
     
    20552059  relaxedTolerance = relaxedTolerance +  error;
    20562060  sumOfRelaxedPrimalInfeasibilities_ = 0.0;
    2057 
    20582061  for (iRow=0;iRow<numberRows_;iRow++) {
    20592062    //assert (fabs(solution[iRow])<1.0e15||getRowStatus(iRow) == basic);
     
    20712074      numberPrimalInfeasibilities_ ++;
    20722075    }
    2073     if (infeasibility>rowPrimalInfeasibility_) {
    2074       rowPrimalInfeasibility_=infeasibility;
    2075       rowPrimalSequence_=iRow;
    2076     }
    20772076    infeasibility = fabs(rowActivities[iRow]-solution[iRow]);
    2078     if (infeasibility>largestSolutionError_)
    2079       largestSolutionError_=infeasibility;
    20802077  }
    20812078  // Check any infeasibilities from dynamic rows
     
    20922089        infeasibility=columnLowerWork_[iColumn]-solution[iColumn];
    20932090      }
    2094       if (infeasibility>columnPrimalInfeasibility_) {
    2095         columnPrimalInfeasibility_=infeasibility;
    2096         columnPrimalSequence_=iColumn;
    2097       }
    20982091      if (infeasibility>primalTolerance) {
    20992092        sumPrimalInfeasibilities_ += infeasibility-primalTolerance_;
     
    21032096      }
    21042097      infeasibility = fabs(columnActivities[iColumn]-solution[iColumn]);
    2105       if (infeasibility>largestSolutionError_)
    2106         largestSolutionError_=infeasibility;
    21072098    }
    21082099  } else {
     
    21192110        infeasibility=columnLowerWork_[iColumn]-solution[iColumn];
    21202111      }
    2121       if (infeasibility>columnPrimalInfeasibility_) {
    2122         columnPrimalInfeasibility_=infeasibility;
    2123         columnPrimalSequence_=iColumn;
    2124       }
    21252112      if (infeasibility>primalTolerance) {
    21262113        sumPrimalInfeasibilities_ += infeasibility-primalTolerance_;
     
    21302117      }
    21312118      infeasibility = fabs(columnActivities[iColumn]-solution[iColumn]);
    2132       if (infeasibility>largestSolutionError_)
    2133         largestSolutionError_=infeasibility;
    21342119    }
    21352120  }
     
    21452130  numberDualInfeasibilities_=0;
    21462131  numberDualInfeasibilitiesWithoutFree_=0;
    2147   columnDualInfeasibility_=0.0;
    2148   columnDualSequence_=-1;
    21492132  if (matrix_->skipDualCheck()&&algorithm_>0&&problemStatus_==-2) {
    21502133    // pretend we found dual infeasibilities
     
    21542137    return;
    21552138  }
    2156   rowDualInfeasibility_=0.0;
    2157   rowDualSequence_=-1;
    21582139  int firstFreePrimal = -1;
    21592140  int firstFreeDual = -1;
    21602141  int numberSuperBasicWithDj=0;
    2161   primalToleranceToGetOptimal_=CoinMax(rowPrimalInfeasibility_,
    2162                                    columnPrimalInfeasibility_);
    2163   remainingDualInfeasibility_=0.0;
    21642142  double relaxedTolerance=dualTolerance_;
    21652143  // we can't really trust infeasibilities if there is dual error
     
    21942172        if (value<0.0) {
    21952173          value = - value;
    2196           if (value>columnDualInfeasibility_) {
    2197             columnDualInfeasibility_=value;
    2198             columnDualSequence_=iColumn;
    2199           }
    22002174          if (value>dualTolerance_) {
    22012175            if (getColumnStatus(iColumn) != isFree) {
     
    22152189              }
    22162190            }
    2217             // maybe we can make feasible by increasing tolerance
    2218             if (distanceUp<largeValue_) {
    2219               if (distanceUp>primalToleranceToGetOptimal_)
    2220                 primalToleranceToGetOptimal_=distanceUp;
    2221             } else {
    2222               //gap too big for any tolerance
    2223               remainingDualInfeasibility_=
    2224                 CoinMax(remainingDualInfeasibility_,value);
    2225             }
    22262191          }
    22272192        }
     
    22312196        // should not be positive
    22322197        if (value>0.0) {
    2233           if (value>columnDualInfeasibility_) {
    2234             columnDualInfeasibility_=value;
    2235             columnDualSequence_=iColumn;
    2236           }
    22372198          if (value>dualTolerance_) {
    22382199            sumDualInfeasibilities_ += value-dualTolerance_;
     
    22432204              numberDualInfeasibilitiesWithoutFree_ ++;
    22442205            // maybe we can make feasible by increasing tolerance
    2245             if (distanceDown<largeValue_&&
    2246                 distanceDown>primalToleranceToGetOptimal_)
    2247               primalToleranceToGetOptimal_=distanceDown;
    22482206          }
    22492207        }
     
    22712229        if (value<0.0) {
    22722230          value = - value;
    2273           if (value>rowDualInfeasibility_) {
    2274             rowDualInfeasibility_=value;
    2275             rowDualSequence_=iRow;
     2231          if (value>dualTolerance_) {
     2232            sumDualInfeasibilities_ += value-dualTolerance_;
     2233            if (value>relaxedTolerance)
     2234              sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
     2235            numberDualInfeasibilities_ ++;
     2236            if (getRowStatus(iRow) != isFree)
     2237              numberDualInfeasibilitiesWithoutFree_ ++;
    22762238          }
     2239        }
     2240      }
     2241      if (distanceDown>primalTolerance_) {
     2242        double value = rowReducedCost_[iRow];
     2243        // should not be positive
     2244        if (value>0.0) {
    22772245          if (value>dualTolerance_) {
    22782246            sumDualInfeasibilities_ += value-dualTolerance_;
     
    22832251              numberDualInfeasibilitiesWithoutFree_ ++;
    22842252            // maybe we can make feasible by increasing tolerance
    2285             if (distanceUp<largeValue_) {
    2286               if (distanceUp>primalToleranceToGetOptimal_)
    2287                 primalToleranceToGetOptimal_=distanceUp;
    2288             } else {
    2289               //gap too big for any tolerance
    2290               remainingDualInfeasibility_=
    2291                 CoinMax(remainingDualInfeasibility_,value);
    2292             }
    2293           }
    2294         }
    2295       }
    2296       if (distanceDown>primalTolerance_) {
    2297         double value = rowReducedCost_[iRow];
    2298         // should not be positive
    2299         if (value>0.0) {
    2300           if (value>rowDualInfeasibility_) {
    2301             rowDualInfeasibility_=value;
    2302             rowDualSequence_=iRow;
    2303           }
    2304           if (value>dualTolerance_) {
    2305             sumDualInfeasibilities_ += value-dualTolerance_;
    2306             if (value>relaxedTolerance)
    2307               sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
    2308             numberDualInfeasibilities_ ++;
    2309             if (getRowStatus(iRow) != isFree)
    2310               numberDualInfeasibilitiesWithoutFree_ ++;
    2311             // maybe we can make feasible by increasing tolerance
    2312             if (distanceDown<largeValue_&&
    2313                 distanceDown>primalToleranceToGetOptimal_)
    2314               primalToleranceToGetOptimal_=distanceDown;
    23152253          }
    23162254        }
     
    25202458  }
    25212459}
     2460//static int scale_times[]={0,0,0,0};
    25222461bool
    25232462ClpSimplex::createRim(int what,bool makeRowCopy, int startFinishOptions)
     
    25392478  if (auxiliaryModel_) {
    25402479    if (auxiliaryModel_->numberRows_==numberRows_&&
    2541         auxiliaryModel_->numberColumns_==numberColumns_)
     2480        auxiliaryModel_->numberColumns_==numberColumns_&&
     2481        (whatsChanged_&511)==511) {
    25422482      oldMatrix=true;
    2543     else
     2483    } else {
    25442484      deleteAuxiliaryModel();
     2485      oldMatrix=false;
     2486    }
    25452487  }
    25462488  if (what==63) {
     
    26042546    if (oldMatrix)
    26052547      checkType = 14;
     2548    if ((specialOptions_&0x1000000)!=0)
     2549      checkType -= 4; // don't check for duplicates
    26062550    if (!matrix_->allElementsInRange(this,smallElement_,1.0e20,checkType)) {
    26072551      problemStatus_=4;
     2552      secondaryStatus_=8;
    26082553      //goodMatrix= false;
    26092554      return false;
     
    26142559      rowCopy_ = matrix_->reverseOrderedCopy();
    26152560    }
     2561#if 0
     2562    if (what==63&&(specialOptions_&131072)!=0) {
     2563      int k=rowScale_ ? 1 : 0;
     2564      if (oldMatrix)
     2565        k+=2;
     2566      scale_times[k]++;
     2567      if ((scale_times[0]+scale_times[1]+scale_times[2]+scale_times[3])%1000==0)
     2568        printf("scale counts %d %d %d %d\n",
     2569               scale_times[0],scale_times[1],scale_times[2],scale_times[3]);
     2570    }
     2571#endif
    26162572    // do scaling if needed
    2617     if (!oldMatrix) {
     2573    if (!oldMatrix&&scalingFlag_<0) {
    26182574      if (scalingFlag_<0&&rowScale_) {
    2619         if (handler_->logLevel()>0)
     2575        //if (handler_->logLevel()>0)
    26202576          printf("How did we get scalingFlag_ %d and non NULL rowScale_? - switching off scaling\n",
    26212577                 scalingFlag_);
     
    27532709      if (clpMatrix&&numberThreads_)
    27542710        clpMatrix->specialRowCopy(this,rowCopy_);
     2711      if (clpMatrix)
     2712        clpMatrix->specialColumnCopy(this);
    27552713    }
    27562714  }
     
    28462804        // If scaled then do all columns later in one loop
    28472805        if (what!=63) {
    2848           for (i=0;i<numberColumns_;i++) {
    2849             double multiplier = rhsScale_/columnScale[i];
    2850             double lowerValue = columnLower_[i];
    2851             double upperValue = columnUpper_[i];
    2852             if (lowerValue>-1.0e20) {
    2853               columnLowerWork_[i]=lowerValue*multiplier;
    2854               if (upperValue>=1.0e20) {
    2855                 columnUpperWork_[i]=COIN_DBL_MAX;
    2856               } else {
    2857                 columnUpperWork_[i]=upperValue*multiplier;
    2858                 if (fabs(columnUpperWork_[i]-columnLowerWork_[i])<=primalTolerance) {
    2859                   if (columnLowerWork_[i]>=0.0) {
    2860                     columnUpperWork_[i] = columnLowerWork_[i];
    2861                   } else if (columnUpperWork_[i]<=0.0) {
    2862                     columnLowerWork_[i] = columnUpperWork_[i];
    2863                   } else {
    2864                     columnUpperWork_[i] = 0.0;
    2865                     columnLowerWork_[i] = 0.0;
    2866                   }
    2867                 }
    2868               }
    2869             } else if (upperValue<1.0e20) {
    2870               columnLowerWork_[i]=-COIN_DBL_MAX;
    2871               columnUpperWork_[i]=upperValue*multiplier;
    2872             } else {
    2873               // free
    2874               columnLowerWork_[i]=-COIN_DBL_MAX;
    2875               columnUpperWork_[i]=COIN_DBL_MAX;
    2876             }
    2877           }
    2878         }
     2806          if ((specialOptions_&131072)==0) {
     2807            for (i=0;i<numberColumns_;i++) {
     2808              double multiplier = rhsScale_/columnScale[i];
     2809              double lowerValue = columnLower_[i];
     2810              double upperValue = columnUpper_[i];
     2811              if (lowerValue>-1.0e20) {
     2812                columnLowerWork_[i]=lowerValue*multiplier;
     2813                if (upperValue>=1.0e20) {
     2814                  columnUpperWork_[i]=COIN_DBL_MAX;
     2815                } else {
     2816                  columnUpperWork_[i]=upperValue*multiplier;
     2817                  if (fabs(columnUpperWork_[i]-columnLowerWork_[i])<=primalTolerance) {
     2818                    if (columnLowerWork_[i]>=0.0) {
     2819                      columnUpperWork_[i] = columnLowerWork_[i];
     2820                    } else if (columnUpperWork_[i]<=0.0) {
     2821                      columnLowerWork_[i] = columnUpperWork_[i];
     2822                    } else {
     2823                      columnUpperWork_[i] = 0.0;
     2824                      columnLowerWork_[i] = 0.0;
     2825                    }
     2826                  }
     2827                }
     2828              } else if (upperValue<1.0e20) {
     2829                columnLowerWork_[i]=-COIN_DBL_MAX;
     2830                columnUpperWork_[i]=upperValue*multiplier;
     2831              } else {
     2832                // free
     2833                columnLowerWork_[i]=-COIN_DBL_MAX;
     2834                columnUpperWork_[i]=COIN_DBL_MAX;
     2835              }
     2836            }
     2837          } else {
     2838            const double * inverseScale = columnScale_+numberColumns_;
     2839            for (i=0;i<numberColumns_;i++) {
     2840              double multiplier = rhsScale_*inverseScale[i];
     2841              double lowerValue = columnLower_[i];
     2842              double upperValue = columnUpper_[i];
     2843              if (lowerValue>-1.0e20) {
     2844                columnLowerWork_[i]=lowerValue*multiplier;
     2845                if (upperValue>=1.0e20) {
     2846                  columnUpperWork_[i]=COIN_DBL_MAX;
     2847                } else {
     2848                  columnUpperWork_[i]=upperValue*multiplier;
     2849                  if (fabs(columnUpperWork_[i]-columnLowerWork_[i])<=primalTolerance) {
     2850                    if (columnLowerWork_[i]>=0.0) {
     2851                      columnUpperWork_[i] = columnLowerWork_[i];
     2852                    } else if (columnUpperWork_[i]<=0.0) {
     2853                      columnLowerWork_[i] = columnUpperWork_[i];
     2854                    } else {
     2855                      columnUpperWork_[i] = 0.0;
     2856                      columnLowerWork_[i] = 0.0;
     2857                    }
     2858                  }
     2859                }
     2860              } else if (upperValue<1.0e20) {
     2861                columnLowerWork_[i]=-COIN_DBL_MAX;
     2862                columnUpperWork_[i]=upperValue*multiplier;
     2863              } else {
     2864                // free
     2865                columnLowerWork_[i]=-COIN_DBL_MAX;
     2866                columnUpperWork_[i]=COIN_DBL_MAX;
     2867              }
     2868            }
     2869          }
     2870        }
    28792871        for (i=0;i<numberRows_;i++) {
    28802872          double multiplier = rhsScale_*rowScale[i];
     
    31633155        double primalTolerance=dblParam_[ClpPrimalTolerance];
    31643156        // on entry
    3165         for (i=0;i<numberColumns_;i++) {
    3166           CoinAssert(fabs(obj[i])<1.0e25);
    3167           double scaleFactor = columnScale_[i];
    3168           double multiplier = rhsScale_/scaleFactor;
    3169           scaleFactor *= direction;
    3170           objectiveWork_[i] = obj[i]*scaleFactor;
    3171           reducedCostWork_[i] = reducedCost_[i]*scaleFactor;
    3172           double lowerValue = columnLower_[i];
    3173           double upperValue = columnUpper_[i];
    3174           if (lowerValue>-1.0e20) {
    3175             columnLowerWork_[i]=lowerValue*multiplier;
    3176             if (upperValue>=1.0e20) {
    3177               columnUpperWork_[i]=COIN_DBL_MAX;
    3178             } else {
    3179               columnUpperWork_[i]=upperValue*multiplier;
    3180               if (fabs(columnUpperWork_[i]-columnLowerWork_[i])<=primalTolerance) {
    3181                 if (columnLowerWork_[i]>=0.0) {
    3182                   columnUpperWork_[i] = columnLowerWork_[i];
    3183                 } else if (columnUpperWork_[i]<=0.0) {
    3184                   columnLowerWork_[i] = columnUpperWork_[i];
    3185                 } else {
    3186                   columnUpperWork_[i] = 0.0;
    3187                   columnLowerWork_[i] = 0.0;
    3188                 }
    3189               }
    3190             }
    3191           } else if (upperValue<1.0e20) {
    3192             columnLowerWork_[i]=-COIN_DBL_MAX;
    3193             columnUpperWork_[i]=upperValue*multiplier;
    3194           } else {
    3195             // free
    3196             columnLowerWork_[i]=-COIN_DBL_MAX;
    3197             columnUpperWork_[i]=COIN_DBL_MAX;
    3198           }
    3199           double value = columnActivity_[i] * multiplier;
    3200           if (fabs(value)>1.0e20) {
    3201             //printf("bad value of %g for column %d\n",value,i);
    3202             setColumnStatus(i,superBasic);
    3203             if (columnUpperWork_[i]<0.0) {
    3204               value=columnUpperWork_[i];
    3205             } else if (columnLowerWork_[i]>0.0) {
    3206               value=columnLowerWork_[i];
    3207             } else {
    3208               value=0.0;
    3209             }
    3210           }
    3211           columnActivityWork_[i] = value;
    3212         }
    3213         for (i=0;i<numberRows_;i++) {
    3214           dual_[i] /= rowScale_[i];
    3215           dual_[i] *= objectiveScale_;
    3216           rowReducedCost_[i] = dual_[i];
    3217           double multiplier = rhsScale_*rowScale_[i];
    3218           double value = rowActivity_[i] * multiplier;
    3219           if (fabs(value)>1.0e20) {
    3220             //printf("bad value of %g for row %d\n",value,i);
    3221             setRowStatus(i,superBasic);
    3222             if (rowUpperWork_[i]<0.0) {
    3223               value=rowUpperWork_[i];
    3224             } else if (rowLowerWork_[i]>0.0) {
    3225               value=rowLowerWork_[i];
    3226             } else {
    3227               value=0.0;
    3228             }
    3229           }
    3230           rowActivityWork_[i] = value;
    3231         }
     3157        if ((specialOptions_&131072)==0) {
     3158          for (i=0;i<numberColumns_;i++) {
     3159            CoinAssert(fabs(obj[i])<1.0e25);
     3160            double scaleFactor = columnScale_[i];
     3161            double multiplier = rhsScale_/scaleFactor;
     3162            scaleFactor *= direction;
     3163            objectiveWork_[i] = obj[i]*scaleFactor;
     3164            reducedCostWork_[i] = reducedCost_[i]*scaleFactor;
     3165            double lowerValue = columnLower_[i];
     3166            double upperValue = columnUpper_[i];
     3167            if (lowerValue>-1.0e20) {
     3168              columnLowerWork_[i]=lowerValue*multiplier;
     3169              if (upperValue>=1.0e20) {
     3170                columnUpperWork_[i]=COIN_DBL_MAX;
     3171              } else {
     3172                columnUpperWork_[i]=upperValue*multiplier;
     3173                if (fabs(columnUpperWork_[i]-columnLowerWork_[i])<=primalTolerance) {
     3174                  if (columnLowerWork_[i]>=0.0) {
     3175                    columnUpperWork_[i] = columnLowerWork_[i];
     3176                  } else if (columnUpperWork_[i]<=0.0) {
     3177                    columnLowerWork_[i] = columnUpperWork_[i];
     3178                  } else {
     3179                    columnUpperWork_[i] = 0.0;
     3180                    columnLowerWork_[i] = 0.0;
     3181                  }
     3182                }
     3183              }
     3184            } else if (upperValue<1.0e20) {
     3185              columnLowerWork_[i]=-COIN_DBL_MAX;
     3186              columnUpperWork_[i]=upperValue*multiplier;
     3187            } else {
     3188              // free
     3189              columnLowerWork_[i]=-COIN_DBL_MAX;
     3190              columnUpperWork_[i]=COIN_DBL_MAX;
     3191            }
     3192            double value = columnActivity_[i] * multiplier;
     3193            if (fabs(value)>1.0e20) {
     3194              //printf("bad value of %g for column %d\n",value,i);
     3195              setColumnStatus(i,superBasic);
     3196              if (columnUpperWork_[i]<0.0) {
     3197                value=columnUpperWork_[i];
     3198              } else if (columnLowerWork_[i]>0.0) {
     3199                value=columnLowerWork_[i];
     3200              } else {
     3201                value=0.0;
     3202              }
     3203            }
     3204            columnActivityWork_[i] = value;
     3205          }
     3206          for (i=0;i<numberRows_;i++) {
     3207            dual_[i] /= rowScale_[i];
     3208            dual_[i] *= objectiveScale_;
     3209            rowReducedCost_[i] = dual_[i];
     3210            double multiplier = rhsScale_*rowScale_[i];
     3211            double value = rowActivity_[i] * multiplier;
     3212            if (fabs(value)>1.0e20) {
     3213              //printf("bad value of %g for row %d\n",value,i);
     3214              setRowStatus(i,superBasic);
     3215              if (rowUpperWork_[i]<0.0) {
     3216                value=rowUpperWork_[i];
     3217              } else if (rowLowerWork_[i]>0.0) {
     3218                value=rowLowerWork_[i];
     3219              } else {
     3220                value=0.0;
     3221              }
     3222            }
     3223            rowActivityWork_[i] = value;
     3224          }
     3225        } else {
     3226          const double * inverseScale = columnScale_+numberColumns_;
     3227          for (i=0;i<numberColumns_;i++) {
     3228            CoinAssert(fabs(obj[i])<1.0e25);
     3229            double scaleFactor = columnScale_[i];
     3230            double multiplier = rhsScale_*inverseScale[i];
     3231            scaleFactor *= direction;
     3232            objectiveWork_[i] = obj[i]*scaleFactor;
     3233            reducedCostWork_[i] = reducedCost_[i]*scaleFactor;
     3234            double lowerValue = columnLower_[i];
     3235            double upperValue = columnUpper_[i];
     3236            if (lowerValue>-1.0e20) {
     3237              columnLowerWork_[i]=lowerValue*multiplier;
     3238              if (upperValue>=1.0e20) {
     3239                columnUpperWork_[i]=COIN_DBL_MAX;
     3240              } else {
     3241                columnUpperWork_[i]=upperValue*multiplier;
     3242                if (fabs(columnUpperWork_[i]-columnLowerWork_[i])<=primalTolerance) {
     3243                  if (columnLowerWork_[i]>=0.0) {
     3244                    columnUpperWork_[i] = columnLowerWork_[i];
     3245                  } else if (columnUpperWork_[i]<=0.0) {
     3246                    columnLowerWork_[i] = columnUpperWork_[i];
     3247                  } else {
     3248                    columnUpperWork_[i] = 0.0;
     3249                    columnLowerWork_[i] = 0.0;
     3250                  }
     3251                }
     3252              }
     3253            } else if (upperValue<1.0e20) {
     3254              columnLowerWork_[i]=-COIN_DBL_MAX;
     3255              columnUpperWork_[i]=upperValue*multiplier;
     3256            } else {
     3257              // free
     3258              columnLowerWork_[i]=-COIN_DBL_MAX;
     3259              columnUpperWork_[i]=COIN_DBL_MAX;
     3260            }
     3261            double value = columnActivity_[i] * multiplier;
     3262            if (fabs(value)>1.0e20) {
     3263              //printf("bad value of %g for column %d\n",value,i);
     3264              setColumnStatus(i,superBasic);
     3265              if (columnUpperWork_[i]<0.0) {
     3266                value=columnUpperWork_[i];
     3267              } else if (columnLowerWork_[i]>0.0) {
     3268                value=columnLowerWork_[i];
     3269              } else {
     3270                value=0.0;
     3271              }
     3272            }
     3273            columnActivityWork_[i] = value;
     3274          }
     3275          inverseScale = rowScale_+numberRows_;
     3276          for (i=0;i<numberRows_;i++) {
     3277            dual_[i] *= inverseScale[i];
     3278            dual_[i] *= objectiveScale_;
     3279            rowReducedCost_[i] = dual_[i];
     3280            double multiplier = rhsScale_*rowScale_[i];
     3281            double value = rowActivity_[i] * multiplier;
     3282            if (fabs(value)>1.0e20) {
     3283              //printf("bad value of %g for row %d\n",value,i);
     3284              setRowStatus(i,superBasic);
     3285              if (rowUpperWork_[i]<0.0) {
     3286                value=rowUpperWork_[i];
     3287              } else if (rowLowerWork_[i]>0.0) {
     3288                value=rowLowerWork_[i];
     3289              } else {
     3290                value=0.0;
     3291              }
     3292            }
     3293            rowActivityWork_[i] = value;
     3294          }
     3295        }
    32323296      } else if (objectiveScale_!=1.0||rhsScale_!=1.0) {
    32333297        // on entry
     
    33823446      *********************************************************/
    33833447      for (iRow=0;iRow<4;iRow++) {
    3384         delete rowArray_[iRow];
    3385         rowArray_[iRow]=new CoinIndexedVector();
    33863448        int length =numberRows2+factorization_->maximumPivots();
    33873449        if (iRow==3||objective_->type()>1)
    33883450          length += numberColumns_;
     3451        if ((specialOptions_&65536)==0||!rowArray_[iRow]) {
     3452          delete rowArray_[iRow];
     3453          rowArray_[iRow]=new CoinIndexedVector();
     3454        }
    33893455        rowArray_[iRow]->reserve(length);
    33903456      }
    33913457     
    33923458      for (iColumn=0;iColumn<2;iColumn++) {
    3393         delete columnArray_[iColumn];
    3394         columnArray_[iColumn]=new CoinIndexedVector();
     3459        if ((specialOptions_&65536)==0||!columnArray_[iColumn]) {
     3460          delete columnArray_[iColumn];
     3461          columnArray_[iColumn]=new CoinIndexedVector();
     3462        }
    33953463        if (!iColumn)
    33963464          columnArray_[iColumn]->reserve(numberColumns_);
     
    34253493        if (iRow==3||objective_->type()>1)
    34263494          length += numberColumns_;
    3427         assert(rowArray_[iRow]->capacity()==length);
     3495        assert(rowArray_[iRow]->capacity()>=length);
    34283496        rowArray_[iRow]->checkClear();
    34293497#endif
     
    34363504        if (iColumn)
    34373505          length=CoinMax(numberRows2,numberColumns_);
    3438         assert(columnArray_[iColumn]->capacity()==length);
     3506        assert(columnArray_[iColumn]->capacity()>=length);
    34393507        columnArray_[iColumn]->checkClear();
    34403508#endif
     
    34533521ClpSimplex::deleteRim(int getRidOfFactorizationData)
    34543522{
     3523  // Just possible empty problem
     3524  int numberRows=numberRows_;
     3525  int numberColumns=numberColumns_;
     3526  if (!numberRows||!numberColumns) {
     3527    numberRows=0;
     3528    if (objective_->type()<2)
     3529      numberColumns=0;
     3530  }
    34553531  if (!auxiliaryModel_) {
    34563532    int i;
     
    34643540    {
    34653541      int nBad=0;
    3466       for (i=0;i<numberColumns_;i++) {
     3542      for (i=0;i<numberColumns;i++) {
    34673543        if (lower_[i]==upper_[i]&&getColumnStatus(i)==basic)
    34683544          nBad++;
     
    34813557      double scaleC = 1.0/objectiveScale_;
    34823558      double scaleR = 1.0/rhsScale_;
    3483       for (i=0;i<numberColumns_;i++) {
    3484         double scaleFactor = columnScale_[i];
    3485         double valueScaled = columnActivityWork_[i];
    3486         double lowerScaled = columnLowerWork_[i];
    3487         double upperScaled = columnUpperWork_[i];
    3488         if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
    3489           if (valueScaled<lowerScaled-primalTolerance_||
    3490               valueScaled>upperScaled+primalTolerance_)
    3491             numberPrimalScaled++;
    3492           else
    3493             upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
    3494         }
    3495         columnActivity_[i] = valueScaled*scaleFactor*scaleR;
    3496         double value = columnActivity_[i];
    3497         if (value<columnLower_[i]-primalTolerance_)
    3498           numberPrimalUnscaled++;
    3499         else if (value>columnUpper_[i]+primalTolerance_)
    3500           numberPrimalUnscaled++;
    3501         double valueScaledDual = reducedCostWork_[i];
    3502         if (valueScaled>columnLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
    3503           numberDualScaled++;
    3504         if (valueScaled<columnUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
    3505           numberDualScaled++;
    3506         reducedCost_[i] = (valueScaledDual*scaleC)/scaleFactor;
    3507         double valueDual = reducedCost_[i];
    3508         if (value>columnLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
    3509           numberDualUnscaled++;
    3510         if (value<columnUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
    3511           numberDualUnscaled++;
    3512       }
    3513       for (i=0;i<numberRows_;i++) {
    3514         double scaleFactor = rowScale_[i];
    3515         double valueScaled = rowActivityWork_[i];
    3516         double lowerScaled = rowLowerWork_[i];
    3517         double upperScaled = rowUpperWork_[i];
    3518         if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
    3519           if (valueScaled<lowerScaled-primalTolerance_||
    3520               valueScaled>upperScaled+primalTolerance_)
    3521             numberPrimalScaled++;
    3522           else
    3523             upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
    3524         }
    3525         rowActivity_[i] = (valueScaled*scaleR)/scaleFactor;
    3526         double value = rowActivity_[i];
    3527         if (value<rowLower_[i]-primalTolerance_)
    3528           numberPrimalUnscaled++;
    3529         else if (value>rowUpper_[i]+primalTolerance_)
    3530           numberPrimalUnscaled++;
    3531         double valueScaledDual = dual_[i]+rowObjectiveWork_[i];;
    3532         if (valueScaled>rowLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
    3533           numberDualScaled++;
    3534         if (valueScaled<rowUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
    3535           numberDualScaled++;
    3536         dual_[i] *= scaleFactor*scaleC;
    3537         double valueDual = dual_[i];
    3538         if (rowObjective_)
    3539           valueDual += rowObjective_[i];
    3540         if (value>rowLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
    3541           numberDualUnscaled++;
    3542         if (value<rowUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
    3543           numberDualUnscaled++;
     3559      if ((specialOptions_&131072)==0) {
     3560        for (i=0;i<numberColumns;i++) {
     3561          double scaleFactor = columnScale_[i];
     3562          double valueScaled = columnActivityWork_[i];
     3563          double lowerScaled = columnLowerWork_[i];
     3564          double upperScaled = columnUpperWork_[i];
     3565          if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
     3566            if (valueScaled<lowerScaled-primalTolerance_||
     3567                valueScaled>upperScaled+primalTolerance_)
     3568              numberPrimalScaled++;
     3569            else
     3570              upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
     3571          }
     3572          columnActivity_[i] = valueScaled*scaleFactor*scaleR;
     3573          double value = columnActivity_[i];
     3574          if (value<columnLower_[i]-primalTolerance_)
     3575            numberPrimalUnscaled++;
     3576          else if (value>columnUpper_[i]+primalTolerance_)
     3577            numberPrimalUnscaled++;
     3578          double valueScaledDual = reducedCostWork_[i];
     3579          if (valueScaled>columnLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
     3580            numberDualScaled++;
     3581          if (valueScaled<columnUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
     3582            numberDualScaled++;
     3583          reducedCost_[i] = (valueScaledDual*scaleC)/scaleFactor;
     3584          double valueDual = reducedCost_[i];
     3585          if (value>columnLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
     3586            numberDualUnscaled++;
     3587          if (value<columnUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
     3588            numberDualUnscaled++;
     3589        }
     3590        for (i=0;i<numberRows;i++) {
     3591          double scaleFactor = rowScale_[i];
     3592          double valueScaled = rowActivityWork_[i];
     3593          double lowerScaled = rowLowerWork_[i];
     3594          double upperScaled = rowUpperWork_[i];
     3595          if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
     3596            if (valueScaled<lowerScaled-primalTolerance_||
     3597                valueScaled>upperScaled+primalTolerance_)
     3598              numberPrimalScaled++;
     3599            else
     3600              upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
     3601          }
     3602          rowActivity_[i] = (valueScaled*scaleR)/scaleFactor;
     3603          double value = rowActivity_[i];
     3604          if (value<rowLower_[i]-primalTolerance_)
     3605            numberPrimalUnscaled++;
     3606          else if (value>rowUpper_[i]+primalTolerance_)
     3607            numberPrimalUnscaled++;
     3608          double valueScaledDual = dual_[i]+rowObjectiveWork_[i];;
     3609          if (valueScaled>rowLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
     3610            numberDualScaled++;
     3611          if (valueScaled<rowUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
     3612            numberDualScaled++;
     3613          dual_[i] *= scaleFactor*scaleC;
     3614          double valueDual = dual_[i];
     3615          if (rowObjective_)
     3616            valueDual += rowObjective_[i];
     3617          if (value>rowLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
     3618            numberDualUnscaled++;
     3619          if (value<rowUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
     3620            numberDualUnscaled++;
     3621        }
     3622      } else {
     3623        const double * inverseScale = columnScale_+numberColumns;
     3624        for (i=0;i<numberColumns;i++) {
     3625          double scaleFactor = columnScale_[i];
     3626          double valueScaled = columnActivityWork_[i];
     3627          double lowerScaled = columnLowerWork_[i];
     3628          double upperScaled = columnUpperWork_[i];
     3629          if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
     3630            if (valueScaled<lowerScaled-primalTolerance_||
     3631                valueScaled>upperScaled+primalTolerance_)
     3632              numberPrimalScaled++;
     3633            else
     3634              upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
     3635          }
     3636          columnActivity_[i] = valueScaled*scaleFactor*scaleR;
     3637          double value = columnActivity_[i];
     3638          if (value<columnLower_[i]-primalTolerance_)
     3639            numberPrimalUnscaled++;
     3640          else if (value>columnUpper_[i]+primalTolerance_)
     3641            numberPrimalUnscaled++;
     3642          double valueScaledDual = reducedCostWork_[i];
     3643          if (valueScaled>columnLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
     3644            numberDualScaled++;
     3645          if (valueScaled<columnUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
     3646            numberDualScaled++;
     3647          reducedCost_[i] = (valueScaledDual*scaleC)*inverseScale[i];
     3648          double valueDual = reducedCost_[i];
     3649          if (value>columnLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
     3650            numberDualUnscaled++;
     3651          if (value<columnUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
     3652            numberDualUnscaled++;
     3653        }
     3654        inverseScale = rowScale_+numberRows;
     3655        for (i=0;i<numberRows;i++) {
     3656          double scaleFactor = rowScale_[i];
     3657          double valueScaled = rowActivityWork_[i];
     3658          double lowerScaled = rowLowerWork_[i];
     3659          double upperScaled = rowUpperWork_[i];
     3660          if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
     3661            if (valueScaled<lowerScaled-primalTolerance_||
     3662                valueScaled>upperScaled+primalTolerance_)
     3663              numberPrimalScaled++;
     3664            else
     3665              upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
     3666          }
     3667          rowActivity_[i] = (valueScaled*scaleR)*inverseScale[i];
     3668          double value = rowActivity_[i];
     3669          if (value<rowLower_[i]-primalTolerance_)
     3670            numberPrimalUnscaled++;
     3671          else if (value>rowUpper_[i]+primalTolerance_)
     3672            numberPrimalUnscaled++;
     3673          double valueScaledDual = dual_[i]+rowObjectiveWork_[i];;
     3674          if (valueScaled>rowLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
     3675            numberDualScaled++;
     3676          if (valueScaled<rowUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
     3677            numberDualScaled++;
     3678          dual_[i] *= scaleFactor*scaleC;
     3679          double valueDual = dual_[i];
     3680          if (rowObjective_)
     3681            valueDual += rowObjective_[i];
     3682          if (value>rowLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
     3683            numberDualUnscaled++;
     3684          if (value<rowUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
     3685            numberDualUnscaled++;
     3686        }
    35443687      }
    35453688      if (!problemStatus_&&!secondaryStatus_) {
     
    35563699      }
    35573700      if (problemStatus_==2) {
    3558         for (i=0;i<numberColumns_;i++) {
     3701        for (i=0;i<numberColumns;i++) {
    35593702          ray_[i] *= columnScale_[i];
    35603703        }
    35613704      } else if (problemStatus_==1&&ray_) {
    3562         for (i=0;i<numberRows_;i++) {
     3705        for (i=0;i<numberRows;i++) {
    35633706          ray_[i] *= rowScale_[i];
    35643707        }
     
    35723715      double scaleC = 1.0/objectiveScale_;
    35733716      double scaleR = 1.0/rhsScale_;
    3574       for (i=0;i<numberColumns_;i++) {
     3717      for (i=0;i<numberColumns;i++) {
    35753718        double valueScaled = columnActivityWork_[i];
    35763719        double lowerScaled = columnLowerWork_[i];
     
    36013744          numberDualUnscaled++;
    36023745      }
    3603       for (i=0;i<numberRows_;i++) {
     3746      for (i=0;i<numberRows;i++) {
    36043747        double valueScaled = rowActivityWork_[i];
    36053748        double lowerScaled = rowLowerWork_[i];
     
    36463789    } else {
    36473790      if (columnActivityWork_) {
    3648         for (i=0;i<numberColumns_;i++) {
     3791        for (i=0;i<numberColumns;i++) {
    36493792          double value = columnActivityWork_[i];
    36503793          double lower = columnLowerWork_[i];
     
    36573800          reducedCost_[i] = reducedCostWork_[i];
    36583801        }
    3659         for (i=0;i<numberRows_;i++) {
     3802        for (i=0;i<numberRows;i++) {
    36603803          double value = rowActivityWork_[i];
    36613804          double lower = rowLowerWork_[i];
     
    36763819    if (optimizationDirection_!=1.0) {
    36773820      // and modify all dual signs
    3678       for (i=0;i<numberColumns_;i++)
     3821      for (i=0;i<numberColumns;i++)
    36793822        reducedCost_[i] *= optimizationDirection_;
    3680       for (i=0;i<numberRows_;i++)
     3823      for (i=0;i<numberRows;i++)
    36813824        dual_[i] *= optimizationDirection_;
    36823825    }
     
    37123855      if (auxiliaryModel_->rowScale_) {
    37133856        const double * scale = auxiliaryModel_->columnScale_;
    3714         const double * inverseScale = scale + numberColumns_;
    3715         for (i=0;i<numberColumns_;i++) {
     3857        const double * inverseScale = scale + numberColumns;
     3858        for (i=0;i<numberColumns;i++) {
    37163859          columnActivity_[i] = auxiliaryModel_->columnActivityWork_[i]*scale[i];
    37173860          reducedCost_[i] = auxiliaryModel_->reducedCostWork_[i]*inverseScale[i];
    37183861        }
    37193862        scale = auxiliaryModel_->rowScale_;
    3720         inverseScale = scale + numberRows_;
    3721         for (i=0;i<numberRows_;i++) {
     3863        inverseScale = scale + numberRows;
     3864        for (i=0;i<numberRows;i++) {
    37223865          rowActivity_[i] = auxiliaryModel_->rowActivityWork_[i]*inverseScale[i];
    37233866        }
    37243867      } else {
    3725         for (i=0;i<numberColumns_;i++) {
     3868        for (i=0;i<numberColumns;i++) {
    37263869          columnActivity_[i] = auxiliaryModel_->columnActivityWork_[i];
    37273870          reducedCost_[i] = auxiliaryModel_->reducedCostWork_[i];
    37283871        }
    3729         for (i=0;i<numberRows_;i++) {
     3872        for (i=0;i<numberRows;i++) {
    37303873          rowActivity_[i] = auxiliaryModel_->rowActivityWork_[i];
    37313874        }
     
    37333876      if (optimizationDirection_!=1.0) {
    37343877        // and modify reduced costs
    3735         for (i=0;i<numberColumns_;i++)
     3878        for (i=0;i<numberColumns;i++)
    37363879          reducedCost_[i] *= optimizationDirection_;
    37373880      }
     
    37403883      if (auxiliaryModel_->rowScale_) {
    37413884        const double * scale = auxiliaryModel_->columnScale_;
    3742         const double * inverseScale = scale + numberColumns_;
    3743         for (i=0;i<numberColumns_;i++) {
     3885        const double * inverseScale = scale + numberColumns;
     3886        for (i=0;i<numberColumns;i++) {
    37443887          double lower = auxiliaryModel_->columnLowerWork_[i];
    37453888          double upper = auxiliaryModel_->columnUpperWork_[i];
     
    37513894        }
    37523895        scale = auxiliaryModel_->rowScale_;
    3753         inverseScale = scale + numberRows_;
    3754         for (i=0;i<numberRows_;i++) {
     3896        inverseScale = scale + numberRows;
     3897        for (i=0;i<numberRows;i++) {
    37553898          double lower = auxiliaryModel_->rowLowerWork_[i];
    37563899          double upper = auxiliaryModel_->rowUpperWork_[i];
     
    37623905        }
    37633906      } else {
    3764         for (i=0;i<numberColumns_;i++) {
     3907        for (i=0;i<numberColumns;i++) {
    37653908          double lower = auxiliaryModel_->columnLowerWork_[i];
    37663909          double upper = auxiliaryModel_->columnUpperWork_[i];
     
    37713914          columnActivity_[i] = value;
    37723915        }
    3773         for (i=0;i<numberRows_;i++) {
     3916        for (i=0;i<numberRows;i++) {
    37743917          double lower = auxiliaryModel_->rowLowerWork_[i];
    37753918          double upper = auxiliaryModel_->rowUpperWork_[i];
     
    40874230                if (newBound > nowLower + 1.0e-12&&newBound>-large) {
    40884231                  // Tighten the lower bound
    4089                   columnLower_[iColumn] = newBound;
    40904232                  numberChanged++;
    40914233                  // check infeasible (relaxed)
    4092                   if (nowUpper - newBound <
    4093                       -100.0*tolerance) {
    4094                     numberInfeasible++;
     4234                  if (nowUpper < newBound) {
     4235                    if (nowUpper - newBound <
     4236                        -100.0*tolerance)
     4237                      numberInfeasible++;
     4238                    else
     4239                      newBound=nowUpper;
    40954240                  }
     4241                  columnLower_[iColumn] = newBound;
    40964242                  // adjust
    40974243                  double now;
     
    41344280                if (newBound < nowUpper - 1.0e-12&&newBound<large) {
    41354281                  // Tighten the upper bound
    4136                   columnUpper_[iColumn] = newBound;
    41374282                  numberChanged++;
    41384283                  // check infeasible (relaxed)
    4139                   if (newBound - nowLower <
    4140                       -100.0*tolerance) {
    4141                     numberInfeasible++;
     4284                  if (nowLower > newBound) {
     4285                    if (newBound - nowLower <
     4286                        -100.0*tolerance)
     4287                      numberInfeasible++;
     4288                    else
     4289                      newBound=nowLower;
    41424290                  }
     4291                  columnUpper_[iColumn] = newBound;
    41434292                  // adjust
    41444293                  double now;
     
    41834332                if (newBound < nowUpper - 1.0e-12&&newBound<large) {
    41844333                  // Tighten the upper bound
    4185                   columnUpper_[iColumn] = newBound;
    41864334                  numberChanged++;
    41874335                  // check infeasible (relaxed)
    4188                   if (newBound - nowLower <
    4189                       -100.0*tolerance) {
    4190                     numberInfeasible++;
     4336                  if (nowLower > newBound) {
     4337                    if (newBound - nowLower <
     4338                        -100.0*tolerance)
     4339                      numberInfeasible++;
     4340                    else
     4341                      newBound=nowLower;
    41914342                  }
     4343                  columnUpper_[iColumn] = newBound;
    41924344                  // adjust
    41934345                  double now;
     
    42304382                if (newBound > nowLower + 1.0e-12&&newBound>-large) {
    42314383                  // Tighten the lower bound
    4232                   columnLower_[iColumn] = newBound;
    42334384                  numberChanged++;
    42344385                  // check infeasible (relaxed)
    4235                   if (nowUpper - newBound <
    4236                       -100.0*tolerance) {
    4237                     numberInfeasible++;
     4386                  if (nowUpper < newBound) {
     4387                    if (nowUpper - newBound <
     4388                        -100.0*tolerance)
     4389                      numberInfeasible++;
     4390                    else
     4391                      newBound=nowUpper;
    42384392                  }
     4393                  columnLower_[iColumn] = newBound;
    42394394                  // adjust
    42404395                  double now;
     
    46824837{
    46834838  return ((ClpSimplexNonlinear *) this)->primalSLP(numberPasses,deltaTolerance);
     4839}
     4840/* Solves problem with nonlinear constraints using SLP - may be used as crash
     4841   for other algorithms when number of iterations small.
     4842   Also exits if all problematical variables are changing
     4843   less than deltaTolerance
     4844*/
     4845int
     4846ClpSimplex::nonlinearSLP(int numberConstraints, ClpConstraint ** constraints,
     4847                   int numberPasses,double deltaTolerance)
     4848{
     4849  return ((ClpSimplexNonlinear *) this)->primalSLP(numberConstraints,constraints,numberPasses,deltaTolerance);
    46844850}
    46854851// Solves non-linear using reduced gradient
     
    52015367    }
    52025368#endif
     5369    // integers
     5370    if (integerType_) {
     5371      int marker=1;
     5372      fwrite(&marker,sizeof(int),1,fp);
     5373      numberWritten = fwrite(integerType_,1,numberColumns_,fp);
     5374      if (numberWritten!=numberColumns_)
     5375        return 1;
     5376    } else {
     5377      int marker=0;
     5378      fwrite(&marker,sizeof(int),1,fp);
     5379    }
    52035380    // just standard type at present
    52045381    assert (matrix_->type()==1);
     
    52725449    }
    52735450    // get an empty factorization so we can set tolerances etc
    5274     factorization_ = new ClpFactorization();
     5451    getEmptyFactorization();
    52755452    // Say sparse
    52765453    factorization_->sparseThreshold(1);
     
    54055582    }
    54065583#endif
     5584    // integers
     5585    int ifInteger;
     5586    delete [] integerType_;
     5587    numberRead = fread(&ifInteger,sizeof(int),1,fp);
     5588    // But try and stay compatible with previous version
     5589    bool alreadyGotLength=false;
     5590    if (numberRead!=1)
     5591      return 1;
     5592    if (ifInteger==1) {
     5593      integerType_ = new char [numberColumns_];
     5594      numberRead = fread(integerType_,1,numberColumns_,fp);
     5595      if (numberRead!=numberColumns_)
     5596        return 1;
     5597    } else {
     5598      integerType_=NULL;
     5599      if (ifInteger) {
     5600        // probably old style save
     5601        alreadyGotLength=true;
     5602        length=ifInteger;
     5603      }
     5604    }
    54075605    // Pivot choices
    54085606    assert(scalars.dualPivotChoice>0&&(scalars.dualPivotChoice&63)<3);
     
    54395637    delete matrix_;
    54405638    // get arrays
    5441     numberRead = fread(&length,sizeof(int),1,fp);
    5442     if (numberRead!=1)
    5443       return 1;
     5639    if (!alreadyGotLength) {
     5640      numberRead = fread(&length,sizeof(int),1,fp);
     5641      if (numberRead!=1)
     5642        return 1;
     5643    }
    54445644    double * elements = new double[length];
    54455645    int * indices = new int[length];
     
    57115911      }
    57125912    }
     5913    if (solution_) {
     5914      // do that as well
     5915      if (!columnScale_) {
     5916        for (i=0;i<numberColumns_;i++) {
     5917          solution_[i] = columnActivity_[i];
     5918        }
     5919      } else {
     5920        for (i=0;i<numberColumns_;i++) {
     5921          solution_[i] = columnActivity_[i]*(rhsScale_/columnScale_[i]);
     5922        }
     5923      }
     5924    }
    57135925  }
    57145926}
     
    57835995  }
    57845996  int returnCode = ClpModel::loadProblem(modelObject);
     5997  const int * integerType = modelObject.integerTypeArray();
     5998  if (integerType) {
     5999    for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
     6000      if (integerType[iColumn])
     6001        setInteger(iColumn);
     6002    }
     6003  }
    57856004  createStatus();
    57866005  if (status) {
     
    58316050  createStatus();
    58326051  return status;
     6052}
     6053// Read file in LP format from file with name filename.
     6054int
     6055ClpSimplex::readLp(const char *filename, const double epsilon )
     6056{
     6057  FILE *fp = fopen(filename, "r");
     6058
     6059  if(!fp) {
     6060    printf("### ERROR: ClpSimplex::readLp():  Unable to open file %s for reading\n",
     6061           filename);
     6062    return(1);
     6063  }
     6064  CoinLpIO m;
     6065  m.readLp(fp, epsilon);
     6066  fclose(fp);
     6067 
     6068  // set problem name
     6069  setStrParam(ClpProbName, m.getProblemName());
     6070  // no errors
     6071  loadProblem(*m.getMatrixByRow(), m.getColLower(), m.getColUpper(),
     6072              m.getObjCoefficients(), m.getRowLower(), m.getRowUpper());
     6073 
     6074  if (m.integerColumns()) {
     6075    integerType_ = new char[numberColumns_];
     6076    CoinMemcpyN(m.integerColumns(),numberColumns_,integerType_);
     6077  } else {
     6078    integerType_ = NULL;
     6079  }
     6080  createStatus();
     6081  unsigned int maxLength=0;
     6082  int iRow;
     6083  rowNames_ = std::vector<std::string> ();
     6084  columnNames_ = std::vector<std::string> ();
     6085  rowNames_.reserve(numberRows_);
     6086  for (iRow=0;iRow<numberRows_;iRow++) {
     6087    const char * name = m.rowName(iRow);
     6088    if (name) {
     6089      maxLength = CoinMax(maxLength,(unsigned int) strlen(name));
     6090      rowNames_.push_back(name);
     6091    } else {
     6092      rowNames_.push_back("");
     6093    }
     6094  }
     6095 
     6096  int iColumn;
     6097  columnNames_.reserve(numberColumns_);
     6098  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     6099    const char * name = m.columnName(iColumn);
     6100    if (name) {
     6101      maxLength = CoinMax(maxLength,(unsigned int) strlen(name));
     6102      columnNames_.push_back(name);
     6103    } else {
     6104      columnNames_.push_back("");
     6105    }
     6106  }
     6107  lengthNames_=(int) maxLength;
     6108 
     6109  return 0;
    58336110}
    58346111#endif
     
    58936170            // set to nearest
    58946171            if (fabs(newValue-rowLower_[i])
    5895                 <fabs(newValue-rowLower_[i])) {
     6172                <fabs(newValue-rowUpper_[i])) {
    58966173              newValue=rowLower_[i];
    58976174              setRowStatus(i,atLowerBound);
     
    59646241            // set to nearest
    59656242            if (fabs(newValue-columnLower_[i])
    5966                 <fabs(newValue-columnLower_[i])) {
     6243                <fabs(newValue-columnUpper_[i])) {
    59676244              newValue=columnLower_[i];
    59686245              setColumnStatus(i,atLowerBound);
     
    65646841    if (sequenceOut_<0) {
    65656842      if (directionIn_<0) {
    6566         assert (fabs(solution_[sequenceIn_]-upperIn_)<1.0e-7);
     6843        assert (fabs(solution_[sequenceIn_]-upperIn_)<5.0e-7);
    65676844        solution_[sequenceIn_]=upperIn_;
    65686845      } else {
    6569         assert (fabs(solution_[sequenceIn_]-lowerIn_)<1.0e-7);
     6846        assert (fabs(solution_[sequenceIn_]-lowerIn_)<5.0e-7);
    65706847        solution_[sequenceIn_]=lowerIn_;
    65716848      }
    65726849    } else {
    65736850      if (directionOut_<0) {
    6574         assert (fabs(solution_[sequenceOut_]-upperOut_)<1.0e-7);
     6851        assert (fabs(solution_[sequenceOut_]-upperOut_)<5.0e-7);
    65756852        solution_[sequenceOut_]=upperOut_;
    65766853      } else {
    6577         assert (fabs(solution_[sequenceOut_]-lowerOut_)<1.0e-7);
     6854        assert (fabs(solution_[sequenceOut_]-lowerOut_)<5.0e-7);
    65786855        solution_[sequenceOut_]=lowerOut_;
    65796856      }
     
    67627039  // sanity check
    67637040  // bad if empty (trap here to avoid using bad matrix_)
     7041#if 0
     7042  // but also check bounds
     7043  {
     7044    int badProblem = 0;
     7045    int i;
     7046    for (i=0;i<numberColumns_;i++) {
     7047      if (columnLower_[i]>columnUpper_[i])
     7048        badProblem++;
     7049    }
     7050    for (i=0;i<numberRows_;i++) {
     7051      if (rowLower_[i]>rowUpper_[i])
     7052        badProblem++;
     7053    }
     7054    if (badProblem) {
     7055      numberDualInfeasibilities_=0;
     7056      sumDualInfeasibilities_=0.0;
     7057      numberPrimalInfeasibilities_=badProblem;
     7058      sumPrimalInfeasibilities_=badProblem;
     7059      secondaryStatus_=6; // so user can see something odd
     7060      problemStatus_=1;
     7061      bool printIt = (specialOptions_&32768)==0 ? true : false; // no message if from Osi
     7062      if (printIt)
     7063        handler_->message(CLP_INFEASIBLE,messages_)
     7064          <<CoinMessageEol;
     7065      return 2;
     7066    }
     7067  }
     7068#endif
    67647069  if (!matrix_||(!matrix_->getNumElements()&&objective_->type()<2)) {
    67657070    int infeasNumber[2];
     
    68737178    // number of times we have declared optimality
    68747179    numberTimesOptimal_=0;
     7180    if (disasterArea_)
     7181      disasterArea_->intoSimplex();
    68757182
    68767183    return 0;
     
    70217328  otherModel.largestPrimalError_ = largestPrimalError_;
    70227329  otherModel.largestDualError_ = largestDualError_;
    7023   otherModel.largestSolutionError_ = largestSolutionError_;
     7330  otherModel.alphaAccuracy_ = alphaAccuracy_;
    70247331  otherModel.alpha_ = alpha_;
    70257332  otherModel.theta_ = theta_;
     
    70447351  otherModel.numberPrimalInfeasibilities_ = numberPrimalInfeasibilities_;
    70457352  otherModel.numberTimesOptimal_ = numberTimesOptimal_;
     7353  otherModel.disasterArea_ = NULL;
    70467354  otherModel.sumOfRelaxedDualInfeasibilities_ = sumOfRelaxedDualInfeasibilities_;
    70477355  otherModel.sumOfRelaxedPrimalInfeasibilities_ = sumOfRelaxedPrimalInfeasibilities_;
     
    71107418    iterationNumber_[i]=-1;
    71117419  }
     7420#ifdef CLP_PROGRESS_WEIGHT
     7421  for (i=0;i<CLP_PROGRESS_WEIGHT;i++) {
     7422    objectiveWeight_[i] = COIN_DBL_MAX;
     7423    infeasibilityWeight_[i] = -1.0; // set to an impossible value
     7424    realInfeasibilityWeight_[i] = COIN_DBL_MAX;
     7425    numberInfeasibilitiesWeight_[i]=-1;
     7426    iterationNumberWeight_[i]=-1;
     7427  }
     7428  drop_ =0.0;
     7429  best_ =0.0;
     7430#endif
     7431  initialWeight_=0.0;
    71127432  for (i=0;i<CLP_CYCLE;i++) {
    71137433    //obj_[i]=COIN_DBL_MAX;
     
    71397459    iterationNumber_[i]=rhs.iterationNumber_[i];
    71407460  }
     7461#ifdef CLP_PROGRESS_WEIGHT
     7462  for (i=0;i<CLP_PROGRESS_WEIGHT;i++) {
     7463    objectiveWeight_[i] = rhs.objectiveWeight_[i];
     7464    infeasibilityWeight_[i] = rhs.infeasibilityWeight_[i];
     7465    realInfeasibilityWeight_[i] = rhs.realInfeasibilityWeight_[i];
     7466    numberInfeasibilitiesWeight_[i]=rhs.numberInfeasibilitiesWeight_[i];
     7467    iterationNumberWeight_[i]=rhs.iterationNumberWeight_[i];
     7468  }
     7469  drop_ = rhs.drop_;
     7470  best_ = rhs.best_;
     7471#endif
     7472  initialWeight_ = rhs.initialWeight_;
    71417473  for (i=0;i<CLP_CYCLE;i++) {
    71427474    //obj_[i]=rhs.obj_[i];
     
    71547486{
    71557487  model_ = model;
    7156   int i;
    7157   for (i=0;i<CLP_PROGRESS;i++) {
    7158     if (model_->algorithm()>=0)
    7159       objective_[i] = COIN_DBL_MAX;
    7160     else
    7161       objective_[i] = -COIN_DBL_MAX;
    7162     infeasibility_[i] = -1.0; // set to an impossible value
    7163     realInfeasibility_[i] = COIN_DBL_MAX;
    7164     numberInfeasibilities_[i]=-1;
    7165     iterationNumber_[i]=-1;
    7166   }
    7167   for (i=0;i<CLP_CYCLE;i++) {
    7168     //obj_[i]=COIN_DBL_MAX;
    7169     in_[i]=-1;
    7170     out_[i]=-1;
    7171     way_[i]=0;
    7172   }
    7173   numberTimes_ = 0;
    7174   numberBadTimes_ = 0;
    7175   oddState_=0;
     7488  reset();
     7489  initialWeight_=0.0;
    71767490}
    71777491// Assignment operator. This copies the data
     
    71887502      iterationNumber_[i]=rhs.iterationNumber_[i];
    71897503    }
     7504#ifdef CLP_PROGRESS_WEIGHT
     7505    for (i=0;i<CLP_PROGRESS_WEIGHT;i++) {
     7506      objectiveWeight_[i] = rhs.objectiveWeight_[i];
     7507      infeasibilityWeight_[i] = rhs.infeasibilityWeight_[i];
     7508      realInfeasibilityWeight_[i] = rhs.realInfeasibilityWeight_[i];
     7509      numberInfeasibilitiesWeight_[i]=rhs.numberInfeasibilitiesWeight_[i];
     7510      iterationNumberWeight_[i]=rhs.iterationNumberWeight_[i];
     7511    }
     7512    drop_ = rhs.drop_;
     7513    best_ = rhs.best_;
     7514#endif
     7515    initialWeight_ = rhs.initialWeight_;
    71907516    for (i=0;i<CLP_CYCLE;i++) {
    71917517      //obj_[i]=rhs.obj_[i];
     
    73627688  return -1;
    73637689}
     7690// Resets as much as possible
     7691void
     7692ClpSimplexProgress::reset()
     7693{
     7694  int i;
     7695  for (i=0;i<CLP_PROGRESS;i++) {
     7696    if (model_->algorithm()>=0)
     7697      objective_[i] = COIN_DBL_MAX;
     7698    else
     7699      objective_[i] = -COIN_DBL_MAX;
     7700    infeasibility_[i] = -1.0; // set to an impossible value
     7701    realInfeasibility_[i] = COIN_DBL_MAX;
     7702    numberInfeasibilities_[i]=-1;
     7703    iterationNumber_[i]=-1;
     7704  }
     7705#ifdef CLP_PROGRESS_WEIGHT
     7706  for (i=0;i<CLP_PROGRESS_WEIGHT;i++) {
     7707    objectiveWeight_[i] = COIN_DBL_MAX;
     7708    infeasibilityWeight_[i] = -1.0; // set to an impossible value
     7709    realInfeasibilityWeight_[i] = COIN_DBL_MAX;
     7710    numberInfeasibilitiesWeight_[i]=-1;
     7711    iterationNumberWeight_[i]=-1;
     7712  }
     7713  drop_ =0.