Changeset 1458 for stable


Ignore:
Timestamp:
Nov 5, 2009 7:34:07 AM (10 years ago)
Author:
forrest
Message:

Creating new stable branch 1.11 from trunk (rev 1457)

Location:
stable/1.11
Files:
112 edited
1 copied

Legend:

Unmodified
Added
Removed
  • stable/1.11/Clp/src/CbcOrClpParam.cpp

    r1404 r1458  
     1/* $Id$ */
    12// Copyright (C) 2002, International Business Machines
    23// Corporation and others.  All Rights Reserved.
     
    3334#else
    3435static char coin_prompt[]="Clp:";
     36#endif
     37#ifdef CLP_CILK
     38#ifndef CBC_THREAD
     39#define CBC_THREAD
     40#endif
    3541#endif
    3642#if WSSMP_BARRIER>=2
     
    6773    action_(INVALID),
    6874    currentKeyWord_(-1),
    69     display_(false),
     75    display_(0),
    7076    intValue_(-1),
    7177    doubleValue_(-1.0),
     
    7783CbcOrClpParam::CbcOrClpParam (std::string name, std::string help,
    7884           double lower, double upper, CbcOrClpParameterType type,
    79                     bool display)
     85                    int display)
    8086  : type_(type),
    8187    lowerIntValue_(0),
     
    99105CbcOrClpParam::CbcOrClpParam (std::string name, std::string help,
    100106           int lower, int upper, CbcOrClpParameterType type,
    101                     bool display)
     107                    int display)
    102108  : type_(type),
    103109    lowerDoubleValue_(0.0),
     
    123129                    std::string firstValue,
    124130                    CbcOrClpParameterType type,int whereUsed,
    125                     bool display)
     131                    int display)
    126132  : type_(type),
    127133    lowerDoubleValue_(0.0),
     
    147153CbcOrClpParam::CbcOrClpParam (std::string name, std::string help,
    148154                    CbcOrClpParameterType type,int whereUsed,
    149                     bool display)
     155                    int display)
    150156  : type_(type),
    151157    lowerDoubleValue_(0.0),
     
    232238{
    233239  std::string::size_type  shriekPos = name_.find('!');
    234   lengthName_ = (unsigned int)name_.length();
     240  lengthName_ = name_.length();
    235241  if ( shriekPos==std::string::npos ) {
    236242    //does not contain '!'
    237243    lengthMatch_= lengthName_;
    238244  } else {
    239     lengthMatch_=(unsigned int)shriekPos;
     245    lengthMatch_=shriekPos;
    240246    name_ = name_.substr(0,shriekPos)+name_.substr(shriekPos+1);
    241247    lengthName_--;
    242248  }
     249}
     250// Returns length of name for printing
     251int
     252CbcOrClpParam::lengthMatchName (  ) const
     253{
     254  if (lengthName_==lengthMatch_)
     255    return lengthName_;
     256  else
     257    return lengthName_+2;
    243258}
    244259// Insert string (only valid for keywords)
     
    285300CbcOrClpParam::parameterOption ( std::string check ) const
    286301{
    287   size_t numberItems = definedKeyWords_.size();
     302  int numberItems = definedKeyWords_.size();
    288303  if (!numberItems) {
    289304    return -1;
    290305  } else {
    291     size_t whichItem=0;
     306    int whichItem=0;
    292307    unsigned int it;
    293308    for (it=0;it<definedKeyWords_.size();it++) {
    294309      std::string thisOne = definedKeyWords_[it];
    295       std::string::size_type shriekPos = thisOne.find('!');
    296       std::string::size_type length1 = thisOne.length();
    297       std::string::size_type length2 = length1;
     310      std::string::size_type  shriekPos = thisOne.find('!');
     311      unsigned int length1 = thisOne.length();
     312      unsigned int length2 = length1;
    298313      if ( shriekPos!=std::string::npos ) {
    299314        //contains '!'
     
    304319      }
    305320      if (check.length()<=length1&&length2<=check.length()) {
    306         std::string::size_type i;
     321        unsigned int i;
    307322        for (i=0;i<check.length();i++) {
    308323          if (tolower(thisOne[i])!=tolower(check[i]))
     
    319334    }
    320335    if (whichItem<numberItems)
    321       return (int)whichItem;
     336      return whichItem;
    322337    else
    323338      return -1;
     
    364379void CoinReadPrintit(const char * input)
    365380{
    366   size_t length =strlen(input);
     381  int length =strlen(input);
    367382  char temp[101];
    368   size_t i;
     383  int i;
    369384  int n=0;
    370385  for (i=0;i<length;i++) {
     
    644659#ifdef COIN_HAS_CBC
    645660double
    646 CbcOrClpParam::doubleParameter (OsiSolverInterface * model) const
     661CbcOrClpParam::doubleParameter (OsiSolverInterface *
     662#ifndef NDEBUG
     663model
     664#endif
     665) const
    647666{
    648667  double value=0.0;
     
    858877      model.setNumberAnalyzeIterations(value);
    859878      break;
    860     case NUMBERMINI:
    861       oldValue=model.sizeMiniTree();
    862       model.setSizeMiniTree(value);
    863       break;
    864879    case CUTPASSINTREE:
    865880      oldValue=model.getMaximumCutPasses();
     
    912927  case NUMBERANALYZE:
    913928    value=model.numberAnalyzeIterations();
    914     break;
    915   case NUMBERMINI:
    916     value=model.sizeMiniTree();
    917929    break;
    918930  case CUTPASSINTREE:
     
    9991011static char * where=NULL;
    10001012extern int CbcOrClpRead_mode;
     1013int CbcOrClpEnvironmentIndex=-1;
     1014static int fillEnv()
     1015{
     1016  char * environ = getenv("CBC_CLP_ENVIRONMENT");
     1017  int length=0;
     1018  if (environ) {
     1019    length = strlen(environ);
     1020    if (CbcOrClpEnvironmentIndex<length) {
     1021      // find next non blank
     1022      char * whereEnv = environ+ CbcOrClpEnvironmentIndex;
     1023      // munch white space
     1024      while(*whereEnv==' '||*whereEnv=='\t'||*whereEnv<' ')
     1025        whereEnv++;
     1026      // copy
     1027      char * put = line;
     1028      while ( *whereEnv != '\0' ) {
     1029        if ( *whereEnv == ' '||*whereEnv == '\t' || *whereEnv < ' ' ) {
     1030          break;
     1031        }
     1032        *put=*whereEnv;
     1033        put++;
     1034        assert (put-line<1000);
     1035        whereEnv++;
     1036      }
     1037      CbcOrClpEnvironmentIndex=whereEnv-environ;
     1038      *put='\0';
     1039      length=strlen(line);
     1040    } else {
     1041      length=0;
     1042    }
     1043  }
     1044  if (!length)
     1045    CbcOrClpEnvironmentIndex=-1;
     1046  return length;
     1047}
    10011048extern FILE * CbcOrClpReadCommand;
    10021049// Simple read stuff
     
    10721119  while (field=="EOL") {
    10731120    if (CbcOrClpRead_mode>0) {
    1074       if (CbcOrClpRead_mode<argc&&argv[CbcOrClpRead_mode]) {
    1075         field = argv[CbcOrClpRead_mode++];
     1121      if ((CbcOrClpRead_mode<argc&&argv[CbcOrClpRead_mode])||
     1122          CbcOrClpEnvironmentIndex>=0) {
     1123        if(CbcOrClpEnvironmentIndex<0) {
     1124          field = argv[CbcOrClpRead_mode++];
     1125        } else {
     1126          if (fillEnv()) {
     1127            field=line;
     1128          } else {
     1129            // not there
     1130            continue;
     1131          }
     1132        }
    10761133        if (field=="-") {
    10771134          std::cout<<"Switching to line mode"<<std::endl;
     
    10821139            // now allow std::cout<<"skipping non-command "<<field<<std::endl;
    10831140            // field="EOL"; // skip
    1084           } else {
     1141          } else if (CbcOrClpEnvironmentIndex<0) {
    10851142            // special dispensation - taken as -import name
    10861143            CbcOrClpRead_mode--;
     
    11191176  if (afterEquals=="") {
    11201177    if (CbcOrClpRead_mode>0) {
    1121       if (CbcOrClpRead_mode<argc) {
    1122         if (argv[CbcOrClpRead_mode][0]!='-') {
    1123           field = argv[CbcOrClpRead_mode++];
    1124         } else if (!strcmp(argv[CbcOrClpRead_mode],"--")) {
    1125           field = argv[CbcOrClpRead_mode++];
    1126           // -- means import from stdin
    1127           field = "-";
    1128         }
     1178      if (CbcOrClpRead_mode<argc||CbcOrClpEnvironmentIndex>=0) {
     1179        if(CbcOrClpEnvironmentIndex<0) {
     1180          if (argv[CbcOrClpRead_mode][0]!='-') {
     1181            field = argv[CbcOrClpRead_mode++];
     1182          } else if (!strcmp(argv[CbcOrClpRead_mode],"--")) {
     1183            field = argv[CbcOrClpRead_mode++];
     1184            // -- means import from stdin
     1185            field = "-";
     1186          }
     1187        } else {
     1188          fillEnv();
     1189          field=line;
     1190        }
    11291191      }
    11301192    } else {
     
    11451207  if (afterEquals=="") {
    11461208    if (CbcOrClpRead_mode>0) {
    1147       if (CbcOrClpRead_mode<argc) {
    1148         // may be negative value so do not check for -
    1149         field = argv[CbcOrClpRead_mode++];
     1209      if (CbcOrClpRead_mode<argc||CbcOrClpEnvironmentIndex>=0) {
     1210        if(CbcOrClpEnvironmentIndex<0) {
     1211          // may be negative value so do not check for -
     1212          field = argv[CbcOrClpRead_mode++];
     1213        } else {
     1214          fillEnv();
     1215          field=line;
     1216        }
    11501217      }
    11511218    } else {
     
    11561223    afterEquals = "";
    11571224  }
    1158   size_t value=0;
     1225  int value=0;
    11591226  //std::cout<<field<<std::endl;
    11601227  if (field!="EOL") {
     
    11721239    *valid=2;
    11731240  }
    1174   return (int)value;
     1241  return value;
    11751242}
    11761243double
     
    11801247  if (afterEquals=="") {
    11811248    if (CbcOrClpRead_mode>0) {
    1182       if (CbcOrClpRead_mode<argc) {
    1183         // may be negative value so do not check for -
    1184         field = argv[CbcOrClpRead_mode++];
     1249      if (CbcOrClpRead_mode<argc||CbcOrClpEnvironmentIndex>=0) {
     1250        if(CbcOrClpEnvironmentIndex<0) {
     1251          // may be negative value so do not check for -
     1252          field = argv[CbcOrClpRead_mode++];
     1253        } else {
     1254          fillEnv();
     1255          field=line;
     1256        }
    11851257      }
    11861258    } else {
     
    12181290  numberParameters=0;
    12191291  parameters[numberParameters++]=
    1220     CbcOrClpParam("?","For help",GENERALQUERY,7,false);
    1221   parameters[numberParameters++]=
    1222     CbcOrClpParam("???","For help",FULLGENERALQUERY,7,false);
     1292    CbcOrClpParam("?","For help",GENERALQUERY,7,0);
     1293  parameters[numberParameters++]=
     1294    CbcOrClpParam("???","For help",FULLGENERALQUERY,7,0);
    12231295  parameters[numberParameters++]=
    12241296    CbcOrClpParam("-","From stdin",
    1225                   STDIN,3,false);
     1297                  STDIN,3,0);
     1298  parameters[numberParameters++]=
     1299    CbcOrClpParam("allC!ommands","Whether to print less used commands",
     1300                  "no",ALLCOMMANDS);
     1301  parameters[numberParameters-1].append("more");
     1302  parameters[numberParameters-1].append("all");
     1303  parameters[numberParameters-1].setLonghelp
     1304    (
     1305     "For the sake of your sanity, only the more useful and simple commands \
     1306are printed out on ?."
     1307     );
    12261308#ifdef COIN_HAS_CBC
    12271309    parameters[numberParameters++]=
     
    12451327basis anyway."
    12461328     );
     1329#endif
     1330#ifdef COIN_HAS_CBC
     1331  parameters[numberParameters++]=
     1332    CbcOrClpParam("artif!icialCost","Costs >= this treated as artificials in feasibility pump",
     1333                  0.0,COIN_DBL_MAX,ARTIFICIALCOST,1);
     1334  parameters[numberParameters-1].setDoubleValue(0.0);
     1335    parameters[numberParameters-1].setLonghelp
     1336    (
     1337     "0.0 off - otherwise variables with costs >= this are treated as artificials and fixed to lower bound in feasibility pump"
     1338     );
     1339#endif
     1340#ifdef COIN_HAS_CLP
    12471341  parameters[numberParameters++]=
    12481342    CbcOrClpParam("auto!Scale","Whether to scale objective, rhs and bounds of problem if they look odd",
    1249                   "off",AUTOSCALE,7,false);
     1343                  "off",AUTOSCALE,7,0);
    12501344  parameters[numberParameters-1].append("on");
    12511345  parameters[numberParameters-1].setLonghelp
     
    12531347     "If you think you may get odd objective values or large equality rows etc then\
    12541348 it may be worth setting this true.  It is still experimental and you may prefer\
    1255  to use objective!Scale and rhs!Scale"
     1349 to use objective!Scale and rhs!Scale."
    12561350     );
    12571351  parameters[numberParameters++]=
     
    12611355    (
    12621356     "This command solves the current model using the  primal dual predictor \
    1263 corrector algorithm.  You will probably want to link in and choose a better \
    1264 ordering and factorization than the default ones provided.  It will also solve models \
     1357corrector algorithm.  You may want to link in an alternative \
     1358ordering and factorization.  It will also solve models \
    12651359with quadratic objectives."
    12661360     
     
    12741368 directory given by 'directory'.  A name of '$' will use the previous value for the name.  This\
    12751369 is initialized to '', i.e. it must be set.  If you have libz then it can read compressed\
    1276  files 'xxxxxxxx.gz'.."
     1370 files 'xxxxxxxx.gz' or xxxxxxxx.bz2."
    12771371     );
    12781372  parameters[numberParameters++]=
     
    12871381  parameters[numberParameters++]=
    12881382    CbcOrClpParam("biasLU","Whether factorization biased towards U",
    1289                   "UU",BIASLU,2,false);
     1383                  "UU",BIASLU,2,0);
    12901384  parameters[numberParameters-1].append("UX");
    12911385  parameters[numberParameters-1].append("LX");
     
    13121406were obtained fairly early in the search so the important point is to select the best variable to branch on.  \
    13131407See whether strong branching did a good job - or did it just take a lot of iterations.  Adjust the strongBranching \
    1314 and trustPseudoCosts parameters."
    1315      );
    1316 #endif
    1317   parameters[numberParameters++]=
    1318     CbcOrClpParam("bscale","Whether to scale in barrier",
    1319                   "off",BARRIERSCALE,7,false);
    1320   parameters[numberParameters-1].append("on");
     1408and trustPseudoCosts parameters.  If cuts did a good job, then you may wish to \
     1409have more rounds of cuts - see passC!uts and passT!ree."
     1410     );
     1411#endif
     1412  parameters[numberParameters++]=
     1413    CbcOrClpParam("bscale","Whether to scale in barrier (and ordering speed)",
     1414                  "off",BARRIERSCALE,7,0);
     1415  parameters[numberParameters-1].append("on");
     1416  parameters[numberParameters-1].append("off1");
     1417  parameters[numberParameters-1].append("on1");
     1418  parameters[numberParameters-1].append("off2");
     1419  parameters[numberParameters-1].append("on2");
    13211420  parameters[numberParameters++]=
    13221421    CbcOrClpParam("chol!esky","Which cholesky algorithm",
     
    13271426  parameters[numberParameters-1].append("fudge!Long");
    13281427  parameters[numberParameters-1].append("wssmp");
    1329 #define REAL_BARRIER
    13301428#else
    13311429  parameters[numberParameters-1].append("fudge!Long_dummy");
     
    13341432#ifdef UFL_BARRIER
    13351433  parameters[numberParameters-1].append("Uni!versityOfFlorida");
    1336 #define REAL_BARRIER
    13371434#else
    13381435  parameters[numberParameters-1].append("Uni!versityOfFlorida_dummy");   
     
    13401437#ifdef TAUCS_BARRIER
    13411438  parameters[numberParameters-1].append("Taucs");
    1342 #define REAL_BARRIER
    13431439#else
    13441440  parameters[numberParameters-1].append("Taucs_dummy");
     
    13461442#ifdef MUMPS_BARRIER
    13471443  parameters[numberParameters-1].append("Mumps");
    1348 #define REAL_BARRIER
    13491444#else
    13501445  parameters[numberParameters-1].append("Mumps_dummy");   
     
    13531448    (
    13541449     "For a barrier code to be effective it needs a good Cholesky ordering and factorization.  \
    1355 You will need to link in one from another source.  See Makefile.locations for some \
     1450The native ordering and factorization is not state of the art, although acceptable.  \
     1451You may want to link in one from another source.  See Makefile.locations for some \
    13561452possibilities."
    13571453     );
     
    13661462  parameters[numberParameters-1].append("forceOn");
    13671463  parameters[numberParameters-1].append("onglobal");
    1368   parameters[numberParameters-1].append("rootglobal");
    13691464  parameters[numberParameters-1].setLonghelp
    13701465    (
     
    13821477     "This switches on a heuristic which does branch and cut on the problem given by just \
    13831478using variables which have appeared in one or more solutions. \
     1479It obviously only tries after two or more solutions. \
     1480See Rounding for meaning of on,both,before"
     1481     );
     1482  parameters[numberParameters++]=
     1483    CbcOrClpParam("combine2!Solutions","Whether to use crossover solution heuristic",
     1484                  "off",CROSSOVER2);
     1485  parameters[numberParameters-1].append("on");
     1486  parameters[numberParameters-1].append("both");
     1487  parameters[numberParameters-1].append("before");
     1488  parameters[numberParameters-1].setLonghelp
     1489    (
     1490     "This switches on a heuristic which does branch and cut on the problem given by \
     1491fixing variables which have same value in two or more solutions. \
    13841492It obviously only tries after two or more solutions. \
    13851493See Rounding for meaning of on,both,before"
     
    14051513  parameters[numberParameters-1].setLonghelp
    14061514    (
    1407      " If the user has Cplex, but wants to use some of Cbc'sheuristics \
     1515     " If the user has Cplex, but wants to use some of Cbc's heuristics \
    14081516then you can!  If this is on, then Cbc will get to the root node and then \
    14091517hand over to Cplex.  If heuristics find a solution this can be significantly \
     
    14161524  parameters[numberParameters++]=
    14171525    CbcOrClpParam("cpp!Generate","Generates C++ code",
    1418                   -1,50000,CPP);
     1526                  -1,50000,CPP,1);
    14191527  parameters[numberParameters-1].setLonghelp
    14201528    (
     
    142315310 gives simplest driver, 1 generates saves and restores, 2 \
    14241532generates saves and restores even for variables at default value. \
    1425 4 bit in cbc generates size dependent code rather than computed values."
     15334 bit in cbc generates size dependent code rather than computed values.  \
     1534This is now deprecated as you can call stand-alone solver - see \
     1535Cbc/examples/driver4.cpp."
    14261536     );
    14271537#ifdef COIN_HAS_CLP
     
    14381548     "If crash is set on and there is an all slack basis then Clp will flip or put structural\
    14391549 variables into basis with the aim of getting dual feasible.  On the whole dual seems to be\
    1440  better without it and there alernative types of 'crash' for primal e.g. 'idiot' or 'sprint'. \
     1550 better without it and there are alternative types of 'crash' for primal e.g. 'idiot' or 'sprint'. \
    14411551I have also added a variant due to Solow and Halim which is as on but just flip.");
    14421552  parameters[numberParameters++]=
     
    14581568  parameters[numberParameters++]=
    14591569    CbcOrClpParam("csv!Statistics","Create one line of statistics",
    1460                   CSVSTATISTICS,2,false);
     1570                  CSVSTATISTICS,2,1);
    14611571  parameters[numberParameters-1].setLonghelp
    14621572    (
     
    14741584every so often.  The original method was every so many nodes but it is more logical \
    14751585to do it whenever depth in tree is a multiple of K.  This option does that and defaults \
    1476 to -1 (off)."
     1586to -1 (off -> code decides)."
     1587     );
     1588  parameters[numberParameters-1].setIntValue(-1);
     1589  parameters[numberParameters++]=
     1590    CbcOrClpParam("cutL!ength","Length of a cut",
     1591                  -1,COIN_INT_MAX,CUTLENGTH);
     1592  parameters[numberParameters-1].setLonghelp
     1593    (
     1594     "At present this only applies to Gomory cuts. -1 (default) leaves as is. \
     1595Any value >0 says that all cuts <= this length can be generated both at \
     1596root node and in tree. 0 says to use some dynamic lengths.  If value >=10,000,000 \
     1597then the length in tree is value%10000000 - so 10000100 means unlimited length \
     1598at root and 100 in tree."
    14771599     );
    14781600  parameters[numberParameters-1].setIntValue(-1);
     
    15021624  parameters[numberParameters++]=
    15031625    CbcOrClpParam("debug!In","read valid solution from file",
    1504                   DEBUG,7,false);
     1626                  DEBUG,7,1);
    15051627  parameters[numberParameters-1].setLonghelp
    15061628    (
     
    15181640  parameters[numberParameters++]=
    15191641    CbcOrClpParam("dense!Threshold","Whether to use dense factorization",
    1520                   -1,10000,DENSE,false);
     1642                  -1,10000,DENSE,1);
    15211643  parameters[numberParameters-1].setLonghelp
    15221644    (
     
    15281650#ifdef COIN_HAS_CBC
    15291651  parameters[numberParameters++]=
    1530     CbcOrClpParam("dextra1","Extra double parameter 1",
    1531                   -COIN_DBL_MAX,COIN_DBL_MAX,DEXTRA1,false);
     1652    CbcOrClpParam("depth!MiniBab","Depth at which to try mini BAB",
     1653                  -COIN_INT_MAX,COIN_INT_MAX,DEPTHMINIBAB);
     1654  parameters[numberParameters-1].setIntValue(-1);
     1655  parameters[numberParameters-1].setLonghelp
     1656    (
     1657     "Rather a complicated parameter but can be useful. -1 means off for large problems but on as if -12 for problems where rows+columns<500, -2 \
     1658means use Cplex if it is linked in.  Otherwise if negative then go into depth first complete search fast branch and bound when depth>= -value-2 (so -3 will use this at depth>=1).  This mode is only switched on after 500 nodes.  If you really want to switch it off for small problems then set this to -999.  If >=0 the value doesn't matter very much.  The code will do approximately 100 nodes of fast branch and bound every now and then at depth>=5.  The actual logic is too twisted to describe here."
     1659     );
     1660  parameters[numberParameters++]=
     1661    CbcOrClpParam("dextra3","Extra double parameter 3",
     1662                  -COIN_DBL_MAX,COIN_DBL_MAX,DEXTRA3,0);
    15321663  parameters[numberParameters-1].setDoubleValue(0.0);
    15331664  parameters[numberParameters++]=
    1534     CbcOrClpParam("dextra2","Extra double parameter 2",
    1535                   -COIN_DBL_MAX,COIN_DBL_MAX,DEXTRA2,false);
     1665    CbcOrClpParam("dextra4","Extra double parameter 4",
     1666                  -COIN_DBL_MAX,COIN_DBL_MAX,DEXTRA4,0);
    15361667  parameters[numberParameters-1].setDoubleValue(0.0);
    15371668  parameters[numberParameters++]=
    1538     CbcOrClpParam("dextra3","Extra double parameter 3",
    1539                   -COIN_DBL_MAX,COIN_DBL_MAX,DEXTRA3,false);
    1540   parameters[numberParameters-1].setDoubleValue(0.0);
    1541   parameters[numberParameters++]=
    1542     CbcOrClpParam("dextra4","Extra double parameter 4",
    1543                   -COIN_DBL_MAX,COIN_DBL_MAX,DEXTRA4,false);
    1544   parameters[numberParameters-1].setDoubleValue(0.0);
    1545   parameters[numberParameters++]=
    15461669    CbcOrClpParam("dextra5","Extra double parameter 5",
    1547                   -COIN_DBL_MAX,COIN_DBL_MAX,DEXTRA5,false);
     1670                  -COIN_DBL_MAX,COIN_DBL_MAX,DEXTRA5,0);
    15481671  parameters[numberParameters-1].setDoubleValue(0.0);
    15491672  parameters[numberParameters++]=
     
    15801703  parameters[numberParameters++]=
    15811704    CbcOrClpParam("dirSample","Set directory where the COIN-OR sample problems are.",
    1582                   DIRSAMPLE);
     1705                  DIRSAMPLE,7,1);
    15831706  parameters[numberParameters-1].setLonghelp
    15841707    (
     
    15901713  parameters[numberParameters++]=
    15911714    CbcOrClpParam("dirNetlib","Set directory where the netlib problems are.",
    1592                   DIRNETLIB);
     1715                  DIRNETLIB,7,1);
    15931716  parameters[numberParameters-1].setLonghelp
    15941717    (
     
    16021725  parameters[numberParameters++]=
    16031726    CbcOrClpParam("dirMiplib","Set directory where the miplib 2003 problems are.",
    1604                   DIRMIPLIB);
     1727                  DIRMIPLIB,7,1);
    16051728  parameters[numberParameters-1].setLonghelp
    16061729    (
     
    16151738  parameters[numberParameters++]=
    16161739    CbcOrClpParam("diveO!pt","Diving options",
    1617                   -1,200,DIVEOPT,false);
    1618   parameters[numberParameters-1].setLonghelp
    1619     (
    1620      "If >2 && <7 then modify diving options"
     1740                  -1,200000,DIVEOPT,1);
     1741  parameters[numberParameters-1].setLonghelp
     1742    (
     1743     "If >2 && <8 then modify diving options - \
     1744         \t3 only at root and if no solution,  \
     1745         \t4 only at root and if this heuristic has not got solution, \
     1746         \t5 only at depth <4, \
     1747         \t6 decay, \
     1748         \t7 run up to 2 times if solution found 4 otherwise."
    16211749     );
    16221750  parameters[numberParameters-1].setIntValue(-1);
     
    16761804    (
    16771805     "Normally heuristics are done in branch and bound.  It may be useful to do them outside. \
    1678 Doing this may also set cutoff."
     1806Only those heuristics with 'both' or 'before' set will run.  \
     1807Doing this may also set cutoff, which can help with preprocessing."
    16791808     );
    16801809#endif
     
    16981827  parameters[numberParameters++]=
    16991828    CbcOrClpParam("dualize","Solves dual reformulation",
    1700                   0,3,DUALIZE,false);
     1829                  0,3,DUALIZE,1);
    17011830  parameters[numberParameters-1].setLonghelp
    17021831    (
     
    17051834  parameters[numberParameters++]=
    17061835    CbcOrClpParam("dualP!ivot","Dual pivot choice algorithm",
    1707                   "auto!matic",DUALPIVOT);
     1836                  "auto!matic",DUALPIVOT,7,1);
    17081837  parameters[numberParameters-1].append("dant!zig");
    17091838  parameters[numberParameters-1].append("partial");
     
    17571886     );
    17581887  parameters[numberParameters++]=
     1888    CbcOrClpParam("environ!ment","Read commands from environment",
     1889                  ENVIRONMENT,7,0);
     1890  parameters[numberParameters-1].setLonghelp
     1891    (
     1892     "This starts reading from environment variable CBC_CLP_ENVIRONMENT."
     1893     );
     1894  parameters[numberParameters++]=
    17591895    CbcOrClpParam("error!sAllowed","Whether to allow import errors",
    17601896                  "off",ERRORSALLOWED,3);
     
    17771913  parameters[numberParameters++]=
    17781914    CbcOrClpParam("exp!eriment","Whether to use testing features",
    1779                   -1,200,EXPERIMENT,false);
     1915                  -1,200,EXPERIMENT,0);
    17801916  parameters[numberParameters-1].setLonghelp
    17811917    (
     
    17921928     "This will write an MPS format file to the given file name.  It will use the default\
    17931929 directory given by 'directory'.  A name of '$' will use the previous value for the name.  This\
    1794  is initialized to 'default.mps'."
     1930 is initialized to 'default.mps'.  \
     1931It can be useful to get rid of the original names and go over to using Rnnnnnnn and Cnnnnnnn.  This can be done by setting 'keepnames' off before importing mps file."
    17951932     );
    17961933#ifdef COIN_HAS_CBC
    17971934  parameters[numberParameters++]=
    17981935    CbcOrClpParam("extra1","Extra integer parameter 1",
    1799                   -COIN_INT_MAX,COIN_INT_MAX,EXTRA1,false);
     1936                  -COIN_INT_MAX,COIN_INT_MAX,EXTRA1,0);
    18001937  parameters[numberParameters-1].setIntValue(-1);
    18011938  parameters[numberParameters++]=
    18021939    CbcOrClpParam("extra2","Extra integer parameter 2",
    1803                   -100,COIN_INT_MAX,EXTRA2,false);
     1940                  -100,COIN_INT_MAX,EXTRA2,0);
    18041941  parameters[numberParameters-1].setIntValue(-1);
    18051942  parameters[numberParameters++]=
    18061943    CbcOrClpParam("extra3","Extra integer parameter 3",
    1807                   -1,COIN_INT_MAX,EXTRA3,false);
     1944                  -1,COIN_INT_MAX,EXTRA3,0);
    18081945  parameters[numberParameters-1].setIntValue(-1);
    18091946  parameters[numberParameters++]=
    18101947    CbcOrClpParam("extra4","Extra integer parameter 4",
    1811                   -COIN_INT_MAX,COIN_INT_MAX,EXTRA4,false);
     1948                  -1,COIN_INT_MAX,EXTRA4,0);
    18121949  parameters[numberParameters-1].setIntValue(-1);
     1950  parameters[numberParameters-1].setLonghelp
     1951    (
     1952     "This switches on yet more special options!! \
     1953The bottom digit is a strategy when to used shadow price stuff e.g. 3 \
     1954means use until a solution is found.  The next two digits say what sort \
     1955of dual information to use.  After that it goes back to powers of 2 so -\n\
     1956\t1000 - switches on experimental hotstart\n\
     1957\t2,4,6000 - switches on experimental methods of stopping cuts\n\
     1958\t8000 - increase minimum drop gradually\n\
     1959\t16000 - switches on alternate gomory criterion"
     1960     );
    18131961#endif
    18141962#ifdef COIN_HAS_CLP
    18151963  parameters[numberParameters++]=
     1964    CbcOrClpParam("fact!orization","Which factorization to use",
     1965                  "normal",FACTORIZATION);
     1966  parameters[numberParameters-1].append("dense");
     1967  parameters[numberParameters-1].append("simple");
     1968  parameters[numberParameters-1].append("osl");
     1969  parameters[numberParameters-1].setLonghelp
     1970    (
     1971     "The default is to use the normal CoinFactorization, but \
     1972other choices are a dense one, osl's or one designed for small problems."
     1973     );
     1974  parameters[numberParameters++]=
    18161975    CbcOrClpParam("fakeB!ound","All bounds <= this value - DEBUG",
    1817                   1.0,1.0e15,FAKEBOUND,false);
     1976                  1.0,1.0e15,FAKEBOUND,0);
    18181977#ifdef COIN_HAS_CBC
    18191978  parameters[numberParameters++]=
     
    18331992    CbcOrClpParam("fix!OnDj","Try heuristic based on fixing variables with \
    18341993reduced costs greater than this",
    1835                   -1.0e20,1.0e20,DJFIX,false);
     1994                  -1.0e20,1.0e20,DJFIX,1);
    18361995  parameters[numberParameters-1].setLonghelp
    18371996    (
     
    18462005    parameters[numberParameters-1].append("ifmove");
    18472006    parameters[numberParameters-1].append("forceOn");
    1848   parameters[numberParameters-1].append("onglobal");
    1849   parameters[numberParameters-1].append("rootglobal");
     2007    parameters[numberParameters-1].append("onglobal");
    18502008    parameters[numberParameters-1].setLonghelp
    18512009    (
     
    18612019     "-1 off.  If 1 then tries to branch to solution given by AMPL or priorities file. \
    18622020If 0 then just tries to set as best solution \
    1863 If 1 then also does that many nodes on fixed problem."
     2021If >1 then also does that many nodes on fixed problem."
     2022     );
     2023  parameters[numberParameters++]=
     2024    CbcOrClpParam("fraction!forBAB","Fraction in feasibility pump",
     2025                  1.0e-5,1.1,SMALLBAB,1);
     2026  parameters[numberParameters-1].setDoubleValue(0.5);
     2027  parameters[numberParameters-1].setLonghelp
     2028    (
     2029     "After a pass in feasibility pump, variables which have not moved \
     2030about are fixed and if the preprocessed model is small enough a few nodes \
     2031of branch and bound are done on reduced problem.  Small problem has to be less than this fraction of original."
    18642032     );
    18652033#endif
    18662034  parameters[numberParameters++]=
    18672035    CbcOrClpParam("gamma!(Delta)","Whether to regularize barrier",
    1868                   "off",GAMMA,7,false);
     2036                  "off",GAMMA,7,1);
    18692037  parameters[numberParameters-1].append("on");
    18702038  parameters[numberParameters-1].append("gamma");
     
    18832051  parameters[numberParameters-1].append("forceOn");
    18842052  parameters[numberParameters-1].append("onglobal");
    1885   parameters[numberParameters-1].append("rootglobal");
     2053  parameters[numberParameters-1].append("forceandglobal");
    18862054  parameters[numberParameters-1].append("forceLongOn");
    18872055  parameters[numberParameters-1].append("long");
     
    19292097  parameters[numberParameters++]=
    19302098    CbcOrClpParam("hOp!tions","Heuristic options",
    1931                   -9999999,9999999,HOPTIONS);
     2099                  -9999999,9999999,HOPTIONS,1);
    19322100  parameters[numberParameters-1].setLonghelp
    19332101    (
     
    19412109  parameters[numberParameters++]=
    19422110    CbcOrClpParam("hot!StartMaxIts","Maximum iterations on hot start",
    1943                   0,COIN_INT_MAX,MAXHOTITS,false);
     2111                  0,COIN_INT_MAX,MAXHOTITS);
    19442112#endif
    19452113#ifdef COIN_HAS_CLP
    19462114  parameters[numberParameters++]=
    19472115    CbcOrClpParam("idiot!Crash","Whether to try idiot crash",
    1948                   -1,999999,IDIOT);
     2116                  -1,99999999,IDIOT);
    19492117  parameters[numberParameters-1].setLonghelp
    19502118    (
     
    19632131 directory given by 'directory'.  A name of '$' will use the previous value for the name.  This\
    19642132 is initialized to '', i.e. it must be set.  If you have libgz then it can read compressed\
    1965  files 'xxxxxxxx.gz'.."
     2133 files 'xxxxxxxx.gz' or 'xxxxxxxx.bz2'.  \
     2134If 'keepnames' is off, then names are dropped -> Rnnnnnnn and Cnnnnnnn."
    19662135     );
    19672136#ifdef COIN_HAS_CBC
     
    19802149    CbcOrClpParam("inf!easibilityWeight","Each integer infeasibility is expected \
    19812150to cost this much",
    1982                   0.0,1.0e20,INFEASIBILITYWEIGHT);
     2151                  0.0,1.0e20,INFEASIBILITYWEIGHT,1);
    19832152  parameters[numberParameters-1].setLonghelp
    19842153    (
     
    20142183  parameters[numberParameters++]=
    20152184    CbcOrClpParam("KKT","Whether to use KKT factorization",
    2016                   "off",KKT,7,false);
     2185                  "off",KKT,7,1);
    20172186  parameters[numberParameters-1].append("on");
    20182187#endif
     
    20262195  parameters[numberParameters-1].append("forceOn");
    20272196  parameters[numberParameters-1].append("onglobal");
    2028   parameters[numberParameters-1].append("rootglobal");
     2197  parameters[numberParameters-1].append("forceandglobal");
    20292198  parameters[numberParameters-1].setLonghelp
    20302199    (
     
    20982267     "This can be used for testing purposes.  The corresponding library call\n\
    20992268      \tsetMaximumIterations(value)\n can be useful.  If the code stops on\
    2100  seconds or by an interrupt this will be treated as stopping on maximum iterations"
     2269 seconds or by an interrupt this will be treated as stopping on maximum iterations.  This is ignored in branchAndCut - use maxN!odes."
    21012270     );
    21022271#endif
     
    21152284  parameters[numberParameters-1].setLonghelp
    21162285    (
    2117      "You may want to stop after (say) two solutions or an hour."
     2286     "You may want to stop after (say) two solutions or an hour.  \
     2287This is checked every node in tree, so it is possible to get more solutions from heuristics."
    21182288     );
    21192289#endif
     
    21302300  parameters[numberParameters++]=
    21312301    CbcOrClpParam("mipO!ptions","Dubious options for mip",
    2132                   0,COIN_INT_MAX,MIPOPTIONS,false);
     2302                  0,COIN_INT_MAX,MIPOPTIONS,0);
    21332303  parameters[numberParameters++]=
    21342304    CbcOrClpParam("more!MipOptions","More dubious options for mip",
    2135                   -1,COIN_INT_MAX,MOREMIPOPTIONS,false);
     2305                  -1,COIN_INT_MAX,MOREMIPOPTIONS,0);
    21362306  parameters[numberParameters++]=
    21372307    CbcOrClpParam("mixed!IntegerRoundingCuts","Whether to use Mixed Integer Rounding cuts",
     
    21422312  parameters[numberParameters-1].append("forceOn");
    21432313  parameters[numberParameters-1].append("onglobal");
    2144   parameters[numberParameters-1].append("rootglobal");
    21452314  parameters[numberParameters-1].setLonghelp
    21462315    (
     
    21582327but this program turns this off to make it look more friendly.  It can be useful\
    21592328 to turn them back on if you want to be able to 'grep' for particular messages or if\
    2160  you intend to override the behavior of a particular message."
     2329 you intend to override the behavior of a particular message.  This only affects Clp not Cbc."
    21612330     );
    21622331#ifdef COIN_HAS_CBC
    21632332  parameters[numberParameters++]=
    2164     CbcOrClpParam("miniT!ree","Size of fast mini tree",
    2165                   0,COIN_INT_MAX,NUMBERMINI,false);
    2166   parameters[numberParameters-1].setLonghelp
    2167     (
    2168      "The idea is that I can do a small tree fast. \
    2169 This is a first try and will hopefully become more sophisticated."
    2170      );
     2333    CbcOrClpParam("moreT!une","Yet more dubious ideas for feasibility pump",
     2334                  0,100000000,FPUMPTUNE2,0);
     2335  parameters[numberParameters-1].setLonghelp
     2336    (
     2337     "Yet more ideas for Feasibility Pump \n\
     2338\t/100000 == 1 use box constraints and original obj in cleanup\n\
     2339\t/1000 == 1 Pump will run twice if no solution found\n\
     2340\t/1000 == 2 Pump will only run after root cuts if no solution found\n\
     2341\t/1000 >10 as above but even if solution found\n\
     2342\t/100 == 1,3.. exact 1.0 for objective values\n\
     2343\t/100 == 2,3.. allow more iterations per pass\n\
     2344\t n fix if value of variable same for last n iterations."
     2345     );
     2346  parameters[numberParameters-1].setIntValue(0);
    21712347  parameters[numberParameters++]=
    21722348    CbcOrClpParam("miplib","Do some of miplib test set",
    2173                   MIPLIB,3);
     2349                  MIPLIB,3,1);
    21742350  parameters[numberParameters++]=
    21752351      CbcOrClpParam("naive!Heuristics","Whether to try some stupid heuristic",
    2176                     "off",NAIVE);
     2352                    "off",NAIVE,7,1);
    21772353  parameters[numberParameters-1].append("on");
    21782354  parameters[numberParameters-1].append("both");
     
    21812357    (
    21822358     "Really silly stuff e.g. fix all integers with costs to zero!. \
    2183 Do options does heuristic before preprocessing"     );
     2359Doh option does heuristic before preprocessing"     );
    21842360#endif
    21852361#ifdef COIN_HAS_CLP
    21862362  parameters[numberParameters++]=
    21872363    CbcOrClpParam("netlib","Solve entire netlib test set",
    2188                   NETLIB_EITHER,3);
     2364                  NETLIB_EITHER,3,1);
    21892365  parameters[numberParameters-1].setLonghelp
    21902366    (
     
    21922368The user can set options before e.g. clp -presolve off -netlib"
    21932369     );
    2194 #ifdef REAL_BARRIER
    21952370  parameters[numberParameters++]=
    21962371    CbcOrClpParam("netlibB!arrier","Solve entire netlib test set with barrier",
    2197                   NETLIB_BARRIER,3);
     2372                  NETLIB_BARRIER,3,1);
    21982373  parameters[numberParameters-1].setLonghelp
    21992374    (
     
    22012376The user can set options before e.g. clp -kkt on -netlib"
    22022377     );
    2203 #endif
    22042378  parameters[numberParameters++]=
    22052379    CbcOrClpParam("netlibD!ual","Solve entire netlib test set (dual)",
    2206                   NETLIB_DUAL,3);
     2380                  NETLIB_DUAL,3,1);
    22072381  parameters[numberParameters-1].setLonghelp
    22082382    (
     
    22122386  parameters[numberParameters++]=
    22132387    CbcOrClpParam("netlibP!rimal","Solve entire netlib test set (primal)",
    2214                   NETLIB_PRIMAL,3);
     2388                  NETLIB_PRIMAL,3,1);
    22152389  parameters[numberParameters-1].setLonghelp
    22162390    (
     
    22202394  parameters[numberParameters++]=
    22212395    CbcOrClpParam("netlibT!une","Solve entire netlib test set with 'best' algorithm",
    2222                   NETLIB_TUNE,3);
     2396                  NETLIB_TUNE,3,1);
    22232397  parameters[numberParameters-1].setLonghelp
    22242398    (
     
    22302404  parameters[numberParameters++]=
    22312405    CbcOrClpParam("network","Tries to make network matrix",
    2232                   NETWORK,7,false);
     2406                  NETWORK,7,0);
    22332407  parameters[numberParameters-1].setLonghelp
    22342408    (
     
    22572431  parameters[numberParameters++]=
    22582432    CbcOrClpParam("numberA!nalyze","Number of analysis iterations",
    2259                   -COIN_INT_MAX,COIN_INT_MAX,NUMBERANALYZE,false);
     2433                  -COIN_INT_MAX,COIN_INT_MAX,NUMBERANALYZE,0);
    22602434  parameters[numberParameters-1].setLonghelp
    22612435    (
     
    22662440  parameters[numberParameters++]=
    22672441    CbcOrClpParam("objective!Scale","Scale factor to apply to objective",
    2268                   -1.0e20,1.0e20,OBJSCALE,false);
     2442                  -1.0e20,1.0e20,OBJSCALE,1);
    22692443  parameters[numberParameters-1].setLonghelp
    22702444    (
    22712445     "If the objective function has some very large values, you may wish to scale them\
    2272  internally by this amount.  It can also be set by autoscale.  It is applied after scaling"
     2446 internally by this amount.  It can also be set by autoscale.  It is applied after scaling.  You are unlikely to need this."
    22732447     );
    22742448  parameters[numberParameters-1].setDoubleValue(1.0);
     
    22772451  parameters[numberParameters++]=
    22782452    CbcOrClpParam("outDup!licates","takes duplicate rows etc out of integer model",
    2279                   OUTDUPROWS,7,false);
     2453                  OUTDUPROWS,7,0);
    22802454#endif
    22812455  parameters[numberParameters++]=
     
    23122486  parameters[numberParameters++]=
    23132487    CbcOrClpParam("passP!resolve","How many passes in presolve",
    2314                   -200,100,PRESOLVEPASS,false);
     2488                  -200,100,PRESOLVEPASS,1);
    23152489  parameters[numberParameters-1].setLonghelp
    23162490    (
     
    23322506  parameters[numberParameters++]=
    23332507    CbcOrClpParam("pertV!alue","Method of perturbation",
    2334                   -5000,102,PERTVALUE,false);
     2508                  -5000,102,PERTVALUE,1);
    23352509  parameters[numberParameters++]=
    23362510    CbcOrClpParam("perturb!ation","Whether to perturb problem",
     
    23462520  parameters[numberParameters++]=
    23472521    CbcOrClpParam("PFI","Whether to use Product Form of Inverse in simplex",
    2348                   "off",PFI,7,false);
     2522                  "off",PFI,7,0);
    23492523  parameters[numberParameters-1].append("on");
    23502524  parameters[numberParameters-1].setLonghelp
     
    23552529#ifdef COIN_HAS_CBC
    23562530  parameters[numberParameters++]=
    2357       CbcOrClpParam("pivot!AndFix","Whether to try Pivot and Fix heuristic",
    2358                     "off",PIVOTANDFIX);
     2531      CbcOrClpParam("pivotAndC!omplement","Whether to try Pivot and Complement heuristic",
     2532                    "off",PIVOTANDCOMPLEMENT);
    23592533  parameters[numberParameters-1].append("on");
    23602534  parameters[numberParameters-1].append("both");
     
    23632537    (
    23642538     "stuff needed. \
    2365 Do options does heuristic before preprocessing"     );
     2539Doh option does heuristic before preprocessing"     );
     2540  parameters[numberParameters++]=
     2541      CbcOrClpParam("pivotAndF!ix","Whether to try Pivot and Fix heuristic",
     2542                    "off",PIVOTANDFIX);
     2543  parameters[numberParameters-1].append("on");
     2544  parameters[numberParameters-1].append("both");
     2545  parameters[numberParameters-1].append("before");
     2546  parameters[numberParameters-1].setLonghelp
     2547    (
     2548     "stuff needed. \
     2549Doh option does heuristic before preprocessing"     );
    23662550#endif
    23672551#ifdef COIN_HAS_CLP
    23682552  parameters[numberParameters++]=
    23692553    CbcOrClpParam("plus!Minus","Tries to make +- 1 matrix",
    2370                   PLUSMINUS,7,false);
     2554                  PLUSMINUS,7,0);
    23712555  parameters[numberParameters-1].setLonghelp
    23722556    (
     
    23772561  parameters[numberParameters++]=
    23782562    CbcOrClpParam("pO!ptions","Dubious print options",
    2379                   0,COIN_INT_MAX,PRINTOPTIONS,false);
     2563                  0,COIN_INT_MAX,PRINTOPTIONS,1);
    23802564  parameters[numberParameters-1].setIntValue(0);
    23812565  parameters[numberParameters-1].setLonghelp
     
    23852569  parameters[numberParameters++]=
    23862570    CbcOrClpParam("preO!pt","Presolve options",
    2387                   0,COIN_INT_MAX,PRESOLVEOPTIONS,false);
     2571                  0,COIN_INT_MAX,PRESOLVEOPTIONS,0);
    23882572#endif
    23892573  parameters[numberParameters++]=
     
    24362620  parameters[numberParameters++]=
    24372621    CbcOrClpParam("primalP!ivot","Primal pivot choice algorithm",
    2438                   "auto!matic",PRIMALPIVOT);
     2622                  "auto!matic",PRIMALPIVOT,7,1);
    24392623  parameters[numberParameters-1].append("exa!ct");
    24402624  parameters[numberParameters-1].append("dant!zig");
     
    25362720  parameters[numberParameters-1].append("forceOn");
    25372721  parameters[numberParameters-1].append("onglobal");
    2538   parameters[numberParameters-1].append("rootglobal");
     2722  parameters[numberParameters-1].append("forceonglobal");
    25392723  parameters[numberParameters-1].append("forceOnBut");
    25402724  parameters[numberParameters-1].append("forceOnStrong");
     
    25482732     );
    25492733  parameters[numberParameters++]=
     2734    CbcOrClpParam("pumpC!utoff","Fake cutoff for use in feasibility pump",
     2735                  -COIN_DBL_MAX,COIN_DBL_MAX,FAKECUTOFF);
     2736  parameters[numberParameters-1].setDoubleValue(0.0);
     2737  parameters[numberParameters-1].setLonghelp
     2738    (
     2739     "0.0 off - otherwise add a constraint forcing objective below this value\
     2740 in feasibility pump"
     2741     );
     2742  parameters[numberParameters++]=
     2743    CbcOrClpParam("pumpI!ncrement","Fake increment for use in feasibility pump",
     2744                  -COIN_DBL_MAX,COIN_DBL_MAX,FAKEINCREMENT,1);
     2745  parameters[numberParameters-1].setDoubleValue(0.0);
     2746  parameters[numberParameters-1].setLonghelp
     2747    (
     2748     "0.0 off - otherwise use as absolute increment to cutoff \
     2749when solution found in feasibility pump"
     2750     );
     2751  parameters[numberParameters++]=
    25502752    CbcOrClpParam("pumpT!une","Dubious ideas for feasibility pump",
    25512753                  0,100000000,FPUMPTUNE);
     
    25532755    (
    25542756     "This fine tunes Feasibility Pump \n\
     2757\t>=10000000 use as objective weight switch\n\
    25552758\t>=1000000 use as accumulate switch\n\
    25562759\t>=1000 use index+1 as number of large loops\n\
    2557 \t==100 use objvalue +0.05*fabs(objvalue) as cutoff OR dextra1 if set\n\
    2558 \t%100 == 10,20 etc for experimentation\n\
    2559 \t1 == fix ints at bounds, 2 fix all integral ints, 3 and continuous at bounds"
     2760\t==100 use objvalue +0.05*fabs(objvalue) as cutoff OR fakeCutoff if set\n\
     2761\t%100 == 10,20 affects how each solve is done\n\
     2762\t1 == fix ints at bounds, 2 fix all integral ints, 3 and continuous at bounds. \
     2763If accumulate is on then after a major pass, variables which have not moved \
     2764are fixed and a small branch and bound is tried."
    25602765     );
    25612766  parameters[numberParameters-1].setIntValue(0);
     
    25782783    (
    25792784     "stuff needed. \
    2580 Do options does heuristic before preprocessing"     );
     2785Doh option does heuristic before preprocessing"     );
    25812786  parameters[numberParameters++]=
    25822787    CbcOrClpParam("ratio!Gap","Stop when gap between best possible and \
     
    25922797  parameters[numberParameters++]=
    25932798    CbcOrClpParam("readS!tored","Import stored cuts from file",
    2594                   STOREDFILE,3,false);
     2799                  STOREDFILE,3,0);
    25952800#endif
    25962801#ifdef COIN_HAS_CLP
    25972802  parameters[numberParameters++]=
    25982803    CbcOrClpParam("reallyO!bjectiveScale","Scale factor to apply to objective in place",
    2599                   -1.0e20,1.0e20,OBJSCALE2,false);
     2804                  -1.0e20,1.0e20,OBJSCALE2,0);
    26002805  parameters[numberParameters-1].setLonghelp
    26012806    (
     
    26052810  parameters[numberParameters++]=
    26062811    CbcOrClpParam("reallyS!cale","Scales model in place",
    2607                   REALLY_SCALE,7,false);
     2812                  REALLY_SCALE,7,0);
    26082813#endif
    26092814#ifdef COIN_HAS_CBC
     
    26362841  parameters[numberParameters++]=
    26372842    CbcOrClpParam("restore!Model","Restore model from binary file",
    2638                   RESTORE);
     2843                  RESTORE,7,1);
    26392844  parameters[numberParameters-1].setLonghelp
    26402845    (
     
    26452850  parameters[numberParameters++]=
    26462851    CbcOrClpParam("reverse","Reverses sign of objective",
    2647                   REVERSE,7,false);
     2852                  REVERSE,7,0);
    26482853  parameters[numberParameters-1].setLonghelp
    26492854    (
     
    26522857  parameters[numberParameters++]=
    26532858    CbcOrClpParam("rhs!Scale","Scale factor to apply to rhs and bounds",
    2654                   -1.0e20,1.0e20,RHSSCALE,false);
     2859                  -1.0e20,1.0e20,RHSSCALE,0);
    26552860  parameters[numberParameters-1].setLonghelp
    26562861    (
    26572862     "If the rhs or bounds have some very large meaningful values, you may wish to scale them\
    2658  internally by this amount.  It can also be set by autoscale"
     2863 internally by this amount.  It can also be set by autoscale.  This should not be needed."
    26592864     );
    26602865  parameters[numberParameters-1].setDoubleValue(1.0);
     
    26732878    (
    26742879     "This switches on Relaxation enforced neighborhood Search. \
    2675 on just does feasibility pump \
     2880on just does 50 nodes \
    26762881200 or 1000 does that many nodes. \
    2677 Do options does heuristic before preprocessing"     );
     2882Doh option does heuristic before preprocessing"     );
    26782883  parameters[numberParameters++]=
    26792884      CbcOrClpParam("Rins","Whether to try Relaxed Induced Neighborhood Search",
     
    26862891    (
    26872892     "This switches on Relaxed induced neighborhood Search. \
    2688 Do options does heuristic before preprocessing"     );
     2893Doh option does heuristic before preprocessing"     );
    26892894  parameters[numberParameters++]=
    26902895    CbcOrClpParam("round!ingHeuristic","Whether to use Rounding heuristic",
     
    27042909  parameters[numberParameters++]=
    27052910    CbcOrClpParam("saveM!odel","Save model to binary file",
    2706                   SAVE);
     2911                  SAVE,7,1);
    27072912  parameters[numberParameters-1].setLonghelp
    27082913    (
     
    27302935  parameters[numberParameters-1].append("geo!metric");
    27312936  parameters[numberParameters-1].append("auto!matic");
     2937  parameters[numberParameters-1].append("dynamic");
     2938  parameters[numberParameters-1].append("rows!only");
    27322939  parameters[numberParameters-1].setLonghelp
    27332940    (
     
    27452952    (
    27462953     "After this many seconds clp will act as if maximum iterations had been reached \
    2747 (if value >=0).  \
    2748 In this program it is really only useful for testing but the library function\n\
    2749       \tsetMaximumSeconds(value)\n can be useful."
     2954(if value >=0)."
    27502955     );
    27512956#else
     
    27602965  parameters[numberParameters++]=
    27612966    CbcOrClpParam("sleep","for debug",
    2762                   DUMMY,7,false);
     2967                  DUMMY,7,0);
    27632968  parameters[numberParameters-1].setLonghelp
    27642969    (
     
    27682973  parameters[numberParameters++]=
    27692974    CbcOrClpParam("slp!Value","Number of slp passes before primal",
    2770                   -1,50000,SLPVALUE,false);
     2975                  -1,50000,SLPVALUE,1);
    27712976  parameters[numberParameters-1].setLonghelp
    27722977    (
     
    27772982  parameters[numberParameters++]=
    27782983    CbcOrClpParam("small!Factorization","Whether to use small factorization",
    2779                   -1,10000,SMALLFACT,false);
     2984                  -1,10000,SMALLFACT,1);
    27802985  parameters[numberParameters-1].setLonghelp
    27812986    (
     
    28153020     );
    28163021  parameters[numberParameters++]=
    2817     CbcOrClpParam("slog!Level","Level of detail in Solver output",
     3022    CbcOrClpParam("slog!Level","Level of detail in (LP) Solver output",
    28183023                  -1,63,SOLVERLOGLEVEL);
    28193024  parameters[numberParameters-1].setLonghelp
    28203025    (
    28213026     "If 0 then there should be no output in normal circumstances.  1 is probably the best\
    2822  value for most uses, while 2 and 3 give more information."
     3027 value for most uses, while 2 and 3 give more information.  This parameter is only used inside MIP - for Clp use 'log'"
    28233028     );
    28243029#else
     
    28363041  parameters[numberParameters++]=
    28373042    CbcOrClpParam("spars!eFactor","Whether factorization treated as sparse",
    2838                   "on",SPARSEFACTOR,7,false);
     3043                  "on",SPARSEFACTOR,7,0);
    28393044  parameters[numberParameters-1].append("off");
    28403045  parameters[numberParameters++]=
    28413046    CbcOrClpParam("special!Options","Dubious options for Simplex - see ClpSimplex.hpp",
    2842                   0,COIN_INT_MAX,SPECIALOPTIONS,false);
     3047                  0,COIN_INT_MAX,SPECIALOPTIONS,0);
    28433048  parameters[numberParameters++]=
    28443049    CbcOrClpParam("sprint!Crash","Whether to try sprint crash",
     
    28833088     );
    28843089  parameters[numberParameters-1].setIntValue(1);
     3090#ifdef CBC_KEEP_DEPRECATED
    28853091  parameters[numberParameters++]=
    28863092    CbcOrClpParam("strengthen","Create strengthened problem",
     
    28903096     "This creates a new problem by applying the root node cuts.  All tight constraints \
    28913097will be in resulting problem"
    2892      );
     3098     );
     3099#endif
    28933100  parameters[numberParameters++]=
    28943101    CbcOrClpParam("strong!Branching","Number of variables to look at in strong branching",
     
    29053112  parameters[numberParameters++]=
    29063113    CbcOrClpParam("subs!titution","How long a column to substitute for in presolve",
    2907                   0,10000,SUBSTITUTION,false);
     3114                  0,10000,SUBSTITUTION,0);
    29083115  parameters[numberParameters-1].setLonghelp
    29093116    (
     
    29163123  parameters[numberParameters++]=
    29173124    CbcOrClpParam("testO!si","Test OsiObject stuff",
    2918                   -1,COIN_INT_MAX,TESTOSI,false);
     3125                  -1,COIN_INT_MAX,TESTOSI,0);
    29193126#endif
    29203127#ifdef CBC_THREAD
    29213128  parameters[numberParameters++]=
    29223129    CbcOrClpParam("thread!s","Number of threads to try and use",
    2923                   -100,100000,THREADS,false);
     3130                  -100,100000,THREADS,1);
    29243131  parameters[numberParameters-1].setLonghelp
    29253132    (
     
    29343141    CbcOrClpParam("tighten!Factor","Tighten bounds using this times largest \
    29353142activity at continuous solution",
    2936                   1.0e-3,1.0e20,TIGHTENFACTOR,false);
     3143                  1.0e-3,1.0e20,TIGHTENFACTOR,0);
    29373144  parameters[numberParameters-1].setLonghelp
    29383145    (
     
    29433150  parameters[numberParameters++]=
    29443151    CbcOrClpParam("tightLP","Poor person's preSolve for now",
    2945                   TIGHTEN,7,false);
     3152                  TIGHTEN,7,0);
    29463153#endif
    29473154#ifdef COIN_HAS_CBC
     
    29583165  parameters[numberParameters++]=
    29593166    CbcOrClpParam("tune!PreProcess","Dubious tuning parameters",
    2960                   0,20000000,PROCESSTUNE,false);
     3167                  0,20000000,PROCESSTUNE,1);
    29613168  parameters[numberParameters-1].setLonghelp
    29623169    (
     
    29773184  parameters[numberParameters-1].append("forceOn");
    29783185  parameters[numberParameters-1].append("onglobal");
    2979   parameters[numberParameters-1].append("rootglobal");
     3186  parameters[numberParameters-1].append("forceandglobal");
    29803187  parameters[numberParameters-1].append("forceLongOn");
    29813188  parameters[numberParameters-1].setLonghelp
     
    29873194  parameters[numberParameters++]=
    29883195    CbcOrClpParam("unitTest","Do unit test",
    2989                   UNITTEST,3);
     3196                  UNITTEST,3,1);
    29903197  parameters[numberParameters-1].setLonghelp
    29913198    (
     
    29943201  parameters[numberParameters++]=
    29953202    CbcOrClpParam("userClp","Hand coded Clp stuff",
    2996                   USERCLP);
     3203                  USERCLP,0,0);
    29973204  parameters[numberParameters-1].setLonghelp
    29983205    (
     
    30033210  parameters[numberParameters++]=
    30043211    CbcOrClpParam("userCbc","Hand coded Cbc stuff",
    3005                   USERCBC);
     3212                  USERCBC,0,0);
    30063213  parameters[numberParameters-1].setLonghelp
    30073214    (
    30083215     "There are times e.g. when using AMPL interface when you may wish to do something unusual.  \
    3009 Look for USERCBC in main driver and modify sample code."
    3010      );
     3216Look for USERCBC in main driver and modify sample code. \
     3217It is possible you can get same effect by using example driver4.cpp."
     3218     );
     3219  parameters[numberParameters++]=
     3220      CbcOrClpParam("Vnd!VariableNeighborhoodSearch","Whether to try Variable Neighborhood Search",
     3221                    "off",VND);
     3222    parameters[numberParameters-1].append("on");
     3223    parameters[numberParameters-1].append("both");
     3224    parameters[numberParameters-1].append("before");
     3225    parameters[numberParameters-1].append("intree");
     3226  parameters[numberParameters-1].setLonghelp
     3227    (
     3228     "This switches on variable neighborhood Search. \
     3229Doh option does heuristic before preprocessing"     );
    30113230#endif
    30123231  parameters[numberParameters++]=
    30133232    CbcOrClpParam("vector","Whether to use vector? Form of matrix in simplex",
    3014                   "off",VECTOR,7,false);
    3015   parameters[numberParameters-1].append("on");
    3016   parameters[numberParameters-1].setLonghelp
    3017     (
    3018      "If this is on and ClpPackedMatrix uses extra column copy in odd format."
     3233                  "off",VECTOR,7,0);
     3234  parameters[numberParameters-1].append("on");
     3235  parameters[numberParameters-1].setLonghelp
     3236    (
     3237     "If this is on ClpPackedMatrix uses extra column copy in odd format."
    30193238     );
    30203239  parameters[numberParameters++]=
    30213240    CbcOrClpParam("verbose","Switches on longer help on single ?",
    3022                   0,15,VERBOSE,false);
     3241                  0,31,VERBOSE,0);
    30233242  parameters[numberParameters-1].setLonghelp
    30243243    (
     
    30293248  parameters[numberParameters++]=
    30303249    CbcOrClpParam("vub!heuristic","Type of vub heuristic",
    3031                   -2,20,VUBTRY,false);
     3250                  -2,20,VUBTRY,0);
    30323251  parameters[numberParameters-1].setLonghelp
    30333252    (
     
    30353254     );
    30363255  parameters[numberParameters-1].setIntValue(-1);
     3256#ifdef ZERO_HALF_CUTS
    30373257  parameters[numberParameters++]=
    30383258    CbcOrClpParam("zero!HalfCuts","Whether to use zero half cuts",
     
    30433263  parameters[numberParameters-1].append("forceOn");
    30443264  parameters[numberParameters-1].append("onglobal");
    3045   parameters[numberParameters-1].append("rootglobal");
    3046   parameters[numberParameters-1].setLonghelp
    3047     (
    3048      "This switches on knapsack cuts (either at root or in entire tree) \
     3265  parameters[numberParameters-1].setLonghelp
     3266    (
     3267     "This switches on zero-half cuts (either at root or in entire tree) \
    30493268See branchAndCut for information on options."
    3050      );
     3269     );
     3270#endif
    30513271#endif
    30523272  assert(numberParameters<CBCMAXPARAMETERS);
     
    30783298    int numberColumnsFile;
    30793299    double objectiveValue;
    3080     fread(&numberRowsFile,sizeof(int),1,fp);
    3081     fread(&numberColumnsFile,sizeof(int),1,fp);
    3082     fread(&objectiveValue,sizeof(double),1,fp);
     3300    int nRead;
     3301    nRead=fread(&numberRowsFile,sizeof(int),1,fp);
     3302    if (nRead!=1)
     3303      throw("Error in fread");
     3304    nRead=fread(&numberColumnsFile,sizeof(int),1,fp);
     3305    if (nRead!=1)
     3306      throw("Error in fread");
     3307    nRead=fread(&objectiveValue,sizeof(double),1,fp);
     3308    if (nRead!=1)
     3309      throw("Error in fread");
    30833310    double * dualRowSolution = lpSolver->dualRowSolution();
    30843311    double * primalRowSolution = lpSolver->primalRowSolution();
     
    31033330      lpSolver->setObjectiveValue(objectiveValue);
    31043331      if (numberRows==numberRowsFile&&numberColumns==numberColumnsFile) {
    3105         fread(primalRowSolution,sizeof(double),numberRows,fp);
    3106         fread(dualRowSolution,sizeof(double),numberRows,fp);
    3107         fread(primalColumnSolution,sizeof(double),numberColumns,fp);
    3108         fread(dualColumnSolution,sizeof(double),numberColumns,fp);
     3332        nRead=fread(primalRowSolution,sizeof(double),numberRows,fp);
     3333        if (nRead!=numberRows)
     3334          throw("Error in fread");
     3335        nRead=fread(dualRowSolution,sizeof(double),numberRows,fp);
     3336        if (nRead!=numberRows)
     3337          throw("Error in fread");
     3338        nRead=fread(primalColumnSolution,sizeof(double),numberColumns,fp);
     3339        if (nRead!=numberColumns)
     3340          throw("Error in fread");
     3341        nRead=fread(dualColumnSolution,sizeof(double),numberColumns,fp);
     3342        if (nRead!=numberColumns)
     3343          throw("Error in fread");
    31093344      } else {
    31103345        std::cout<<"Mismatch on rows and/or columns - truncating"<<std::endl;
    31113346        double * temp = new double [CoinMax(numberRowsFile,numberColumnsFile)];
    3112         fread(temp,sizeof(double),numberRowsFile,fp);
    3113  CoinMemcpyN(temp,numberRows,primalRowSolution);
    3114         fread(temp,sizeof(double),numberRowsFile,fp);
    3115  CoinMemcpyN(temp,numberRows,dualRowSolution);
    3116         fread(temp,sizeof(double),numberColumnsFile,fp);
    3117  CoinMemcpyN(temp,numberColumns,primalColumnSolution);
    3118         fread(temp,sizeof(double),numberColumnsFile,fp);
    3119  CoinMemcpyN(temp,numberColumns,dualColumnSolution);
     3347        nRead=fread(temp,sizeof(double),numberRowsFile,fp);
     3348        if (nRead!=numberRowsFile)
     3349          throw("Error in fread");
     3350        CoinMemcpyN(temp,numberRows,primalRowSolution);
     3351        nRead=fread(temp,sizeof(double),numberRowsFile,fp);
     3352        if (nRead!=numberRowsFile)
     3353          throw("Error in fread");
     3354        CoinMemcpyN(temp,numberRows,dualRowSolution);
     3355        nRead=fread(temp,sizeof(double),numberColumnsFile,fp);
     3356        if (nRead!=numberColumnsFile)
     3357          throw("Error in fread");
     3358        CoinMemcpyN(temp,numberColumns,primalColumnSolution);
     3359        nRead=fread(temp,sizeof(double),numberColumnsFile,fp);
     3360        if (nRead!=numberColumnsFile)
     3361          throw("Error in fread");
     3362        CoinMemcpyN(temp,numberColumns,dualColumnSolution);
    31203363        delete [] temp;
    31213364      }
     
    31773420    int numberColumns=lpSolver->numberColumns();
    31783421    double objectiveValue = lpSolver->objectiveValue();
    3179     fwrite(&numberRows,sizeof(int),1,fp);
    3180     fwrite(&numberColumns,sizeof(int),1,fp);
    3181     fwrite(&objectiveValue,sizeof(double),1,fp);
     3422    int nWrite;
     3423    nWrite=fwrite(&numberRows,sizeof(int),1,fp);
     3424    if (nWrite!=1)
     3425      throw("Error in fwrite");
     3426    nWrite=fwrite(&numberColumns,sizeof(int),1,fp);
     3427    if (nWrite!=1)
     3428      throw("Error in fwrite");
     3429    nWrite=fwrite(&objectiveValue,sizeof(double),1,fp);
     3430    if (nWrite!=1)
     3431      throw("Error in fwrite");
    31823432    double * dualRowSolution = lpSolver->dualRowSolution();
    31833433    double * primalRowSolution = lpSolver->primalRowSolution();
    3184     fwrite(primalRowSolution,sizeof(double),numberRows,fp);
    3185     fwrite(dualRowSolution,sizeof(double),numberRows,fp);
     3434    nWrite=fwrite(primalRowSolution,sizeof(double),numberRows,fp);
     3435    if (nWrite!=numberRows)
     3436      throw("Error in fwrite");
     3437    nWrite=fwrite(dualRowSolution,sizeof(double),numberRows,fp);
     3438    if (nWrite!=numberRows)
     3439      throw("Error in fwrite");
    31863440    double * dualColumnSolution = lpSolver->dualColumnSolution();
    31873441    double * primalColumnSolution = lpSolver->primalColumnSolution();
    3188     fwrite(primalColumnSolution,sizeof(double),numberColumns,fp);
    3189     fwrite(dualColumnSolution,sizeof(double),numberColumns,fp);
     3442    nWrite=fwrite(primalColumnSolution,sizeof(double),numberColumns,fp);
     3443    if (nWrite!=numberColumns)
     3444      throw("Error in fwrite");
     3445    nWrite=fwrite(dualColumnSolution,sizeof(double),numberColumns,fp);
     3446    if (nWrite!=numberColumns)
     3447      throw("Error in fwrite");
    31903448    fclose(fp);
    31913449  } else {
  • stable/1.11/Clp/src/CbcOrClpParam.hpp

    r1344 r1458  
     1
     2/* $Id$ */
    13// Copyright (C) 2002, International Business Machines
    24// Corporation and others.  All Rights Reserved.
     
    5355   
    5456    DJFIX = 81, TIGHTENFACTOR,PRESOLVETOLERANCE,OBJSCALE2,
    55     DEXTRA1, DEXTRA2, DEXTRA3, DEXTRA4, DEXTRA5,
     57    FAKEINCREMENT, FAKECUTOFF, ARTIFICIALCOST,DEXTRA3, SMALLBAB, DEXTRA4, DEXTRA5,
    5658
    5759    SOLVERLOGLEVEL=101,
     
    6466
    6567    STRONGBRANCHING=151,CUTDEPTH, MAXNODES,NUMBERBEFORE,NUMBERANALYZE,
    66     NUMBERMINI,MIPOPTIONS,MOREMIPOPTIONS,MAXHOTITS,FPUMPITS,MAXSOLS,
    67     FPUMPTUNE,TESTOSI,EXTRA1,EXTRA2,EXTRA3,EXTRA4,CUTPASSINTREE,
     68    MIPOPTIONS,MOREMIPOPTIONS,MAXHOTITS,FPUMPITS,MAXSOLS,
     69    FPUMPTUNE,TESTOSI,EXTRA1,EXTRA2,EXTRA3,EXTRA4,DEPTHMINIBAB,CUTPASSINTREE,
    6870    THREADS,CUTPASS,VUBTRY,DENSE,EXPERIMENT,DIVEOPT,STRATEGY,SMALLFACT,
    69     HOPTIONS,
     71    HOPTIONS,CUTLENGTH,FPUMPTUNE2,
    7072#ifdef COIN_HAS_CBC
    7173    LOGLEVEL ,
     
    7577    PRIMALPIVOT,PRESOLVE,CRASH,BIASLU,PERTURBATION,MESSAGES,AUTOSCALE,
    7678    CHOLESKY,KKT,BARRIERSCALE,GAMMA,CROSSOVER,PFI,INTPRINT,VECTOR,
     79    FACTORIZATION,ALLCOMMANDS,
    7780   
    7881    NODESTRATEGY = 251,BRANCHSTRATEGY,CUTSSTRATEGY,HEURISTICSTRATEGY,
     
    8285    LANDPCUTS,RINS,RESIDCUTS,RENS,DIVINGS,DIVINGC,DIVINGF,DIVINGG,DIVINGL,
    8386    DIVINGP,DIVINGV,DINS,PIVOTANDFIX,RANDROUND,NAIVE,ZEROHALFCUTS,CPX,
     87    CROSSOVER2,PIVOTANDCOMPLEMENT,VND,
    8488   
    8589    DIRECTORY=301,DIRSAMPLE,DIRNETLIB,DIRMIPLIB,IMPORT,EXPORT,RESTORE,SAVE,DUALSIMPLEX,PRIMALSIMPLEX,EITHERSIMPLEX,
     
    8791    TIGHTEN,FAKEBOUND,HELP,PLUSMINUS,NETWORK,ALLSLACK,REVERSE,BARRIER,NETLIB_BARRIER,NETLIB_TUNE,
    8892    REALLY_SCALE,BASISIN,BASISOUT,SOLVECONTINUOUS,CLEARCUTS,VERSION,STATISTICS,DEBUG,DUMMY,PRINTMASK,
    89     OUTDUPROWS,USERCLP,MODELIN,CSVSTATISTICS,STOREDFILE,
     93    OUTDUPROWS,USERCLP,MODELIN,CSVSTATISTICS,STOREDFILE,ENVIRONMENT,
    9094
    9195    BAB=351,MIPLIB,STRENGTHEN,PRIORITYIN,USERCBC,DOHEURISTIC,
     
    107111  CbcOrClpParam (  );
    108112  CbcOrClpParam (std::string name, std::string help,
    109            double lower, double upper, CbcOrClpParameterType type,bool display=true);
     113           double lower, double upper, CbcOrClpParameterType type,int display=2);
    110114  CbcOrClpParam (std::string name, std::string help,
    111            int lower, int upper, CbcOrClpParameterType type,bool display=true);
     115           int lower, int upper, CbcOrClpParameterType type,int display=2);
    112116  // Other strings will be added by insert
    113117  CbcOrClpParam (std::string name, std::string help, std::string firstValue,
    114            CbcOrClpParameterType type,int whereUsed=7,bool display=true);
     118           CbcOrClpParameterType type,int whereUsed=7,int display=2);
    115119  // Action
    116120  CbcOrClpParam (std::string name, std::string help,
    117            CbcOrClpParameterType type,int whereUsed=7,bool display=true);
     121           CbcOrClpParameterType type,int whereUsed=7,int display=2);
    118122  /// Copy constructor.
    119123  CbcOrClpParam(const CbcOrClpParam &);
     
    178182  /// Returns name which could match
    179183  std::string matchName (  ) const;
     184  /// Returns length of name for ptinting
     185  int lengthMatchName (  ) const;
    180186  /// Returns parameter option which matches (-1 if none)
    181187  int parameterOption ( std::string check ) const;
     
    212218  { return type_;}
    213219  /// whether to display
    214   inline bool displayThis() const
     220  inline int displayThis() const
    215221  { return display_;}
    216222  /// Set Long help
     
    265271  int currentKeyWord_;
    266272  /// Display on ?
    267   bool display_;
     273  int display_;
    268274  /// Integer parameter - current value
    269275  int intValue_;
  • stable/1.11/Clp/src/ClpCholeskyBase.cpp

    r1321 r1458  
     1/* $Id$ */
    12// Copyright (C) 2002, International Business Machines
    23// Corporation and others.  All Rights Reserved.
     4/*----------------------------------------------------------------------------*/
     5/*      Ordering code - courtesy of Anshul Gupta                              */
     6/*      (C) Copyright IBM Corporation 1997, 2009.  All Rights Reserved.       */
     7/*----------------------------------------------------------------------------*/
     8
     9/*  A compact no-frills Approximate Minimum Local Fill ordering code.
     10
     11    References:
     12
     13[1] Ordering Sparse Matrices Using Approximate Minimum Local Fill.
     14    Edward Rothberg, SGI Manuscript, April 1996.
     15[2] An Approximate Minimum Degree Ordering Algorithm.
     16    T. Davis, P. Amestoy, and I. Duff, TR-94-039, CIS Department,
     17    University of Florida, December 1994.
     18*/
     19/*----------------------------------------------------------------------------*/
     20
    321
    422#include "CoinPragma.hpp"
    523
    6 #include <iostream>
     24#include <iostream> 
    725
    826#include "ClpCholeskyBase.hpp"
     
    86104#else
    87105  // actually long double
    88   workDouble_ = (double *) ClpCopyOfArray((longWork *) rhs.workDouble_,numberRows_);
     106  workDouble_ = reinterpret_cast<double *> (ClpCopyOfArray(reinterpret_cast<CoinWorkDouble *> (rhs.workDouble_),numberRows_));
    89107#endif
    90108  link_ = ClpCopyOfArray(rhs.link_,numberRows_);
     
    171189#else
    172190    // actually long double
    173     workDouble_ = (double *) ClpCopyOfArray((longWork *) rhs.workDouble_,numberRows_);
     191    workDouble_ = reinterpret_cast<double *> (ClpCopyOfArray(reinterpret_cast<CoinWorkDouble *> (rhs.workDouble_),numberRows_));
    174192#endif
    175193    link_ = ClpCopyOfArray(rhs.link_,numberRows_);
     
    195213   region1 is rows+columns, region2 is rows */
    196214void
    197 ClpCholeskyBase::solveKKT (double * region1, double * region2, const double * diagonal,
    198                            double diagonalScaleFactor)
     215ClpCholeskyBase::solveKKT (CoinWorkDouble * region1, CoinWorkDouble * region2, const CoinWorkDouble * diagonal,
     216                           CoinWorkDouble diagonalScaleFactor)
    199217{
    200218  if (!doKKT_) {
     
    202220    int numberColumns = model_->numberColumns();
    203221    int numberTotal = numberRows_+numberColumns;
    204     double * region1Save = new double[numberTotal];
     222    CoinWorkDouble * region1Save = new CoinWorkDouble[numberTotal];
    205223    for (iColumn=0;iColumn<numberTotal;iColumn++) {
    206224      region1[iColumn] *= diagonal[iColumn];
     
    209227    multiplyAdd(region1+numberColumns,numberRows_,-1.0,region2,1.0);
    210228    model_->clpMatrix()->times(1.0,region1,region2);
    211     double maximumRHS = maximumAbsElement(region2,numberRows_);
    212     double scale=1.0;
    213     double unscale=1.0;
     229    CoinWorkDouble maximumRHS = maximumAbsElement(region2,numberRows_);
     230    CoinWorkDouble scale=1.0;
     231    CoinWorkDouble unscale=1.0;
    214232    if (maximumRHS>1.0e-30) {
    215233      if (maximumRHS<=0.5) {
    216         double factor=2.0;
     234        CoinWorkDouble factor=2.0;
    217235        while (maximumRHS<=0.5) {
    218236          maximumRHS*=factor;
     
    220238        } /* endwhile */
    221239      } else if (maximumRHS>=2.0&&maximumRHS<=COIN_DBL_MAX) {
    222         double factor=0.5;
     240        CoinWorkDouble factor=0.5;
    223241        while (maximumRHS>=2.0) {
    224242          maximumRHS*=factor;
     
    246264    int numberColumns = model_->numberColumns();
    247265    int numberTotal = numberColumns + numberRowsModel;
    248     double * array = new double [numberRows_];
     266    CoinWorkDouble * array = new CoinWorkDouble [numberRows_];
    249267    CoinMemcpyN(region1,numberTotal,array);
    250268    CoinMemcpyN(region2,numberRowsModel,array+numberTotal);
     
    253271    int iRow;
    254272    for (iRow=0;iRow<numberTotal;iRow++) {
    255       if (rowsDropped_[iRow]&&fabs(array[iRow])>1.0e-8) {
     273      if (rowsDropped_[iRow]&&CoinAbs(array[iRow])>1.0e-8) {
    256274        printf("row region1 %d dropped %g\n",iRow,array[iRow]);
    257275      }
    258276    }
    259277    for (;iRow<numberRows_;iRow++) {
    260       if (rowsDropped_[iRow]&&fabs(array[iRow])>1.0e-8) {
     278      if (rowsDropped_[iRow]&&CoinAbs(array[iRow])>1.0e-8) {
    261279        printf("row region2 %d dropped %g\n",iRow,array[iRow]);
    262280      }
     
    621639{
    622640  model_=model;
     641#define BASE_ORDER 2
     642#if BASE_ORDER>0
     643  if (!doKKT_&&model_->numberRows()>6) {
     644    if (preOrder(false,true,false))
     645      return -1;
     646    //rowsDropped_ = new char [numberRows_];
     647    numberRowsDropped_=0;
     648    memset(rowsDropped_,0,numberRows_);
     649    //rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
     650    // approximate minimum degree
     651    return orderAMD();
     652  }
     653#endif
    623654  int numberRowsModel = model_->numberRows();
    624655  int numberColumns = model_->numberColumns();
     
    636667  rowsDropped_ = new char [numberRows_];
    637668  numberRowsDropped_=0;
     669  memset(rowsDropped_,0,numberRows_);
    638670  rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
    639671  const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
     
    749781  delete [] count;
    750782  permuteInverse_ = new int [numberRows_];
    751   memset(rowsDropped_,0,numberRows_);
    752783  for (iRow=0;iRow<numberRows_;iRow++) {
    753784    //permute_[iRow]=iRow; // force no permute
     
    757788  }
    758789  return 0;
     790}
     791#if BASE_ORDER==1
     792/* Orders rows and saves pointer to matrix.and model */
     793int
     794ClpCholeskyBase::orderAMD()
     795{
     796  permuteInverse_ = new int [numberRows_];
     797  permute_ = new int[numberRows_];
     798  // Do ordering
     799  int returnCode=0;
     800  // get more space and full matrix
     801  int space = 6*sizeFactor_+100000;
     802  int * temp = new int [space];
     803  int * which = new int[2*numberRows_];
     804  CoinBigIndex * tempStart = new CoinBigIndex [numberRows_+1];
     805  memset(which,0,numberRows_*sizeof(int));
     806  for (int iRow=0;iRow<numberRows_;iRow++) {
     807    which[iRow]+= choleskyStart_[iRow+1]-choleskyStart_[iRow]-1;
     808    for (CoinBigIndex j=choleskyStart_[iRow]+1;j<choleskyStart_[iRow+1];j++) {
     809      int jRow=choleskyRow_[j];
     810      which[jRow]++;
     811    }
     812  }
     813  CoinBigIndex sizeFactor =0;
     814  for (int iRow=0;iRow<numberRows_;iRow++) {
     815    int length = which[iRow];
     816    permute_[iRow]=length;
     817    tempStart[iRow]=sizeFactor;
     818    which[iRow]=sizeFactor;
     819    sizeFactor += length;
     820  }
     821  for (int iRow=0;iRow<numberRows_;iRow++) {
     822    assert (choleskyRow_[choleskyStart_[iRow]]==iRow);
     823    for (CoinBigIndex j=choleskyStart_[iRow]+1;j<choleskyStart_[iRow+1];j++) {
     824      int jRow=choleskyRow_[j];
     825      int put = which[iRow];
     826      temp[put]=jRow;
     827      which[iRow]++;
     828      put = which[jRow];
     829      temp[put]=iRow;
     830      which[jRow]++;
     831    }
     832  }
     833  for (int iRow=1;iRow<numberRows_;iRow++)
     834    assert (which[iRow-1]==tempStart[iRow]);
     835  CoinBigIndex lastSpace=sizeFactor;
     836  delete [] choleskyRow_;
     837  choleskyRow_=temp;
     838  delete [] choleskyStart_;
     839  choleskyStart_=tempStart;
     840  // linked lists of sizes and lengths
     841  int * first = new int [numberRows_];
     842  int * next = new int [numberRows_];
     843  int * previous = new int [numberRows_];
     844  char * mark = new char[numberRows_];
     845  memset(mark,0,numberRows_);
     846  CoinBigIndex * sort = new CoinBigIndex [numberRows_];
     847  int * order = new int [numberRows_];
     848  for (int iRow=0;iRow<numberRows_;iRow++) {
     849    first[iRow]=-1;
     850    next[iRow]=-1;
     851    previous[iRow]=-1;
     852    permuteInverse_[iRow]=-1;
     853  }
     854  int large = 1000+2*numberRows_;
     855  for (int iRow=0;iRow<numberRows_;iRow++) {
     856    int n = permute_[iRow];
     857    if (first[n]<0) {
     858      first[n]=iRow;
     859      previous[iRow]=n+large;
     860      next[iRow]=n+2*large;
     861    } else {
     862      int k=first[n];
     863      assert (k<numberRows_);
     864      first[n]=iRow;
     865      previous[iRow]=n+large;
     866      assert (previous[k]==n+large);
     867      next[iRow]=k;
     868      previous[k]=iRow;
     869    }
     870  }
     871  int smallest=0;
     872  int done=0;
     873  int numberCompressions=0;
     874  int numberExpansions=0;
     875  while (smallest<numberRows_) {
     876    if (first[smallest]<0||first[smallest]>numberRows_) {
     877      smallest++;
     878      continue;
     879    }
     880    int iRow=first[smallest];
     881    if (false) {
     882      int kRow=-1;
     883      int ss=999999;
     884      for (int jRow=numberRows_-1;jRow>=0;jRow--) {
     885        if (permuteInverse_[jRow]<0) {
     886          int length = permute_[jRow];
     887          assert (length>0);
     888          for (CoinBigIndex j=choleskyStart_[jRow];
     889               j<choleskyStart_[jRow]+length;j++) {
     890            int jjRow=choleskyRow_[j];
     891            assert (permuteInverse_[jjRow]<0);
     892          }
     893          if (length<ss) {
     894            ss=length;
     895            kRow=jRow;
     896          }
     897        }
     898      }
     899      assert (smallest==ss);
     900      printf("list chose %d - full chose %d - length %d\n",
     901             iRow,kRow,ss);
     902    }
     903    int kNext = next[iRow];
     904    first[smallest]=kNext;
     905    if (kNext<numberRows_)
     906      previous[kNext]=previous[iRow];
     907    previous[iRow]=-1;
     908    next[iRow]=-1;
     909    permuteInverse_[iRow]=done;
     910    done++;
     911    // Now add edges
     912    CoinBigIndex start=choleskyStart_[iRow];
     913    CoinBigIndex end=choleskyStart_[iRow]+permute_[iRow];
     914    int nSave=0;
     915    for (CoinBigIndex k=start;k<end;k++) {
     916      int kRow=choleskyRow_[k];
     917      assert (permuteInverse_[kRow]<0);
     918      //if (permuteInverse_[kRow]<0)
     919      which[nSave++]=kRow;
     920    }
     921    for (int i=0;i<nSave;i++) {
     922      int kRow = which[i];
     923      int length = permute_[kRow];
     924      CoinBigIndex start=choleskyStart_[kRow];
     925      CoinBigIndex end=choleskyStart_[kRow]+length;
     926      for (CoinBigIndex j=start;j<end;j++) {
     927        int jRow=choleskyRow_[j];
     928        mark[jRow]=1;
     929      }
     930      mark[kRow]=1;
     931      int n=nSave;
     932      for (int j=0;j<nSave;j++) {
     933        int jRow=which[j];
     934        if (!mark[jRow]) {
     935          which[n++]=jRow;
     936        }
     937      }
     938      for (CoinBigIndex j=start;j<end;j++) {
     939        int jRow=choleskyRow_[j];
     940        mark[jRow]=0;
     941      }
     942      mark[kRow]=0;
     943      if (n>nSave) {
     944        bool inPlace = (n-nSave==1);
     945        if (!inPlace) {
     946          // extend
     947          int length = n-nSave+end-start;
     948          if (length+lastSpace>space) {
     949            // need to compress
     950            numberCompressions++;
     951            int newN=0;
     952            for (int iRow=0;iRow<numberRows_;iRow++) {
     953              CoinBigIndex start=choleskyStart_[iRow];
     954              if (permuteInverse_[iRow]<0) {
     955                sort[newN]=start;
     956                order[newN++]=iRow;
     957              } else {
     958                choleskyStart_[iRow]=-1;
     959                permute_[iRow]=0;
     960              }
     961            }
     962            CoinSort_2(sort,sort+newN,order);
     963            sizeFactor=0;
     964            for (int k=0;k<newN;k++) {
     965              int iRow=order[k];
     966              int length = permute_[iRow];
     967              CoinBigIndex start=choleskyStart_[iRow];
     968              choleskyStart_[iRow]=sizeFactor;
     969              for (int j=0;j<length;j++)
     970                choleskyRow_[sizeFactor+j]=choleskyRow_[start+j];
     971              sizeFactor += length;
     972            }
     973            lastSpace=sizeFactor;
     974            if ((sizeFactor+numberRows_+1000)*4>3*space) {
     975              // need to expand
     976              numberExpansions++;
     977              space = (3*space)/2;
     978              int * temp = new int [space];
     979              memcpy(temp,choleskyRow_,sizeFactor*sizeof(int));
     980              delete [] choleskyRow_;
     981              choleskyRow_=temp;
     982            }
     983          }
     984        }
     985        // Now add
     986        start=choleskyStart_[kRow];
     987        end=choleskyStart_[kRow]+permute_[kRow];
     988        if (!inPlace)
     989          choleskyStart_[kRow]=lastSpace;
     990        CoinBigIndex put = choleskyStart_[kRow];
     991        for (CoinBigIndex j=start;j<end;j++) {
     992          int jRow=choleskyRow_[j];
     993          if (permuteInverse_[jRow]<0)
     994            choleskyRow_[put++]=jRow;
     995        }
     996        for (int j=nSave;j<n;j++) {
     997          int jRow=which[j];
     998          choleskyRow_[put++]=jRow;
     999        }
     1000        if (!inPlace) {
     1001          permute_[kRow]=put-lastSpace;
     1002          lastSpace=put;
     1003        } else {
     1004          permute_[kRow]=put-choleskyStart_[kRow];
     1005        }
     1006      } else {
     1007        // pack down for new counts
     1008        CoinBigIndex put=start;
     1009        for (CoinBigIndex j=start;j<end;j++) {
     1010          int jRow=choleskyRow_[j];
     1011          if (permuteInverse_[jRow]<0)
     1012            choleskyRow_[put++]=jRow;
     1013        }
     1014        permute_[kRow]=put-start;
     1015      }
     1016      // take out
     1017      int iNext = next[kRow];
     1018      int iPrevious = previous[kRow];
     1019      if (iPrevious<numberRows_) {
     1020        next[iPrevious]=iNext;
     1021      } else {
     1022        assert (iPrevious==length+large);
     1023        first[length]=iNext;
     1024      }
     1025      if (iNext<numberRows_) {
     1026        previous[iNext]=iPrevious;
     1027      } else {
     1028        assert (iNext==length+2*large);
     1029      }
     1030      // put in
     1031      length=permute_[kRow];
     1032      smallest = CoinMin(smallest,length);
     1033      if (first[length]<0||first[length]>numberRows_) {
     1034        first[length]=kRow;
     1035        previous[kRow]=length+large;
     1036        next[kRow]=length+2*large;
     1037      } else {
     1038        int k=first[length];
     1039        assert (k<numberRows_);
     1040        first[length]=kRow;
     1041        previous[kRow]=length+large;
     1042        assert (previous[k]==length+large);
     1043        next[kRow]=k;
     1044        previous[k]=kRow;
     1045      }
     1046    }
     1047  }
     1048  // do rest
     1049  for (int iRow=0;iRow<numberRows_;iRow++) {
     1050    if (permuteInverse_[iRow]<0)
     1051      permuteInverse_[iRow]=done++;
     1052  }
     1053  printf("%d compressions, %d expansions\n",
     1054         numberCompressions,numberExpansions);
     1055  assert (done==numberRows_);
     1056  delete [] sort;
     1057  delete [] order;
     1058  delete [] which;
     1059  delete [] mark;
     1060  delete [] first;
     1061  delete [] next;
     1062  delete [] previous;
     1063  delete [] choleskyRow_;
     1064  choleskyRow_=NULL;
     1065  delete [] choleskyStart_;
     1066  choleskyStart_=NULL;
     1067#ifndef NDEBUG
     1068  for (int iRow=0;iRow<numberRows_;iRow++) {
     1069    permute_[iRow]=-1;
     1070  }
     1071#endif
     1072  for (int iRow=0;iRow<numberRows_;iRow++) {
     1073    permute_[permuteInverse_[iRow]]=iRow;
     1074  }
     1075#ifndef NDEBUG
     1076  for (int iRow=0;iRow<numberRows_;iRow++) {
     1077    assert(permute_[iRow]>=0&&permute_[iRow]<numberRows_);
     1078  }
     1079#endif
     1080  return returnCode;
     1081}
     1082#elif BASE_ORDER==2
     1083/*----------------------------------------------------------------------------*/
     1084/*      (C) Copyright IBM Corporation 1997, 2009.  All Rights Reserved.       */
     1085/*----------------------------------------------------------------------------*/
     1086
     1087/*  A compact no-frills Approximate Minimum Local Fill ordering code.
     1088
     1089    References:
     1090
     1091[1] Ordering Sparse Matrices Using Approximate Minimum Local Fill.
     1092    Edward Rothberg, SGI Manuscript, April 1996.
     1093[2] An Approximate Minimum Degree Ordering Algorithm.
     1094    T. Davis, P. Amestoy, and I. Duff, TR-94-039, CIS Department,
     1095    University of Florida, December 1994.
     1096*/
     1097
     1098#include <math.h>
     1099#include <stdlib.h>
     1100
     1101typedef int             WSI;
     1102
     1103#define NORDTHRESH      7
     1104#define MAXIW           2147000000
     1105
     1106//#define WSSMP_DBG
     1107#ifdef WSSMP_DBG
     1108static void chk (WSI ind, WSI i, WSI lim) {
     1109
     1110  if (i <= 0 || i > lim) {
     1111    printf("%d: chk: myamlf: WAR**: i, lim = %d, %d\n",ind,i,lim);
     1112  }
     1113}
     1114#endif
     1115
     1116static void
     1117myamlf(WSI n, WSI xadj[], WSI adjncy[], WSI dgree[], WSI varbl[],
     1118       WSI snxt[], WSI perm[], WSI invp[], WSI head[], WSI lsize[],
     1119       WSI flag[], WSI erscore[], WSI locaux, WSI adjln, WSI speed)
     1120{
     1121WSI l, i, j, k;
     1122double x, y;
     1123WSI maxmum, fltag, nodeg, scln, nm1, deg, tn,
     1124    locatns, ipp, jpp, nnode, numpiv, node,
     1125    nodeln, currloc, counter, numii, mindeg,
     1126    i0, i1, i2, i4, i5, i6, i7, i9,
     1127    j0, j1, j2, j3, j4, j5, j6, j7, j8, j9;
     1128/*                                             n cannot be less than NORDTHRESH
     1129 if (n < 3) {
     1130    if (n > 0) perm[0] = invp[0] = 1;
     1131    if (n > 1) perm[1] = invp[1] = 2;
     1132    return;
     1133 }
     1134*/
     1135#ifdef WSSMP_DBG
     1136 printf("Myamlf: n, locaux, adjln, speed = %d,%d,%d,%d\n",n,locaux,adjln,speed);
     1137 for (i = 0; i < n; i++) flag[i] = 0;
     1138 k = xadj[0];
     1139 for (i = 1; i <= n; i++) {
     1140   j = k + dgree[i-1];
     1141   if (j < k || k < 1) printf("ERR**: myamlf: i, j, k = %d,%d,%d\n",i,j,k);
     1142   for (l = k; l < j; l++) {
     1143      if (adjncy[l - 1] < 1 || adjncy[l - 1] > n || adjncy[l - 1] == i)
     1144        printf("ERR**: myamlf: i, l, adjj[l] = %d,%d,%d\n",i,l,adjncy[l - 1]);
     1145      if (flag[adjncy[l - 1] - 1] == i)
     1146        printf("ERR**: myamlf: duplicate entry: %d - %d\n",i,adjncy[l - 1]);
     1147      flag[adjncy[l - 1] - 1] = i;
     1148      nm1 = adjncy[l - 1] - 1;
     1149      for (fltag = xadj[nm1]; fltag < xadj[nm1] + dgree[nm1]; fltag++) {
     1150        if (adjncy[fltag - 1] == i) goto L100;
     1151      }
     1152      printf("ERR**: Unsymmetric %d %d\n",i,fltag);
     1153L100:;
     1154   }
     1155   k = j;
     1156   if (k > locaux) printf("ERR**: i, j, k, locaux = %d, %d, %d, %d\n",
     1157                                  i, j, k, locaux);
     1158 }
     1159 if (n*(n-1) < locaux-1) printf("ERR**: myamlf: too many nnz, n, ne = %d, %d\n",
     1160                                n,adjln);
     1161 for (i = 0; i < n; i++) flag[i] = 1;
     1162 if (n > 10000) printf ("Finished error checking in AMF\n");
     1163 for (i = locaux; i <= adjln; i++) adjncy[i - 1] = -22;
     1164#endif
     1165
     1166 maxmum = 0;
     1167 mindeg = 1;
     1168 fltag = 2;
     1169 locatns = locaux - 1;
     1170 nm1 = n - 1;
     1171 counter = 1;
     1172 for (l = 0; l < n; l++) {
     1173     j = erscore[l];
     1174#ifdef WSSMP_DBG
     1175chk(1,j,n);
     1176#endif
     1177     if (j > 0) {
     1178         nnode = head[j - 1];
     1179         if (nnode) perm[nnode - 1] = l + 1;
     1180         snxt[l] = nnode;
     1181         head[j - 1] = l + 1;
     1182     } else {
     1183         invp[l] = -(counter++);
     1184         flag[l] = xadj[l] = 0;
     1185     }
     1186 }
     1187 while (counter <= n) {
     1188   for (deg = mindeg; head[deg - 1] < 1; deg++) {};
     1189     nodeg = 0;
     1190#ifdef WSSMP_DBG
     1191chk(2,deg,n-1);
     1192#endif
     1193     node = head[deg - 1];
     1194#ifdef WSSMP_DBG
     1195chk(3,node,n);
     1196#endif
     1197     mindeg = deg;
     1198     nnode = snxt[node - 1];
     1199     if (nnode) perm[nnode - 1] = 0;
     1200     head[deg - 1] = nnode;
     1201     nodeln = invp[node - 1];
     1202     numpiv = varbl[node - 1];
     1203     invp[node - 1] = -counter;
     1204     counter += numpiv;
     1205     varbl[node - 1] = -numpiv;
     1206     if (nodeln) {
     1207         j4 = locaux;
     1208         i5 = lsize[node - 1] - nodeln;
     1209         i2 = nodeln + 1;
     1210         l = xadj[node - 1];
     1211         for (i6 = 1; i6 <= i2; ++i6) {
     1212             if (i6 == i2) {
     1213                 tn = node, i0 = l, scln = i5;
     1214             } else {
     1215#ifdef WSSMP_DBG
     1216chk(4,l,adjln);
     1217#endif
     1218                 tn = adjncy[l-1];
     1219                 l++;
     1220#ifdef WSSMP_DBG
     1221chk(5,tn,n);
     1222#endif
     1223                 i0 = xadj[tn - 1];
     1224                 scln = lsize[tn - 1];
     1225             }
     1226             for (i7 = 1; i7 <= scln; ++i7) {
     1227#ifdef WSSMP_DBG
     1228chk(6,i0,adjln);
     1229#endif
     1230               i = adjncy[i0 - 1];
     1231               i0++;
     1232#ifdef WSSMP_DBG
     1233chk(7,i,n);
     1234#endif
     1235               numii = varbl[i - 1];
     1236               if (numii > 0) {
     1237                 if (locaux > adjln) {
     1238                   lsize[node - 1] -= i6;
     1239                   xadj[node - 1] = l;
     1240                   if (!lsize[node - 1]) xadj[node - 1] = 0;
     1241                   xadj[tn - 1] = i0;
     1242                   lsize[tn - 1] = scln - i7;
     1243                   if (!lsize[tn - 1]) xadj[tn - 1] = 0;
     1244                   for (j = 1; j <= n; j++) {
     1245                     i4 = xadj[j - 1];
     1246                     if (i4 > 0) {
     1247                       xadj[j - 1] = adjncy[i4 - 1];
     1248#ifdef WSSMP_DBG
     1249chk(8,i4,adjln);
     1250#endif
     1251                       adjncy[i4 - 1] = -j;
     1252                     }
     1253                   }
     1254                   i9 = j4 - 1;
     1255                   j6 = j7 = 1;
     1256                   while (j6 <= i9) {
     1257#ifdef WSSMP_DBG
     1258chk(9,j6,adjln);
     1259#endif
     1260                     j = -adjncy[j6 - 1];
     1261                     j6++;
     1262                     if (j > 0) {
     1263#ifdef WSSMP_DBG
     1264chk(10,j7,adjln);
     1265chk(11,j,n);
     1266#endif
     1267                       adjncy[j7 - 1] = xadj[j - 1];
     1268                       xadj[j - 1] = j7++;
     1269                       j8 = lsize[j - 1] - 1 + j7;
     1270                       for (; j7 < j8; j7++,j6++) adjncy[j7-1] = adjncy[j6-1];
     1271                     }
     1272                   }
     1273                   for (j0=j7; j4 < locaux; j4++,j7++) {
     1274#ifdef WSSMP_DBG
     1275chk(12,j4,adjln);
     1276#endif
     1277                     adjncy[j7 - 1] = adjncy[j4 - 1];
     1278                   }
     1279                   j4 = j0;
     1280                   locaux = j7;
     1281                   i0 = xadj[tn - 1];
     1282                   l = xadj[node - 1];
     1283                 }
     1284                 adjncy[locaux-1] = i;
     1285                 locaux++;
     1286                 varbl[i - 1] = -numii;
     1287                 nodeg += numii;
     1288                 ipp = perm[i - 1];
     1289                 nnode = snxt[i - 1];
     1290#ifdef WSSMP_DBG
     1291if (ipp) chk(13,ipp,n);
     1292if (nnode) chk(14,nnode,n);
     1293#endif
     1294                 if (ipp) snxt[ipp - 1] = nnode;
     1295                 else head[erscore[i - 1] - 1] = nnode;
     1296                 if (nnode) perm[nnode - 1] = ipp;
     1297               }
     1298             }
     1299             if (tn != node) {
     1300                 flag[tn - 1] = 0, xadj[tn - 1] = -node;
     1301             }
     1302         }
     1303         currloc = (j5 = locaux) - j4;
     1304         locatns += currloc;
     1305     } else {
     1306         i1 = (j4 = xadj[node - 1]) + lsize[node - 1];
     1307         for (j = j5 = j4; j < i1; j++) {
     1308#ifdef WSSMP_DBG
     1309chk(15,j,adjln);
     1310#endif
     1311             i = adjncy[j - 1];
     1312#ifdef WSSMP_DBG
     1313chk(16,i,n);
     1314#endif
     1315             numii = varbl[i - 1];
     1316             if (numii > 0) {
     1317                 nodeg += numii;
     1318                 varbl[i - 1] = -numii;
     1319#ifdef WSSMP_DBG
     1320chk(17,j5,adjln);
     1321#endif
     1322                 adjncy[j5-1] = i;
     1323                 ipp = perm[i - 1];
     1324                 nnode = snxt[i - 1];
     1325                 j5++;
     1326#ifdef WSSMP_DBG
     1327if (ipp) chk(18,ipp,n);
     1328if (nnode) chk(19,nnode,n);
     1329#endif
     1330                 if (ipp) snxt[ipp - 1] = nnode;
     1331                 else head[erscore[i - 1] - 1] = nnode;
     1332                 if (nnode) perm[nnode - 1] = ipp;
     1333             }
     1334         }
     1335         currloc = 0;
     1336     }
     1337#ifdef WSSMP_DBG
     1338chk(20,node,n);
     1339#endif
     1340     lsize[node - 1] = j5 - (xadj[node - 1] = j4);
     1341     dgree[node - 1] = nodeg;
     1342     if (fltag+n < 0 || fltag+n > MAXIW) {
     1343       for (i = 1; i <= n; i++) if (flag[i - 1]) flag[i - 1] = 1;
     1344       fltag = 2;
     1345     }
     1346     for (j3 = j4; j3 < j5; j3++) {
     1347#ifdef WSSMP_DBG
     1348chk(21,j3,adjln);
     1349#endif
     1350       i = adjncy[j3 - 1];
     1351#ifdef WSSMP_DBG
     1352chk(22,i,n);
     1353#endif
     1354       j = invp[i - 1];
     1355       if (j > 0) {
     1356         i4 = fltag - (numii = -varbl[i - 1]);
     1357         i2 = xadj[i - 1] + j;
     1358         for (l = xadj[i - 1]; l < i2; l++) {
     1359#ifdef WSSMP_DBG
     1360chk(23,l,adjln);
     1361#endif
     1362           tn = adjncy[l - 1];
     1363#ifdef WSSMP_DBG
     1364chk(24,tn,n);
     1365#endif
     1366           j9 = flag[tn - 1];
     1367           if (j9 >= fltag) j9 -= numii; else if (j9) j9 = dgree[tn - 1] + i4;
     1368           flag[tn - 1] = j9;
     1369         }
     1370       }
     1371     }
     1372     for (j3 = j4; j3 < j5; j3++) {
     1373#ifdef WSSMP_DBG
     1374chk(25,j3,adjln);
     1375#endif
     1376         i = adjncy[j3 - 1];
     1377         i5 = deg = 0;
     1378#ifdef WSSMP_DBG
     1379chk(26,i,n);
     1380#endif
     1381         j1 = (i4 = xadj[i - 1]) + invp[i - 1];
     1382         for (l = j0 = i4; l < j1; l++) {
     1383#ifdef WSSMP_DBG
     1384chk(27,l,adjln);
     1385#endif
     1386             tn = adjncy[l - 1];
     1387#ifdef WSSMP_DBG
     1388chk(70,tn,n);
     1389#endif
     1390             j8 = flag[tn - 1];
     1391             if (j8) {
     1392                 deg += (j8 - fltag);
     1393#ifdef WSSMP_DBG
     1394chk(28,i4,adjln);
     1395#endif
     1396                 adjncy[i4-1] = tn;
     1397                 i5 += tn;
     1398                 i4++;
     1399                 while (i5 >= nm1) i5 -= nm1;
     1400             }
     1401         }
     1402         invp[i - 1] = (j2 = i4) - j0 + 1;
     1403         i2 = j0 + lsize[i - 1];
     1404         for (l = j1; l < i2; l++) {
     1405#ifdef WSSMP_DBG
     1406chk(29,l,adjln);
     1407#endif
     1408             j = adjncy[l - 1];
     1409#ifdef WSSMP_DBG
     1410chk(30,j,n);
     1411#endif
     1412             numii = varbl[j - 1];
     1413             if (numii > 0) {
     1414                 deg += numii;
     1415#ifdef WSSMP_DBG
     1416chk(31,i4,adjln);
     1417#endif
     1418                 adjncy[i4-1] = j;
     1419                 i5 += j;
     1420                 i4++;
     1421                 while (i5 >= nm1) i5 -= nm1;
     1422             }
     1423         }
     1424         if (invp[i - 1] == 1 && j2 == i4) {
     1425             numii = -varbl[i - 1];
     1426             xadj[i - 1] = -node;
     1427             nodeg -= numii;
     1428             counter += numii;
     1429             numpiv += numii;
     1430             invp[i - 1] = varbl[i - 1] = 0;
     1431         } else {
     1432             if (dgree[i - 1] > deg) dgree[i - 1] = deg;
     1433#ifdef WSSMP_DBG
     1434chk(32,j2,adjln);
     1435chk(33,j0,adjln);
     1436#endif
     1437             adjncy[i4 - 1] = adjncy[j2 - 1];
     1438             adjncy[j2 - 1] = adjncy[j0 - 1];
     1439             adjncy[j0 - 1] = node;
     1440             lsize[i - 1] = i4 - j0 + 1;
     1441             i5++;
     1442#ifdef WSSMP_DBG
     1443chk(35,i5,n);
     1444#endif
     1445             j = head[i5 - 1];
     1446             if (j > 0) {
     1447                 snxt[i - 1] = perm[j - 1];
     1448                 perm[j - 1] = i;
     1449             } else {
     1450                 snxt[i - 1] = -j;
     1451                 head[i5 - 1] = -i;
     1452             }
     1453             perm[i - 1] = i5;
     1454         }
     1455     }
     1456#ifdef WSSMP_DBG
     1457chk(36,node,n);
     1458#endif
     1459     dgree[node - 1] = nodeg;
     1460     if (maxmum < nodeg) maxmum = nodeg;
     1461     fltag += maxmum;
     1462#ifdef WSSMP_DBG
     1463     if (fltag+n+1 < 0) printf("ERR**: fltag = %d (A)\n",fltag);
     1464#endif
     1465     for (j3 = j4; j3 < j5; ++j3) {
     1466#ifdef WSSMP_DBG
     1467chk(37,j3,adjln);
     1468#endif
     1469       i = adjncy[j3 - 1];
     1470#ifdef WSSMP_DBG
     1471chk(38,i,n);
     1472#endif
     1473       if (varbl[i - 1] < 0) {
     1474         i5 = perm[i - 1];
     1475#ifdef WSSMP_DBG
     1476chk(39,i5,n);
     1477#endif
     1478         j = head[i5 - 1];
     1479         if (j) {
     1480           if (j < 0) {
     1481                head[i5 - 1] = 0, i = -j;
     1482           } else {
     1483#ifdef WSSMP_DBG
     1484chk(40,j,n);
     1485#endif
     1486                i = perm[j - 1];
     1487                perm[j - 1] = 0;
     1488           }
     1489           while (i) {
     1490#ifdef WSSMP_DBG
     1491chk(41,i,n);
     1492#endif
     1493             if (!snxt[i - 1]) {
     1494                 i = 0;
     1495             } else {
     1496               k = invp[i - 1];
     1497               i2 = xadj[i - 1] + (scln = lsize[i - 1]);
     1498               for (l = xadj[i - 1] + 1; l < i2; l++) {
     1499#ifdef WSSMP_DBG
     1500chk(42,l,adjln);
     1501chk(43,adjncy[l - 1],n);
     1502#endif
     1503                 flag[adjncy[l - 1] - 1] = fltag;
     1504               }
     1505               jpp = i;
     1506               j = snxt[i - 1];
     1507               while (j) {
     1508#ifdef WSSMP_DBG
     1509chk(44,j,n);
     1510#endif
     1511                 if (lsize[j - 1] == scln && invp[j - 1] == k) {
     1512                   i2 = xadj[j - 1] + scln;
     1513                   j8 = 1;
     1514                   for (l = xadj[j - 1] + 1; l < i2; l++) {
     1515#ifdef WSSMP_DBG
     1516chk(45,l,adjln);
     1517chk(46,adjncy[l - 1],n);
     1518#endif
     1519                       if (flag[adjncy[l - 1] - 1] != fltag) {
     1520                         j8 = 0;
     1521                         break;
     1522                       }
     1523                   }
     1524                   if (j8) {
     1525                     xadj[j - 1] = -i;
     1526                     varbl[i - 1] += varbl[j - 1];
     1527                     varbl[j - 1] = invp[j - 1] = 0;
     1528#ifdef WSSMP_DBG
     1529chk(47,j,n);
     1530chk(48,jpp,n);
     1531#endif
     1532                     j = snxt[j - 1];
     1533                     snxt[jpp - 1] = j;
     1534                   } else {
     1535                     jpp = j;
     1536#ifdef WSSMP_DBG
     1537chk(49,j,n);
     1538#endif
     1539                     j = snxt[j - 1];
     1540                   }
     1541                 } else {
     1542                   jpp = j;
     1543#ifdef WSSMP_DBG
     1544chk(50,j,n);
     1545#endif
     1546                   j = snxt[j - 1];
     1547                 }
     1548               }
     1549               fltag++;
     1550#ifdef WSSMP_DBG
     1551               if (fltag+n+1 < 0) printf("ERR**: fltag = %d (B)\n",fltag);
     1552#endif
     1553#ifdef WSSMP_DBG
     1554chk(51,i,n);
     1555#endif
     1556               i = snxt[i - 1];
     1557             }
     1558           }
     1559         }
     1560       }
     1561     }
     1562     j8 = nm1 - counter;
     1563     switch (speed) {
     1564case 1:
     1565       for (j = j3 = j4; j3 < j5; j3++) {
     1566#ifdef WSSMP_DBG
     1567chk(52,j3,adjln);
     1568#endif
     1569         i = adjncy[j3 - 1];
     1570#ifdef WSSMP_DBG
     1571chk(53,i,n);
     1572#endif
     1573         numii = varbl[i - 1];
     1574         if (numii < 0) {
     1575           k = 0;
     1576           i4 = (l = xadj[i - 1]) + invp[i - 1];
     1577           for (; l < i4; l++) {
     1578#ifdef WSSMP_DBG
     1579chk(54,l,adjln);
     1580chk(55,adjncy[l - 1],n);
     1581#endif
     1582              i5 = dgree[adjncy[l - 1] - 1];
     1583              if (k < i5) k = i5;
     1584           }
     1585           x = static_cast<double> (k - 1);
     1586           varbl[i - 1] = abs(numii);
     1587           j9 = dgree[i - 1] + nodeg;
     1588           deg = (j8 > j9 ? j9 : j8) + numii;
     1589           if (deg < 1) deg = 1;
     1590           y = static_cast<double> (deg);
     1591           dgree[i - 1] = deg;
     1592/*
     1593           printf("%i %f should >= %i %f\n",deg,y,k-1,x);
     1594           if (y < x) printf("** \n"); else printf("\n");
     1595*/
     1596           y = y*y - y;
     1597           x = y - (x*x - x);
     1598           if (x < 1.1) x = 1.1;
     1599           deg = static_cast<WSI>(sqrt(x));
     1600/*         if (deg < 1) deg = 1; */
     1601           erscore[i - 1] = deg;
     1602#ifdef WSSMP_DBG
     1603chk(56,deg,n-1);
     1604#endif
     1605           nnode = head[deg - 1];
     1606           if (nnode) perm[nnode - 1] = i;
     1607           snxt[i - 1] = nnode;
     1608           perm[i - 1] = 0;
     1609#ifdef WSSMP_DBG
     1610chk(57,j,adjln);
     1611#endif
     1612           head[deg - 1] = adjncy[j-1] = i;
     1613           j++;
     1614           if (deg < mindeg) mindeg = deg;
     1615         }
     1616       }
     1617       break;
     1618case 2:
     1619       for (j = j3 = j4; j3 < j5; j3++) {
     1620#ifdef WSSMP_DBG
     1621chk(58,j3,adjln);
     1622#endif
     1623         i = adjncy[j3 - 1];
     1624#ifdef WSSMP_DBG
     1625chk(59,i,n);
     1626#endif
     1627         numii = varbl[i - 1];
     1628         if (numii < 0) {
     1629           x = static_cast<double> (dgree[adjncy[xadj[i - 1] - 1] - 1] - 1);
     1630           varbl[i - 1] = abs(numii);
     1631           j9 = dgree[i - 1] + nodeg;
     1632           deg = (j8 > j9 ? j9 : j8) + numii;
     1633           if (deg < 1) deg = 1;
     1634           y = static_cast<double> (deg);
     1635           dgree[i - 1] = deg;
     1636/*
     1637           printf("%i %f should >= %i %f",deg,y,dgree[adjncy[xadj[i - 1] - 1] - 1]-1,x);
     1638           if (y < x) printf("** \n"); else printf("\n");
     1639*/
     1640           y = y*y - y;
     1641           x = y - (x*x - x);
     1642           if (x < 1.1) x = 1.1;
     1643           deg = static_cast<WSI>(sqrt(x));
     1644/*         if (deg < 1) deg = 1; */
     1645           erscore[i - 1] = deg;
     1646#ifdef WSSMP_DBG
     1647chk(60,deg,n-1);
     1648#endif
     1649           nnode = head[deg - 1];
     1650           if (nnode) perm[nnode - 1] = i;
     1651           snxt[i - 1] = nnode;
     1652           perm[i - 1] = 0;
     1653#ifdef WSSMP_DBG
     1654chk(61,j,adjln);
     1655#endif
     1656           head[deg - 1] = adjncy[j-1] = i;
     1657           j++;
     1658           if (deg < mindeg) mindeg = deg;
     1659         }
     1660       }
     1661       break;
     1662default:
     1663       for (j = j3 = j4; j3 < j5; j3++) {
     1664#ifdef WSSMP_DBG
     1665chk(62,j3,adjln);
     1666#endif
     1667         i = adjncy[j3 - 1];
     1668#ifdef WSSMP_DBG
     1669chk(63,i,n);
     1670#endif
     1671         numii = varbl[i - 1];
     1672         if (numii < 0) {
     1673           varbl[i - 1] = abs(numii);
     1674           j9 = dgree[i - 1] + nodeg;
     1675           deg = (j8 > j9 ? j9 : j8) + numii;
     1676           if (deg < 1) deg = 1;
     1677           dgree[i - 1] = deg;
     1678#ifdef WSSMP_DBG
     1679chk(64,deg,n-1);
     1680#endif
     1681           nnode = head[deg - 1];
     1682           if (nnode) perm[nnode - 1] = i;
     1683           snxt[i - 1] = nnode;
     1684           perm[i - 1] = 0;
     1685#ifdef WSSMP_DBG
     1686chk(65,j,adjln);
     1687#endif
     1688           head[deg - 1] = adjncy[j-1] = i;
     1689           j++;
     1690           if (deg < mindeg) mindeg = deg;
     1691         }
     1692       }
     1693       break;
     1694     }
     1695     if (currloc) {
     1696#ifdef WSSMP_DBG
     1697chk(66,node,n);
     1698#endif
     1699       locatns += (lsize[node - 1] - currloc), locaux = j;
     1700     }
     1701     varbl[node - 1] = numpiv + nodeg;
     1702     lsize[node - 1] = j - j4;
     1703     if (!lsize[node - 1]) flag[node - 1] = xadj[node - 1] = 0;
     1704 }
     1705 for (l = 1; l <= n; l++) {
     1706   if (!invp[l - 1]) {
     1707     for (i = -xadj[l - 1]; invp[i - 1] >= 0; i = -xadj[i - 1]){};
     1708       tn = i;
     1709#ifdef WSSMP_DBG
     1710chk(67,tn,n);
     1711#endif
     1712       k = -invp[tn - 1];
     1713       i = l;
     1714#ifdef WSSMP_DBG
     1715chk(68,i,n);
     1716#endif
     1717       while (invp[i - 1] >= 0) {
     1718         nnode = -xadj[i - 1];
     1719         xadj[i - 1] = -tn;
     1720         if (!invp[i - 1]) invp[i - 1] = k++;
     1721         i = nnode;
     1722       }
     1723       invp[tn - 1] = -k;
     1724   }
     1725 }
     1726 for (l = 0; l < n; l++) {
     1727   i = abs(invp[l]);
     1728#ifdef WSSMP_DBG
     1729chk(69,i,n);
     1730#endif
     1731   invp[l] = i;
     1732   perm[i - 1] = l + 1;
     1733 }
     1734 return;
     1735}
     1736// This code is not needed, but left in in case needed sometime
     1737#if 0
     1738/*C--------------------------------------------------------------------------*/
     1739void amlfdr(WSI *n, WSI xadj[], WSI adjncy[], WSI dgree[], WSI *adjln,
     1740            WSI *locaux, WSI varbl[], WSI snxt[], WSI perm[],
     1741            WSI head[], WSI invp[], WSI lsize[], WSI flag[], WSI *ispeed)
     1742{
     1743WSI nn,nlocaux,nadjln,speed,i,j,mx,mxj,*erscore;
     1744
     1745#ifdef WSSMP_DBG
     1746  printf("Calling amlfdr for n, speed = %d, %d\n",*n,*ispeed);
     1747#endif
     1748
     1749  if ((nn = *n) == 0) return;
     1750
     1751#ifdef WSSMP_DBG
     1752  if (nn == 31) {
     1753    printf("n = %d; adjln = %d; locaux = %d; ispeed = %d\n",
     1754            *n,*adjln,*locaux,*ispeed);
     1755  }
     1756#endif
     1757
     1758  if (nn < NORDTHRESH) {
     1759    for (i = 0; i < nn; i++) lsize[i] = i;
     1760    for (i = nn; i > 0; i--) {
     1761      mx = dgree[0];
     1762      mxj = 0;
     1763      for (j = 1; j < i; j++)
     1764        if (dgree[j] > mx) {
     1765          mx = dgree[j];
     1766          mxj = j;
     1767        }
     1768      invp[lsize[mxj]] = i;
     1769      dgree[mxj] = dgree[i-1];
     1770      lsize[mxj] = lsize[i-1];
     1771    }
     1772    return;
     1773  }
     1774  speed = *ispeed;
     1775  if (speed < 3) {
     1776/* 
     1777    erscore = (WSI *)malloc(nn * sizeof(WSI));
     1778    if (erscore == NULL) speed = 3;
     1779*/
     1780    wscmal (&nn, &i, &erscore);
     1781    if (i != 0) speed = 3;
     1782  }
     1783  if (speed > 2) erscore = dgree;
     1784  if (speed < 3) {
     1785    for (i = 0; i < nn; i++) {
     1786      perm[i] = 0;
     1787      invp[i] = 0;
     1788      head[i] = 0;
     1789      flag[i] = 1;
     1790      varbl[i] = 1;
     1791      lsize[i] = dgree[i];
     1792      erscore[i] = dgree[i];
     1793    }
     1794  } else {
     1795    for (i = 0; i < nn; i++) {
     1796      perm[i] = 0;
     1797      invp[i] = 0;
     1798      head[i] = 0;
     1799      flag[i] = 1;
     1800      varbl[i] = 1;
     1801      lsize[i] = dgree[i];
     1802    }
     1803  }
     1804  nlocaux = *locaux;
     1805  nadjln = *adjln;
     1806
     1807  myamlf(nn, xadj, adjncy, dgree, varbl, snxt, perm, invp,
     1808         head, lsize, flag, erscore, nlocaux, nadjln, speed);
     1809/*
     1810  if (speed < 3) free(erscore);
     1811*/
     1812  if (speed < 3) wscfr(&erscore);
     1813  return;
     1814}
     1815#endif // end of taking out amlfdr
     1816/*C--------------------------------------------------------------------------*/
     1817#endif
     1818// Orders rows
     1819int
     1820ClpCholeskyBase::orderAMD()
     1821{
     1822  permuteInverse_ = new int [numberRows_];
     1823  permute_ = new int[numberRows_];
     1824  // Do ordering
     1825  int returnCode=0;
     1826  // get full matrix
     1827  int space = 2*sizeFactor_+10000+4*numberRows_;
     1828  int * temp = new int [space];
     1829  CoinBigIndex * count = new CoinBigIndex [numberRows_];
     1830  CoinBigIndex * tempStart = new CoinBigIndex [numberRows_+1];
     1831  memset(count,0,numberRows_*sizeof(int));
     1832  for (int iRow=0;iRow<numberRows_;iRow++) {
     1833    count[iRow]+= choleskyStart_[iRow+1]-choleskyStart_[iRow]-1;
     1834    for (CoinBigIndex j=choleskyStart_[iRow]+1;j<choleskyStart_[iRow+1];j++) {
     1835      int jRow=choleskyRow_[j];
     1836      count[jRow]++;
     1837    }
     1838  }
     1839#define OFFSET 1
     1840  CoinBigIndex sizeFactor =0;
     1841  for (int iRow=0;iRow<numberRows_;iRow++) {
     1842    int length = count[iRow];
     1843    permute_[iRow]=length;
     1844    // add 1 to starts
     1845    tempStart[iRow]=sizeFactor+OFFSET;
     1846    count[iRow]=sizeFactor;
     1847    sizeFactor += length;
     1848  }
     1849  tempStart[numberRows_]=sizeFactor+OFFSET;
     1850  // add 1 to rows
     1851  for (int iRow=0;iRow<numberRows_;iRow++) {
     1852    assert (choleskyRow_[choleskyStart_[iRow]]==iRow);
     1853    for (CoinBigIndex j=choleskyStart_[iRow]+1;j<choleskyStart_[iRow+1];j++) {
     1854      int jRow=choleskyRow_[j];
     1855      int put = count[iRow];
     1856      temp[put]=jRow+OFFSET;
     1857      count[iRow]++;
     1858      put = count[jRow];
     1859      temp[put]=iRow+OFFSET;
     1860      count[jRow]++;
     1861    }
     1862  }
     1863  for (int iRow=1;iRow<numberRows_;iRow++)
     1864    assert (count[iRow-1]==tempStart[iRow]-OFFSET);
     1865  delete [] choleskyRow_;
     1866  choleskyRow_=temp;
     1867  delete [] choleskyStart_;
     1868  choleskyStart_=tempStart;
     1869  int locaux=sizeFactor+OFFSET;
     1870  delete [] count;
     1871  int speed=integerParameters_[0];
     1872  if (speed<1||speed>2)
     1873    speed=3;
     1874  int * use = new int [((speed<3) ? 7 : 6)*numberRows_];
     1875  int * dgree = use;
     1876  int * varbl = dgree+numberRows_;
     1877  int * snxt = varbl+numberRows_;
     1878  int * head = snxt+numberRows_;
     1879  int * lsize = head+numberRows_;
     1880  int * flag = lsize+numberRows_;
     1881  int * erscore;
     1882  for (int i=0;i<numberRows_;i++) {
     1883    dgree[i]=choleskyStart_[i+1]-choleskyStart_[i];
     1884    head[i]=dgree[i];
     1885    snxt[i]=0;
     1886    permute_[i]=0;
     1887    permuteInverse_[i]=0;
     1888    head[i]=0;
     1889    flag[i]=1;
     1890    varbl[i]=1;
     1891    lsize[i]=dgree[i];
     1892  }
     1893  if (speed<3) {
     1894    erscore = flag+numberRows_;
     1895    for (int i=0;i<numberRows_;i++)
     1896      erscore[i]=dgree[i];
     1897  } else {
     1898    erscore = dgree;
     1899  }
     1900  myamlf(numberRows_, choleskyStart_, choleskyRow_,
     1901         dgree, varbl, snxt, permute_, permuteInverse_,
     1902         head, lsize, flag, erscore, locaux, space, speed);
     1903  for (int iRow=0;iRow<numberRows_;iRow++) {
     1904    permute_[iRow]--;
     1905  }
     1906  for (int iRow=0;iRow<numberRows_;iRow++) {
     1907    permuteInverse_[permute_[iRow]]=iRow;
     1908  }
     1909  for (int iRow=0;iRow<numberRows_;iRow++) {
     1910    assert (permuteInverse_[iRow]>=0&&permuteInverse_[iRow]<numberRows_);
     1911  }
     1912  delete [] use;
     1913  delete [] choleskyRow_;
     1914  choleskyRow_=NULL;
     1915  delete [] choleskyStart_;
     1916  choleskyStart_=NULL;
     1917  return returnCode;
    7591918}
    7601919/* Does Symbolic factorization given permutation.
     
    12372396#else
    12382397    // actually long double
    1239     workDouble_ = (double *) (new longWork[numberRows_]);
     2398    workDouble_ = reinterpret_cast<double *> (new CoinWorkDouble[numberRows_]);
    12402399#endif
    12412400    diagonal_ = new longDouble[numberRows_];
     
    13382497      int k=iRow;
    13392498      int linked = link_[iRow];
     2499#ifndef NDEBUG
     2500      int count=0;
     2501#endif
    13402502      while (linked<=kRow) {
    13412503        k=linked;
    13422504        linked = link_[k];
     2505#ifndef NDEBUG
     2506        count++;
     2507        assert (count<numberRows_);
     2508#endif
    13432509      }
    13442510      nz++;
     
    15212687/* Factorize - filling in rowsDropped and returning number dropped */
    15222688int
    1523 ClpCholeskyBase::factorize(const double * diagonal, int * rowsDropped)
     2689ClpCholeskyBase::factorize(const CoinWorkDouble * diagonal, int * rowsDropped)
    15242690{
    15252691  const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
     
    15332699  int numberColumns=model_->clpMatrix()->getNumCols();
    15342700  //perturbation
    1535   longDouble perturbation=model_->diagonalPerturbation()*model_->diagonalNorm();
     2701  CoinWorkDouble perturbation=model_->diagonalPerturbation()*model_->diagonalNorm();
    15362702  //perturbation=perturbation*perturbation*100000000.0;
    15372703  if (perturbation>1.0) {
     
    15402706      std::cout <<"large perturbation "<<perturbation<<std::endl;
    15412707#endif
    1542     perturbation=sqrt(perturbation);;
     2708    perturbation=CoinSqrt(perturbation);
    15432709    perturbation=1.0;
    15442710  }
     
    15482714  CoinZeroN(work,numberRows_);
    15492715  int newDropped=0;
    1550   double largest=1.0;
    1551   double smallest=COIN_DBL_MAX;
     2716  CoinWorkDouble largest=1.0;
     2717  CoinWorkDouble smallest=COIN_DBL_MAX;
    15522718  int numberDense=0;
    15532719  if (!doKKT_) {
    1554     const double * diagonalSlack = diagonal + numberColumns;
     2720    const CoinWorkDouble * diagonalSlack = diagonal + numberColumns;
    15552721    if (dense_)
    15562722      numberDense=dense_->numberRows();
     
    15772743      }
    15782744    }
    1579     longDouble delta2 = model_->delta(); // add delta*delta to diagonal
     2745    CoinWorkDouble delta2 = model_->delta(); // add delta*delta to diagonal
    15802746    delta2 *= delta2;
    15812747    // largest in initial matrix
    1582     double largest2=1.0e-20;
     2748    CoinWorkDouble largest2=1.0e-20;
    15832749    for (iRow=0;iRow<numberRows_;iRow++) {
    15842750      longDouble * put = sparseFactor_+choleskyStart_[iRow];
     
    15972763            CoinBigIndex start=columnStart[iColumn];
    15982764            CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
    1599             longDouble multiplier = diagonal[iColumn]*elementByRow[k];
     2765            CoinWorkDouble multiplier = diagonal[iColumn]*elementByRow[k];
    16002766            for (CoinBigIndex j=start;j<end;j++) {
    16012767              int jRow=row[j];
    16022768              int jNewRow = permuteInverse_[jRow];
    16032769              if (jNewRow>=iRow&&!rowsDropped_[jRow]) {
    1604                 longDouble value=element[j]*multiplier;
     2770                CoinWorkDouble value=element[j]*multiplier;
    16052771                work[jNewRow] += value;
    16062772              }
     
    16092775        }
    16102776        diagonal_[iRow]=work[iRow];
    1611         largest2 = CoinMax(largest2,fabs(work[iRow]));
     2777        largest2 = CoinMax(largest2,CoinAbs(work[iRow]));
    16122778        work[iRow]=0.0;
    16132779        int j;
     
    16152781          int jRow = which[j];
    16162782          put[j]=work[jRow];
    1617           largest2 = CoinMax(largest2,fabs(work[jRow]));
     2783          largest2 = CoinMax(largest2,CoinAbs(work[jRow]));
    16182784          work[jRow]=0.0;
    16192785        }
     
    16292795    //check sizes
    16302796    largest2*=1.0e-20;
    1631     largest = CoinMin(largest2,1.0e-11);
     2797    largest = CoinMin(largest2,CHOL_SMALL_VALUE);
    16322798    int numberDroppedBefore=0;
    16332799    for (iRow=0;iRow<numberRows_;iRow++) {
     
    16362802      rowsDropped[iRow]=dropped;
    16372803      if (!dropped) {
    1638         longDouble diagonal = diagonal_[iRow];
     2804        CoinWorkDouble diagonal = diagonal_[iRow];
    16392805        if (diagonal>largest2) {
    16402806          diagonal_[iRow]=diagonal+perturbation;
     
    16632829      for ( i=0;i<numberRows_;i++) {
    16642830        assert (diagonal_[i]>=0.0);
    1665         diagonal_[i]=sqrt(diagonal_[i]);
     2831        diagonal_[i]=CoinSqrt(diagonal_[i]);
    16662832      }
    16672833      // Update dense columns (just L)
     
    16732839            a[j]=0.0;
    16742840        }
    1675         solveLong(a,1);
     2841        for (i=0;i<numberRows_;i++) {
     2842          int iRow = permute_[i];
     2843          workDouble_[i] = a[iRow];
     2844        }
     2845        for (i=0;i<numberRows_;i++) {
     2846          CoinWorkDouble value=workDouble_[i];
     2847          CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
     2848          CoinBigIndex j;
     2849          for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
     2850            int iRow = choleskyRow_[j+offset];
     2851            workDouble_[iRow] -= sparseFactor_[j]*value;
     2852          }
     2853        }
     2854        for (i=0;i<numberRows_;i++) {
     2855          int iRow = permute_[i];
     2856          a[iRow]=workDouble_[i]*diagonal_[i];
     2857        }
    16762858      }
    16772859      dense_->resetRowsDropped();
     
    16822864        const longDouble * a = denseColumn_+i*numberRows_;
    16832865        // do diagonal
    1684         longDouble value = denseDiagonal[i];
     2866        CoinWorkDouble value = denseDiagonal[i];
    16852867        const longDouble * b = denseColumn_+i*numberRows_;
    16862868        for (int k=0;k<numberRows_;k++)
     
    16882870        denseDiagonal[i]=value;
    16892871        for (int j=i+1;j<numberDense;j++) {
    1690           longDouble value = 0.0;
     2872          CoinWorkDouble value = 0.0;
    16912873          const longDouble * b = denseColumn_+j*numberRows_;
    16922874          for (int k=0;k<numberRows_;k++)
     
    17782960    }
    17792961    if (permuted) {
    1780       longDouble delta2 = model_->delta(); // add delta*delta to bottom
     2962      CoinWorkDouble delta2 = model_->delta(); // add delta*delta to bottom
    17812963      delta2 *= delta2;
    17822964      // Need to permute - ugly
     
    17882970          if (iOriginalRow<numberColumns) {
    17892971            iColumn=iOriginalRow;
    1790             longDouble value = diagonal[iColumn];
    1791             if (fabs(value)>1.0e-100) {
     2972            CoinWorkDouble value = diagonal[iColumn];
     2973            if (CoinAbs(value)>1.0e-100) {
    17922974              value = 1.0/value;
    1793               largest = CoinMax(largest,fabs(value));
     2975              largest = CoinMax(largest,CoinAbs(value));
    17942976              diagonal_[iRow] = -value;
    17952977              CoinBigIndex start=columnStart[iColumn];
     
    18002982                if (kRow>iRow) {
    18012983                  work[kRow]=element[j];
    1802                   largest = CoinMax(largest,fabs(element[j]));
     2984                  largest = CoinMax(largest,CoinAbs(element[j]));
    18032985                }
    18042986              }
     
    18072989            }
    18082990          } else if (iOriginalRow<numberTotal) {
    1809             longDouble value = diagonal[iOriginalRow];
    1810             if (fabs(value)>1.0e-100) {
     2991            CoinWorkDouble value = diagonal[iOriginalRow];
     2992            if (CoinAbs(value)>1.0e-100) {
    18112993              value = 1.0/value;
    1812               largest = CoinMax(largest,fabs(value));
     2994              largest = CoinMax(largest,CoinAbs(value));
    18132995            } else {
    18142996              value = 1.0e100;
     
    18283010              if (jNewRow>iRow) {
    18293011                work[jNewRow]=elementByRow[j];
    1830                 largest = CoinMax(largest,fabs(elementByRow[j]));
     3012                largest = CoinMax(largest,CoinAbs(elementByRow[j]));
    18313013              }
    18323014            }
     
    18573039            CoinBigIndex j;
    18583040            iColumn=iOriginalRow;
    1859             longDouble value = diagonal[iColumn];
    1860             if (fabs(value)>1.0e-100) {
     3041            CoinWorkDouble value = diagonal[iColumn];
     3042            if (CoinAbs(value)>1.0e-100) {
    18613043              value = 1.0/value;
    18623044              for (j=columnQuadraticStart[iColumn];
     
    18703052                }
    18713053              }
    1872               largest = CoinMax(largest,fabs(value));
     3054              largest = CoinMax(largest,CoinAbs(value));
    18733055              diagonal_[iRow] = -value;
    18743056              CoinBigIndex start=columnStart[iColumn];
     
    18793061                if (kRow>iRow) {
    18803062                  work[kRow]=element[j];
    1881                   largest = CoinMax(largest,fabs(element[j]));
     3063                  largest = CoinMax(largest,CoinAbs(element[j]));
    18823064                }
    18833065              }
     
    18863068            }
    18873069          } else if (iOriginalRow<numberTotal) {
    1888             longDouble value = diagonal[iOriginalRow];
    1889             if (fabs(value)>1.0e-100) {
     3070            CoinWorkDouble value = diagonal[iOriginalRow];
     3071            if (CoinAbs(value)>1.0e-100) {
    18903072              value = 1.0/value;
    1891               largest = CoinMax(largest,fabs(value));
     3073              largest = CoinMax(largest,CoinAbs(value));
    18923074            } else {
    18933075              value = 1.0e100;
     
    19073089              if (jNewRow>iRow) {
    19083090                work[jNewRow]=elementByRow[j];
    1909                 largest = CoinMax(largest,fabs(elementByRow[j]));
     3091                largest = CoinMax(largest,CoinAbs(elementByRow[j]));
    19103092              }
    19113093            }
     
    19313113          longDouble * put = sparseFactor_+choleskyStart_[iColumn];
    19323114          int * which = choleskyRow_+indexStart_[iColumn];
    1933           longDouble value = diagonal[iColumn];
    1934           if (fabs(value)>1.0e-100) {
     3115          CoinWorkDouble value = diagonal[iColumn];
     3116          if (CoinAbs(value)>1.0e-100) {
    19353117            value = 1.0/value;
    1936             largest = CoinMax(largest,fabs(value));
     3118            largest = CoinMax(largest,CoinAbs(value));
    19373119            diagonal_[iColumn] = -value;
    19383120            CoinBigIndex start=columnStart[iColumn];
     
    19423124              //sparseFactor_[numberElements++]=element[j];
    19433125              work[row[j]+numberTotal]=element[j];
    1944               largest = CoinMax(largest,fabs(element[j]));
     3126              largest = CoinMax(largest,CoinAbs(element[j]));
    19453127            }
    19463128          } else {
     
    19653147          int * which = choleskyRow_+indexStart_[iColumn];
    19663148          int number = choleskyStart_[iColumn+1]-choleskyStart_[iColumn];
    1967           longDouble value = diagonal[iColumn];
     3149          CoinWorkDouble value = diagonal[iColumn];
    19683150          CoinBigIndex j;
    1969           if (fabs(value)>1.0e-100) {
     3151          if (CoinAbs(value)>1.0e-100) {
    19703152            value = 1.0/value;
    19713153            for (j=columnQuadraticStart[iColumn];
     
    19783160              }
    19793161            }
    1980             largest = CoinMax(largest,fabs(value));
     3162            largest = CoinMax(largest,CoinAbs(value));
    19813163            diagonal_[iColumn] = -value;
    19823164            CoinBigIndex start=columnStart[iColumn];
     
    19843166            for (j=start;j<end;j++) {
    19853167              work[row[j]+numberTotal]=element[j];
    1986               largest = CoinMax(largest,fabs(element[j]));
     3168              largest = CoinMax(largest,CoinAbs(element[j]));
    19873169            }
    19883170            for (j=0;j<number;j++) {
     
    20053187        longDouble * put = sparseFactor_+choleskyStart_[iColumn];
    20063188        int * which = choleskyRow_+indexStart_[iColumn];
    2007         longDouble value = diagonal[iColumn];
    2008         if (fabs(value)>1.0e-100) {
     3189        CoinWorkDouble value = diagonal[iColumn];
     3190        if (CoinAbs(value)>1.0e-100) {
    20093191          value = 1.0/value;
    2010           largest = CoinMax(largest,fabs(value));
     3192          largest = CoinMax(largest,CoinAbs(value));
    20113193        } else {
    20123194          value = 1.0e100;
     
    20233205      }
    20243206      // Finish diagonal
    2025       longDouble delta2 = model_->delta(); // add delta*delta to bottom
     3207      CoinWorkDouble delta2 = model_->delta(); // add delta*delta to bottom
    20263208      delta2 *= delta2;
    20273209      for (iRow=0;iRow<numberRowsModel;iRow++) {
     
    20373219    //check sizes
    20383220    largest*=1.0e-20;
    2039     largest = CoinMin(largest,1.0e-11);
     3221    largest = CoinMin(largest,CHOL_SMALL_VALUE);
    20403222    doubleParameters_[10]=CoinMax(1.0e-20,largest);
    20413223    integerParameters_[20]=0;
     
    20563238    // Should save adjustments in ..R_
    20573239    int n1=0,n2=0;
    2058     double * primalR = model_->primalR();
    2059     double * dualR = model_->dualR();
     3240    CoinWorkDouble * primalR = model_->primalR();
     3241    CoinWorkDouble * dualR = model_->dualR();
    20603242    for (iRow=0;iRow<numberTotal;iRow++) {
    20613243      if (rowsDropped2[iRow]) {
     
    20933275ClpCholeskyBase::factorizePart2(int * rowsDropped)
    20943276{
    2095   longDouble largest=doubleParameters_[3];
    2096   longDouble smallest=doubleParameters_[4];
     3277  CoinWorkDouble largest=doubleParameters_[3];
     3278  CoinWorkDouble smallest=doubleParameters_[4];
    20973279  int numberDropped=integerParameters_[20];
    20983280  // probably done before
     
    21593341      for ( jRow=lastRow;jRow<iRow;jRow++) {
    21603342        int jCount = jRow-lastRow;
    2161         longDouble diagonalValue = diagonal_[jRow];
     3343        CoinWorkDouble diagonalValue = diagonal_[jRow];
    21623344        CoinBigIndex start=choleskyStart_[jRow];
    21633345        CoinBigIndex end=choleskyStart_[jRow+1];
     
    21653347          jCount--;
    21663348          CoinBigIndex get = choleskyStart_[kRow]+jCount;
    2167           longDouble a_jk = sparseFactor_[get];
    2168           longDouble value1 = d[kRow]*a_jk;
     3349          CoinWorkDouble a_jk = sparseFactor_[get];
     3350          CoinWorkDouble value1 = d[kRow]*a_jk;
    21693351          diagonalValue -= a_jk*value1;
    21703352          for (CoinBigIndex j=start;j<end;j++)
     
    22223404    }
    22233405    // for each column L[*,kRow] that affects L[*,iRow]
    2224     longDouble diagonalValue=diagonal_[iRow];
     3406    CoinWorkDouble diagonalValue=diagonal_[iRow];
    22253407    int nextRow = link_[iRow];
    22263408    int kRow=0;
     
    22343416      CoinBigIndex end=choleskyStart_[kRow+1];
    22353417      assert(k<end);
    2236       longDouble a_ik=sparseFactor_[k++];
    2237       longDouble value1=d[kRow]*a_ik;
     3418      CoinWorkDouble a_ik=sparseFactor_[k++];
     3419      CoinWorkDouble value1=d[kRow]*a_ik;
    22383420      // update first
    22393421      first[kRow]=k;
     
    22593441            CoinBigIndex j=first[kkRow];
    22603442            //int iiRow = choleskyRow_[j+indexStart_[kkRow]-choleskyStart_[kkRow]];
    2261             longDouble a = sparseFactor_[j];
    2262             longDouble dValue = d[kkRow]*a;
     3443            CoinWorkDouble a = sparseFactor_[j];
     3444            CoinWorkDouble dValue = d[kkRow]*a;
    22633445            diagonalValue -= a*dValue;
    22643446            work[kkRow]=dValue;
     
    22713453          for (int i=0;i<length;i++) {
    22723454            int lRow = choleskyRow_[currentIndex++];
    2273             longDouble t0 = work[lRow];
     3455            CoinWorkDouble t0 = work[lRow];
    22743456            for (int kkRow=kRow;kkRow<last;kkRow++) {
    22753457              CoinBigIndex j = first[kkRow]+i;
     
    23403522        for (CoinBigIndex j=start;j<end;j++) {
    23413523          int jRow = choleskyRow_[j+offset];
    2342           longDouble value = sparseFactor_[j]-work[jRow];
     3524          CoinWorkDouble value = sparseFactor_[j]-work[jRow];
    23433525          work[jRow]=0.0;
    23443526          sparseFactor_[j]= diagonalValue*value;
     
    23503532    // do dense
    23513533    // update dense part
    2352     updateDense(d,work,first);
     3534    updateDense(d,/*work,*/first);
    23533535    ClpCholeskyDense dense;
    23543536    // just borrow space
     
    23713553    dense.setIntegerParameter(20,0);
    23723554    dense.setIntegerParameter(34,firstPositive);
     3555    dense.setModel(model_);
    23733556    dense.factorizePart2(dropped);
    23743557    largest=dense.getDoubleParameter(3);
     
    23873570}
    23883571// Updates dense part (broken out for profiling)
    2389 void ClpCholeskyBase::updateDense(longDouble * d, longDouble * work, int * first)
     3572void ClpCholeskyBase::updateDense(longDouble * d, /*longDouble * work,*/ int * first)
    23903573{
    23913574  for (int iRow=0;iRow<firstDense_;iRow++) {
     
    23953578      CoinBigIndex offset = indexStart_[iRow]-choleskyStart_[iRow];
    23963579      if (clique_[iRow]<2) {
    2397         longDouble dValue=d[iRow];
     3580        CoinWorkDouble dValue=d[iRow];
    23983581        for (CoinBigIndex k=start;k<end;k++) {
    23993582          int kRow = choleskyRow_[k+offset];
    24003583          assert(kRow>=firstDense_);
    2401           longDouble a_ik=sparseFactor_[k];
    2402           longDouble value1=dValue*a_ik;
     3584          CoinWorkDouble a_ik=sparseFactor_[k];
     3585          CoinWorkDouble value1=dValue*a_ik;
    24033586          diagonal_[kRow] -= value1*a_ik;
    24043587          CoinBigIndex base = choleskyStart_[kRow]-kRow-1;
    24053588          for (CoinBigIndex j=k+1;j<end;j++) {
    24063589            int jRow=choleskyRow_[j+offset];
    2407             longDouble a_jk = sparseFactor_[j];
     3590            CoinWorkDouble a_jk = sparseFactor_[j];
    24083591            sparseFactor_[base+jRow] -= a_jk*value1;
    24093592          }
     
    24113594      } else if (clique_[iRow]<3) {
    24123595        // do as pair
    2413         longDouble dValue0=d[iRow];
    2414         longDouble dValue1=d[iRow+1];
     3596        CoinWorkDouble dValue0=d[iRow];
     3597        CoinWorkDouble dValue1=d[iRow+1];
    24153598        int offset1 = first[iRow+1]-start;
    24163599        // skip row
     
    24193602          int kRow = choleskyRow_[k+offset];
    24203603          assert(kRow>=firstDense_);
    2421           longDouble a_ik0=sparseFactor_[k];
    2422           longDouble value0=dValue0*a_ik0;
    2423           longDouble a_ik1=sparseFactor_[k+offset1];
    2424           longDouble value1=dValue1*a_ik1;
     3604          CoinWorkDouble a_ik0=sparseFactor_[k];
     3605          CoinWorkDouble value0=dValue0*a_ik0;
     3606          CoinWorkDouble a_ik1=sparseFactor_[k+offset1];
     3607          CoinWorkDouble value1=dValue1*a_ik1;
    24253608          diagonal_[kRow] -= value0*a_ik0 + value1*a_ik1;
    24263609          CoinBigIndex base = choleskyStart_[kRow]-kRow-1;
    24273610          for (CoinBigIndex j=k+1;j<end;j++) {
    24283611            int jRow=choleskyRow_[j+offset];
    2429             longDouble a_jk0 = sparseFactor_[j];
    2430             longDouble a_jk1 = sparseFactor_[j+offset1];
     3612            CoinWorkDouble a_jk0 = sparseFactor_[j];
     3613            CoinWorkDouble a_jk1 = sparseFactor_[j+offset1];
    24313614            sparseFactor_[base+jRow] -= a_jk0*value0+a_jk1*value1;
    24323615          }
     
    24413624        // maybe later get fancy on big cliques and do transpose copy
    24423625        // seems only worth going to 3 on Intel
    2443         longDouble dValue0=d[iRow];
    2444         longDouble dValue1=d[iRow+1];
    2445         longDouble dValue2=d[iRow+2];
     3626        CoinWorkDouble dValue0=d[iRow];
     3627        CoinWorkDouble dValue1=d[iRow+1];
     3628        CoinWorkDouble dValue2=d[iRow+2];
    24463629        // get offsets and skip rows
    24473630        int offset1 = first[++iRow]-start;
     
    24503633          int kRow = choleskyRow_[k+offset];
    24513634          assert(kRow>=firstDense_);
    2452           double diagonalValue=diagonal_[kRow];
    2453           longDouble a_ik0=sparseFactor_[k];
    2454           longDouble value0=dValue0*a_ik0;
    2455           longDouble a_ik1=sparseFactor_[k+offset1];
    2456           longDouble value1=dValue1*a_ik1;
    2457           longDouble a_ik2=sparseFactor_[k+offset2];
    2458           longDouble value2=dValue2*a_ik2;
     3635          CoinWorkDouble diagonalValue=diagonal_[kRow];
     3636          CoinWorkDouble a_ik0=sparseFactor_[k];
     3637          CoinWorkDouble value0=dValue0*a_ik0;
     3638          CoinWorkDouble a_ik1=sparseFactor_[k+offset1];
     3639          CoinWorkDouble value1=dValue1*a_ik1;
     3640          CoinWorkDouble a_ik2=sparseFactor_[k+offset2];
     3641          CoinWorkDouble value2=dValue2*a_ik2;
    24593642          CoinBigIndex base = choleskyStart_[kRow]-kRow-1;
    24603643          diagonal_[kRow] = diagonalValue - value0*a_ik0 - value1*a_ik1 - value2*a_ik2;
    24613644          for (CoinBigIndex j=k+1;j<end;j++) {
    24623645            int jRow=choleskyRow_[j+offset];
    2463             longDouble a_jk0 = sparseFactor_[j];
    2464             longDouble a_jk1 = sparseFactor_[j+offset1];
    2465             longDouble a_jk2 = sparseFactor_[j+offset2];
     3646            CoinWorkDouble a_jk0 = sparseFactor_[j];
     3647            CoinWorkDouble a_jk1 = sparseFactor_[j+offset1];
     3648            CoinWorkDouble a_jk2 = sparseFactor_[j+offset2];
    24663649            sparseFactor_[base+jRow] -= a_jk0*value0+a_jk1*value1+a_jk2*value2;
    24673650          }
     
    24723655        // maybe later get fancy on big cliques and do transpose copy
    24733656        // maybe only worth going to 3 on Intel (but may have hidden registers)
    2474         longDouble dValue0=d[iRow];
    2475         longDouble dValue1=d[iRow+1];
    2476         longDouble dValue2=d[iRow+2];
    2477         longDouble dValue3=d[iRow+3];
     3657        CoinWorkDouble dValue0=d[iRow];
     3658        CoinWorkDouble dValue1=d[iRow+1];
     3659        CoinWorkDouble dValue2=d[iRow+2];
     3660        CoinWorkDouble dValue3=d[iRow+3];
    24783661        // get offsets and skip rows
    24793662        int offset1 = first[++iRow]-start;
     
    24833666          int kRow = choleskyRow_[k+offset];
    24843667          assert(kRow>=firstDense_);
    2485           double diagonalValue=diagonal_[kRow];
    2486           longDouble a_ik0=sparseFactor_[k];
    2487           longDouble value0=dValue0*a_ik0;
    2488           longDouble a_ik1=sparseFactor_[k+offset1];
    2489           longDouble value1=dValue1*a_ik1;
    2490           longDouble a_ik2=sparseFactor_[k+offset2];
    2491           longDouble value2=dValue2*a_ik2;
    2492           longDouble a_ik3=sparseFactor_[k+offset3];
    2493           longDouble value3=dValue3*a_ik3;
     3668          CoinWorkDouble diagonalValue=diagonal_[kRow];
     3669          CoinWorkDouble a_ik0=sparseFactor_[k];
     3670          CoinWorkDouble value0=dValue0*a_ik0;
     3671          CoinWorkDouble a_ik1=sparseFactor_[k+offset1];
     3672          CoinWorkDouble value1=dValue1*a_ik1;
     3673          CoinWorkDouble a_ik2=sparseFactor_[k+offset2];
     3674          CoinWorkDouble value2=dValue2*a_ik2;
     3675          CoinWorkDouble a_ik3=sparseFactor_[k+offset3];
     3676          CoinWorkDouble value3=dValue3*a_ik3;
    24943677          CoinBigIndex base = choleskyStart_[kRow]-kRow-1;
    24953678          diagonal_[kRow] = diagonalValue - (value0*a_ik0 + value1*a_ik1 + value2*a_ik2+value3*a_ik3);
    24963679          for (CoinBigIndex j=k+1;j<end;j++) {
    24973680            int jRow=choleskyRow_[j+offset];
    2498             longDouble a_jk0 = sparseFactor_[j];
    2499             longDouble a_jk1 = sparseFactor_[j+offset1];
    2500             longDouble a_jk2 = sparseFactor_[j+offset2];
    2501             longDouble a_jk3 = sparseFactor_[j+offset3];
     3681            CoinWorkDouble a_jk0 = sparseFactor_[j];
     3682            CoinWorkDouble a_jk1 = sparseFactor_[j+offset1];
     3683            CoinWorkDouble a_jk2 = sparseFactor_[j+offset2];
     3684            CoinWorkDouble a_jk3 = sparseFactor_[j+offset3];
    25023685            sparseFactor_[base+jRow] -= a_jk0*value0+a_jk1*value1+a_jk2*value2+a_jk3*value3;
    25033686          }
     
    25103693/* Uses factorization to solve. */
    25113694void
    2512 ClpCholeskyBase::solve (double * region)
     3695ClpCholeskyBase::solve (CoinWorkDouble * region)
    25133696{
    25143697  if (!whichDense_) {
     
    25203703    // do change;
    25213704    int numberDense = dense_->numberRows();
    2522     double * change = new double[numberDense];
     3705    CoinWorkDouble * change = new CoinWorkDouble[numberDense];
    25233706    for (i=0;i<numberDense;i++) {
    25243707      const longDouble * a = denseColumn_+i*numberRows_;
    2525       longDouble value =0.0;
     3708      CoinWorkDouble value =0.0;
    25263709      for (int iRow=0;iRow<numberRows_;iRow++)
    25273710        value += a[iRow]*region[iRow];
     
    25323715    for (i=0;i<numberDense;i++) {
    25333716      const longDouble * a = denseColumn_+i*numberRows_;
    2534       longDouble value = change[i];
     3717      CoinWorkDouble value = change[i];
    25353718      for (int iRow=0;iRow<numberRows_;iRow++)
    25363719        region[iRow] -= value*a[iRow];
     
    25453728*/
    25463729void
    2547 ClpCholeskyBase::solve(double * region, int type)
     3730ClpCholeskyBase::solve(CoinWorkDouble * region, int type)
    25483731{
    25493732#ifdef CLP_DEBUG
    25503733  double * regionX=NULL;
    2551   if (sizeof(longWork)!=sizeof(double)&&type==3) {
     3734  if (sizeof(CoinWorkDouble)!=sizeof(double)&&type==3) {
    25523735    regionX=ClpCopyOfArray(region,numberRows_);
    25533736  }
    25543737#endif
    2555   longWork * work = reinterpret_cast<longWork *> (workDouble_);
     3738  CoinWorkDouble * work = reinterpret_cast<CoinWorkDouble *> (workDouble_);
    25563739  int i;
    25573740  CoinBigIndex j;
     
    25633746  case 1:
    25643747    for (i=0;i<numberRows_;i++) {
    2565       longDouble value=work[i];
     3748      CoinWorkDouble value=work[i];
    25663749      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
    25673750      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
     
    25783761    for (i=numberRows_-1;i>=0;i--) {
    25793762      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
    2580       longDouble value=work[i]*diagonal_[i];
     3763      CoinWorkDouble value=work[i]*diagonal_[i];
    25813764      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
    25823765        int iRow = choleskyRow_[j+offset];
     
    25913774    for (i=0;i<firstDense_;i++) {
    25923775      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
    2593       longDouble value=work[i];
     3776      CoinWorkDouble value=work[i];
    25943777      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
    25953778        int iRow = choleskyRow_[j+offset];
     
    26033786      int nDense = numberRows_-firstDense_;
    26043787      dense.reserveSpace(this,nDense);
    2605 #if CLP_LONG_CHOLESKY!=1
    2606       dense.solveLong(work+firstDense_);
    2607 #else
    2608       dense.solveLongWork(work+firstDense_);
    2609 #endif
     3788      dense.solve(work+firstDense_);
    26103789      for (i=numberRows_-1;i>=firstDense_;i--) {
    2611         longDouble value=work[i];
     3790        CoinWorkDouble value=work[i];
    26123791        int iRow = permute_[i];
    26133792        region[iRow]=value;
     
    26163795    for (i=firstDense_-1;i>=0;i--) {
    26173796      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
    2618       longDouble value=work[i]*diagonal_[i];
     3797      CoinWorkDouble value=work[i]*diagonal_[i];
    26193798      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
    26203799        int iRow = choleskyRow_[j+offset];
     
    26343813    double largestO=0.0;
    26353814    for (i=0;i<numberRows_;i++) {
    2636       largestO = CoinMax(largestO,fabs(regionX[i]));
     3815      largestO = CoinMax(largestO,CoinAbs(regionX[i]));
    26373816    }
    26383817    for (i=0;i<numberRows_;i++) {
     
    26423821    for (i=0;i<firstDense_;i++) {
    26433822      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
    2644       longDouble value=work[i];
     3823      CoinWorkDouble value=work[i];
    26453824      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
    26463825        int iRow = choleskyRow_[j+offset];
     
    26543833      int nDense = numberRows_-firstDense_;
    26553834      dense.reserveSpace(this,nDense);
    2656       dense.solveLong(work+firstDense_);
     3835      dense.solve(work+firstDense_);
    26573836      for (i=numberRows_-1;i>=firstDense_;i--) {
    2658         longDouble value=work[i];
     3837        CoinWorkDouble value=work[i];
    26593838        int iRow = permute_[i];
    26603839        regionX[iRow]=value;
     
    26633842    for (i=firstDense_-1;i>=0;i--) {
    26643843      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
    2665       longDouble value=work[i]*diagonal_[i];
     3844      CoinWorkDouble value=work[i]*diagonal_[i];
    26663845      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
    26673846        int iRow = choleskyRow_[j+offset];
     
    26753854    double largestV=0.0;
    26763855    for (i=0;i<numberRows_;i++) {
    2677       largest = CoinMax(largest,fabs(region[i]-regionX[i]));
    2678       largestV = CoinMax(largestV,fabs(region[i]));
     3856      largest = CoinMax(largest,CoinAbs(region[i]-regionX[i]));
     3857      largestV = CoinMax(largestV,CoinAbs(region[i]));
    26793858    }
    26803859    printf("largest difference %g, largest %g, largest original %g\n",
     
    26843863#endif
    26853864}
    2686 /* solve - 1 just first half, 2 just second half - 3 both.
    2687    If 1 and 2 then diagonal has sqrt of inverse otherwise inverse
    2688 */
     3865#if 0 //CLP_LONG_CHOLESKY
     3866/* Uses factorization to solve. */
    26893867void
    2690 ClpCholeskyBase::solveLong(longDouble * region, int type)
     3868ClpCholeskyBase::solve (CoinWorkDouble * region)
    26913869{
     3870  assert (!whichDense_) ;
     3871  CoinWorkDouble * work = reinterpret_cast<CoinWorkDouble *> (workDouble_);
    26923872  int i;
    26933873  CoinBigIndex j;
    26943874  for (i=0;i<numberRows_;i++) {
    26953875    int iRow = permute_[i];
    2696     workDouble_[i] = region[iRow];
    2697   }
    2698   switch (type) {
    2699   case 1:
    2700     for (i=0;i<numberRows_;i++) {
    2701       longDouble value=workDouble_[i];
    2702       CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
    2703       for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
    2704         int iRow = choleskyRow_[j+offset];
    2705         workDouble_[iRow] -= sparseFactor_[j]*value;
    2706       }
    2707     }
    2708     for (i=0;i<numberRows_;i++) {
    2709       int iRow = permute_[i];
    2710       region[iRow]=workDouble_[i]*diagonal_[i];
    2711     }
    2712     break;
    2713   case 2:
    2714     for (i=numberRows_-1;i>=0;i--) {
    2715       CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
    2716       longDouble value=workDouble_[i]*diagonal_[i];
    2717       for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
    2718         int iRow = choleskyRow_[j+offset];
    2719         value -= sparseFactor_[j]*workDouble_[iRow];
    2720       }
    2721       workDouble_[i]=value;
     3876    work[i] = region[iRow];
     3877  }
     3878  for (i=0;i<firstDense_;i++) {
     3879    CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
     3880    CoinWorkDouble value=work[i];
     3881    for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
     3882      int iRow = choleskyRow_[j+offset];
     3883      work[iRow] -= sparseFactor_[j]*value;
     3884    }
     3885  }
     3886  if (firstDense_<numberRows_) {
     3887    // do dense
     3888    ClpCholeskyDense dense;
     3889    // just borrow space
     3890    int nDense = numberRows_-firstDense_;
     3891    dense.reserveSpace(this,nDense);
     3892    dense.solve(work+firstDense_);
     3893    for (i=numberRows_-1;i>=firstDense_;i--) {
     3894      CoinWorkDouble value=work[i];
    27223895      int iRow = permute_[i];
    27233896      region[iRow]=value;
    27243897    }
    2725     break;
    2726   case 3:
    2727     for (i=0;i<firstDense_;i++) {
    2728       CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
    2729       longDouble value=workDouble_[i];
    2730       for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
    2731         int iRow = choleskyRow_[j+offset];
    2732         workDouble_[iRow] -= sparseFactor_[j]*value;
    2733       }
    2734     }
    2735     if (firstDense_<numberRows_) {
    2736       // do dense
    2737       ClpCholeskyDense dense;
    2738       // just borrow space
    2739       int nDense = numberRows_-firstDense_;
    2740       dense.reserveSpace(this,nDense);
    2741       dense.solveLong(workDouble_+firstDense_);
    2742       for (i=numberRows_-1;i>=firstDense_;i--) {
    2743         longDouble value=workDouble_[i];
    2744         int iRow = permute_[i];
    2745         region[iRow]=value;
    2746       }
    2747     }
    2748     for (i=firstDense_-1;i>=0;i--) {
    2749       CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
    2750       longDouble value=workDouble_[i]*diagonal_[i];
    2751       for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
    2752         int iRow = choleskyRow_[j+offset];
    2753         value -= sparseFactor_[j]*workDouble_[iRow];
    2754       }
    2755       workDouble_[i]=value;
    2756       int iRow = permute_[i];
    2757       region[iRow]=value;
    2758     }
    2759     break;
     3898  }
     3899  for (i=firstDense_-1;i>=0;i--) {
     3900    CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
     3901    CoinWorkDouble value=work[i]*diagonal_[i];
     3902    for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
     3903      int iRow = choleskyRow_[j+offset];
     3904      value -= sparseFactor_[j]*work[iRow];
     3905    }
     3906    work[i]=value;
     3907    int iRow = permute_[i];
     3908    region[iRow]=value;
    27603909  }
    27613910}
    2762 
     3911#endif
  • stable/1.11/Clp/src/ClpCholeskyBase.hpp

    r1055 r1458  
     1/* $Id$ */
    12// Copyright (C) 2003, International Business Machines
    23// Corporation and others.  All Rights Reserved.
     
    67#include "CoinPragma.hpp"
    78#include "CoinFinite.hpp"
     9//#define CLP_LONG_CHOLESKY 0
     10#ifndef CLP_LONG_CHOLESKY
    811#define CLP_LONG_CHOLESKY 0
     12#endif
     13/* valid combinations are
     14   CLP_LONG_CHOLESKY 0 and COIN_LONG_WORK 0
     15   CLP_LONG_CHOLESKY 1 and COIN_LONG_WORK 1
     16   CLP_LONG_CHOLESKY 2 and COIN_LONG_WORK 1
     17*/
     18#if COIN_LONG_WORK==0
     19#if CLP_LONG_CHOLESKY>0
     20#define CHOLESKY_BAD_COMBINATION
     21#endif
     22#else
     23#if CLP_LONG_CHOLESKY==0
     24#define CHOLESKY_BAD_COMBINATION
     25#endif
     26#endif
     27#ifdef CHOLESKY_BAD_COMBINATION
     28#  warning("Bad combination of CLP_LONG_CHOLESKY and COIN_BIG_DOUBLE/COIN_LONG_WORK");
     29"Bad combination of CLP_LONG_CHOLESKY and COIN_LONG_WORK"
     30#endif
    931#if CLP_LONG_CHOLESKY>1
    1032typedef long double longDouble;
    11 typedef long double longWork;
     33#define CHOL_SMALL_VALUE 1.0e-15
    1234#elif CLP_LONG_CHOLESKY==1
    1335typedef double longDouble;
    14 typedef long double longWork;
     36#define CHOL_SMALL_VALUE 1.0e-11
    1537#else
    1638typedef double longDouble;
    17 typedef double longWork;
     39#define CHOL_SMALL_VALUE 1.0e-11
    1840#endif
    1941class ClpInterior;
     
    4668  /** Factorize - filling in rowsDropped and returning number dropped.
    4769      If return code negative then out of memory */
    48   virtual int factorize(const double * diagonal, int * rowsDropped) ;
     70  virtual int factorize(const CoinWorkDouble * diagonal, int * rowsDropped) ;
    4971  /** Uses factorization to solve. */
    50   virtual void solve (double * region) ;
     72  virtual void solve (CoinWorkDouble * region) ;
    5173  /** Uses factorization to solve. - given as if KKT.
    5274   region1 is rows+columns, region2 is rows */
    53   virtual void solveKKT (double * region1, double * region2, const double * diagonal,
    54                          double diagonalScaleFactor);
     75  virtual void solveKKT (CoinWorkDouble * region1, CoinWorkDouble * region2, const CoinWorkDouble * diagonal,
     76                         CoinWorkDouble diagonalScaleFactor);
     77private:
     78  /// AMD ordering
     79  int orderAMD();
     80public:
    5581  //@}
    5682
     
    141167protected:
    142168  /// Sets type
    143   void setType(int type) {type_=type;}
     169  inline void setType(int type) {type_=type;}
     170  /// model.
     171  inline void setModel(ClpInterior * model)
     172  { model_=model;}
    144173   //@}
    145174   
     
    162191  If 1 and 2 then diagonal has sqrt of inverse otherwise inverse
    163192  */
    164   void solve(double * region, int type);
    165   void solveLong(longDouble * region, int type);
     193  void solve(CoinWorkDouble * region, int type);
    166194  /// Forms ADAT - returns nonzero if not enough memory
    167195  int preOrder(bool lowerTriangular, bool includeDiagonal, bool doKKT);
    168196  /// Updates dense part (broken out for profiling)
    169   void updateDense(longDouble * d, longDouble * work, int * first);
     197  void updateDense(longDouble * d, /*longDouble * work,*/ int * first);
    170198  //@}
    171199   
  • stable/1.11/Clp/src/ClpCholeskyDense.cpp

    r1404 r1458  
    1 // Copyright (C) 2003, International Business Machines
    2 // Corporation and others.  All Rights Reserved.
    3 
    4 
    5 
     1/* $Id$ */
     2/* Copyright (C) 2003, International Business Machines Corporation
     3   and others.  All Rights Reserved. */
    64#include "CoinPragma.hpp"
    75#include "CoinHelperFunctions.hpp"
     6#include "ClpHelperFunctions.hpp"
    87
    98#include "ClpInterior.hpp"
     
    1211#include "ClpQuadraticObjective.hpp"
    1312
    14 //#############################################################################
    15 // Constructors / Destructor / Assignment
    16 //#############################################################################
    17 
    18 //-------------------------------------------------------------------
    19 // Default Constructor
    20 //-------------------------------------------------------------------
     13/*#############################################################################*/
     14/* Constructors / Destructor / Assignment*/
     15/*#############################################################################*/
     16
     17/*-------------------------------------------------------------------*/
     18/* Default Constructor */
     19/*-------------------------------------------------------------------*/
    2120ClpCholeskyDense::ClpCholeskyDense ()
    2221  : ClpCholeskyBase(),
     
    2625}
    2726
    28 //-------------------------------------------------------------------
    29 // Copy constructor
    30 //-------------------------------------------------------------------
     27/*-------------------------------------------------------------------*/
     28/* Copy constructor */
     29/*-------------------------------------------------------------------*/
    3130ClpCholeskyDense::ClpCholeskyDense (const ClpCholeskyDense & rhs)
    3231  : ClpCholeskyBase(rhs),
    3332    borrowSpace_(rhs.borrowSpace_)
    3433{
    35   assert(!rhs.borrowSpace_||!rhs.sizeFactor_); // can't do if borrowing space
    36 }
    37 
    38 
    39 //-------------------------------------------------------------------
    40 // Destructor
    41 //-------------------------------------------------------------------
     34  assert(!rhs.borrowSpace_||!rhs.sizeFactor_); /* can't do if borrowing space*/
     35}
     36
     37
     38/*-------------------------------------------------------------------*/
     39/* Destructor */
     40/*-------------------------------------------------------------------*/
    4241ClpCholeskyDense::~ClpCholeskyDense ()
    4342{
    4443  if (borrowSpace_) {
    45     // set NULL
     44    /* set NULL*/
    4645    sparseFactor_=NULL;
    4746    workDouble_=NULL;
     
    5049}
    5150
    52 //----------------------------------------------------------------
    53 // Assignment operator
    54 //-------------------------------------------------------------------
     51/*----------------------------------------------------------------*/
     52/* Assignment operator */
     53/*-------------------------------------------------------------------*/
    5554ClpCholeskyDense &
    5655ClpCholeskyDense::operator=(const ClpCholeskyDense& rhs)
    5756{
    5857  if (this != &rhs) {
    59     assert(!rhs.borrowSpace_||!rhs.sizeFactor_); // can't do if borrowing space
     58    assert(!rhs.borrowSpace_||!rhs.sizeFactor_); /* can't do if borrowing space*/
    6059    ClpCholeskyBase::operator=(rhs);
    6160    borrowSpace_=rhs.borrowSpace_;
     
    6362  return *this;
    6463}
    65 //-------------------------------------------------------------------
    66 // Clone
    67 //-------------------------------------------------------------------
     64/*-------------------------------------------------------------------*/
     65/* Clone*/
     66/*-------------------------------------------------------------------*/
    6867ClpCholeskyBase * ClpCholeskyDense::clone() const
    6968{
    7069  return new ClpCholeskyDense(*this);
    7170}
    72 // If not power of 2 then need to redo a bit
     71/* If not power of 2 then need to redo a bit*/
    7372#define BLOCK 16
    7473#define BLOCKSHIFT 4
    75 // Block unroll if power of 2 and at least 8
     74/* Block unroll if power of 2 and at least 8*/
    7675#define BLOCKUNROLL
    7776
     
    8786  numberRows_ = numberRows;
    8887  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
    89   // allow one stripe extra
     88  /* allow one stripe extra*/
    9089  numberBlocks = numberBlocks + ((numberBlocks*(numberBlocks+1))/2);
    9190  sizeFactor_=numberBlocks*BLOCKSQ;
    92   //#define CHOL_COMPARE
     91  /*#define CHOL_COMPARE*/
    9392#ifdef CHOL_COMPARE 
    9493  sizeFactor_ += 95000;
     
    115114{
    116115 int numberBlocks = (numberRows+BLOCK-1)>>BLOCKSHIFT;
    117   // allow one stripe extra
     116 /* allow one stripe extra*/
    118117  numberBlocks = numberBlocks + ((numberBlocks*(numberBlocks+1))/2);
    119118  CoinBigIndex sizeFactor=numberBlocks*BLOCKSQ;
     
    151150/* Factorize - filling in rowsDropped and returning number dropped */
    152151int
    153 ClpCholeskyDense::factorize(const double * diagonal, int * rowsDropped)
     152ClpCholeskyDense::factorize(const CoinWorkDouble * diagonal, int * rowsDropped)
    154153{
    155154  assert (!borrowSpace_);
     
    164163  int numberColumns=model_->clpMatrix()->getNumCols();
    165164  CoinZeroN(sparseFactor_,sizeFactor_);
    166   //perturbation
    167   longDouble perturbation=model_->diagonalPerturbation()*model_->diagonalNorm();
     165  /*perturbation*/
     166  CoinWorkDouble perturbation=model_->diagonalPerturbation()*model_->diagonalNorm();
    168167  perturbation=perturbation*perturbation;
    169168  if (perturbation>1.0) {
    170169#ifdef COIN_DEVELOP
    171     //if (model_->model()->logLevel()&4)
     170    /*if (model_->model()->logLevel()&4) */
    172171      std::cout <<"large perturbation "<<perturbation<<std::endl;
    173172#endif
    174     perturbation=sqrt(perturbation);;
     173    perturbation=CoinSqrt(perturbation);;
    175174    perturbation=1.0;
    176175  }
    177176  int iRow;
    178177  int newDropped=0;
    179   double largest=1.0;
    180   double smallest=COIN_DBL_MAX;
    181   longDouble delta2 = model_->delta(); // add delta*delta to diagonal
     178  CoinWorkDouble largest=1.0;
     179  CoinWorkDouble smallest=COIN_DBL_MAX;
     180  CoinWorkDouble delta2 = model_->delta(); /* add delta*delta to diagonal*/
    182181  delta2 *= delta2;
    183182  if (!doKKT_) {
    184183    longDouble * work = sparseFactor_;
    185     work--; // skip diagonal
     184    work--; /* skip diagonal*/
    186185    int addOffset=numberRows_-1;
    187     const double * diagonalSlack = diagonal + numberColumns;
    188     // largest in initial matrix
    189     double largest2=1.0e-20;
     186    const CoinWorkDouble * diagonalSlack = diagonal + numberColumns;
     187    /* largest in initial matrix*/
     188    CoinWorkDouble largest2=1.0e-20;
    190189    for (iRow=0;iRow<numberRows_;iRow++) {
    191190      if (!rowsDropped_[iRow]) {
    192191        CoinBigIndex startRow=rowStart[iRow];
    193192        CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
    194         longDouble diagonalValue = diagonalSlack[iRow]+delta2;
     193        CoinWorkDouble diagonalValue = diagonalSlack[iRow]+delta2;
    195194        for (CoinBigIndex k=startRow;k<endRow;k++) {
    196195          int iColumn=column[k];
    197196          CoinBigIndex start=columnStart[iColumn];
    198197          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
    199           longDouble multiplier = diagonal[iColumn]*elementByRow[k];
     198          CoinWorkDouble multiplier = diagonal[iColumn]*elementByRow[k];
    200199          for (CoinBigIndex j=start;j<end;j++) {
    201200            int jRow=row[j];
    202201            if (!rowsDropped_[jRow]) {
    203202              if (jRow>iRow) {
    204                 longDouble value=element[j]*multiplier;
     203                CoinWorkDouble value=element[j]*multiplier;
    205204                work[jRow] += value;
    206205              } else if (jRow==iRow) {
    207                 longDouble value=element[j]*multiplier;
     206                CoinWorkDouble value=element[j]*multiplier;
    208207                diagonalValue += value;
    209208              }
     
    212211        }
    213212        for (int j=iRow+1;j<numberRows_;j++)
    214           largest2 = CoinMax(largest2,fabs(work[j]));
     213          largest2 = CoinMax(largest2,CoinAbs(work[j]));
    215214        diagonal_[iRow]=diagonalValue;
    216         largest2 = CoinMax(largest2,fabs(diagonalValue));
     215        largest2 = CoinMax(largest2,CoinAbs(diagonalValue));
    217216      } else {
    218         // dropped
     217        /* dropped*/
    219218        diagonal_[iRow]=1.0;
    220219      }
     
    222221      work += addOffset;
    223222    }
    224     //check sizes
     223    /*check sizes*/
    225224    largest2*=1.0e-20;
    226     largest = CoinMin(largest2,1.0e-11);
     225    largest = CoinMin(largest2,CHOL_SMALL_VALUE);
    227226    int numberDroppedBefore=0;
    228227    for (iRow=0;iRow<numberRows_;iRow++) {
    229228      int dropped=rowsDropped_[iRow];
    230       // Move to int array
     229      /* Move to int array*/
    231230      rowsDropped[iRow]=dropped;
    232231      if (!dropped) {
    233         longDouble diagonal = diagonal_[iRow];
     232        CoinWorkDouble diagonal = diagonal_[iRow];
    234233        if (diagonal>largest2) {
    235234          diagonal_[iRow]=diagonal+perturbation;
     
    245244    doubleParameters_[3]=0.0;
    246245    doubleParameters_[4]=COIN_DBL_MAX;
    247     integerParameters_[34]=0; // say all must be positive
     246    integerParameters_[34]=0; /* say all must be positive*/
    248247#ifdef CHOL_COMPARE 
    249248    if (numberRows_<200)
     
    257256      std::cout<<"Cholesky - largest "<<largest<<" smallest "<<smallest<<std::endl;
    258257    choleskyCondition_=largest/smallest;
    259     //drop fresh makes some formADAT easier
     258    /*drop fresh makes some formADAT easier*/
    260259    if (newDropped||numberRowsDropped_) {
    261260      newDropped=0;
     
    264263        rowsDropped_[i]=dropped;
    265264        if (dropped==2) {
    266           //dropped this time
     265          /*dropped this time*/
    267266          rowsDropped[newDropped++]=i;
    268267          rowsDropped_[i]=0;
     
    273272    }
    274273  } else {
    275     // KKT
     274    /* KKT*/
    276275    CoinPackedMatrix * quadratic = NULL;
    277276    ClpQuadraticObjective * quadraticObj =
     
    279278    if (quadraticObj)
    280279      quadratic = quadraticObj->quadraticObjective();
    281     // matrix
     280    /* matrix*/
    282281    int numberRowsModel = model_->numberRows();
    283282    int numberColumns = model_->numberColumns();
    284283    int numberTotal = numberColumns + numberRowsModel;
    285284    longDouble * work = sparseFactor_;
    286     work--; // skip diagonal
     285    work--; /* skip diagonal*/
    287286    int addOffset=numberRows_-1;
    288287    int iColumn;
    289288    if (!quadratic) {
    290289      for (iColumn=0;iColumn<numberColumns;iColumn++) {
    291         longDouble value = diagonal[iColumn];
    292         if (fabs(value)>1.0e-100) {
     290        CoinWorkDouble value = diagonal[iColumn];
     291        if (CoinAbs(value)>1.0e-100) {
    293292          value = 1.0/value;
    294           largest = CoinMax(largest,fabs(value));
     293          largest = CoinMax(largest,CoinAbs(value));
    295294          diagonal_[iColumn] = -value;
    296295          CoinBigIndex start=columnStart[iColumn];
    297296          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
    298297          for (CoinBigIndex j=start;j<end;j++) {
    299             //choleskyRow_[numberElements]=row[j]+numberTotal;
    300             //sparseFactor_[numberElements++]=element[j];
     298            /*choleskyRow_[numberElements]=row[j]+numberTotal;*/
     299            /*sparseFactor_[numberElements++]=element[j];*/
    301300            work[row[j]+numberTotal]=element[j];
    302             largest = CoinMax(largest,fabs(element[j]));
     301            largest = CoinMax(largest,CoinAbs(element[j]));
    303302          }
    304303        } else {
     
    309308      }
    310309    } else {
    311       // Quadratic
     310      /* Quadratic*/
    312311      const int * columnQuadratic = quadratic->getIndices();
    313312      const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
     
    315314      const double * quadraticElement = quadratic->getElements();
    316315      for (iColumn=0;iColumn<numberColumns;iColumn++) {
    317         longDouble value = diagonal[iColumn];
     316        CoinWorkDouble value = diagonal[iColumn];
    318317        CoinBigIndex j;
    319         if (fabs(value)>1.0e-100) {
     318        if (CoinAbs(value)>1.0e-100) {
    320319          value = 1.0/value;
    321320          for (j=columnQuadraticStart[iColumn];
     
    328327            }
    329328          }
    330           largest = CoinMax(largest,fabs(value));
     329          largest = CoinMax(largest,CoinAbs(value));
    331330          diagonal_[iColumn] = -value;
    332331          CoinBigIndex start=columnStart[iColumn];
     
    334333          for (j=start;j<end;j++) {
    335334            work[row[j]+numberTotal]=element[j];
    336             largest = CoinMax(largest,fabs(element[j]));
     335            largest = CoinMax(largest,CoinAbs(element[j]));
    337336          }
    338337        } else {
     
    344343      }
    345344    }
    346     // slacks
     345    /* slacks*/
    347346    for (iColumn=numberColumns;iColumn<numberTotal;iColumn++) {
    348       longDouble value = diagonal[iColumn];
    349       if (fabs(value)>1.0e-100) {
     347      CoinWorkDouble value = diagonal[iColumn];
     348      if (CoinAbs(value)>1.0e-100) {
    350349        value = 1.0/value;
    351         largest = CoinMax(largest,fabs(value));
     350        largest = CoinMax(largest,CoinAbs(value));
    352351      } else {
    353352        value = 1.0e100;
     
    358357      work += addOffset;
    359358    }
    360     // Finish diagonal
     359    /* Finish diagonal*/
    361360    for (iRow=0;iRow<numberRowsModel;iRow++) {
    362361      diagonal_[iRow+numberTotal]=delta2;
    363362    }
    364     //check sizes
     363    /*check sizes*/
    365364    largest*=1.0e-20;
    366     largest = CoinMin(largest,1.0e-11);
     365    largest = CoinMin(largest,CHOL_SMALL_VALUE);
    367366    doubleParameters_[10]=CoinMax(1.0e-20,largest);
    368367    integerParameters_[20]=0;
    369368    doubleParameters_[3]=0.0;
    370369    doubleParameters_[4]=COIN_DBL_MAX;
    371     // Set up LDL cutoff
     370    /* Set up LDL cutoff*/
    372371    integerParameters_[34]=numberTotal;
    373     // KKT
     372    /* KKT*/
    374373    int * rowsDropped2 = new int[numberRows_];
    375374    CoinZeroN(rowsDropped2,numberRows_);
     
    385384      std::cout<<"Cholesky - largest "<<largest<<" smallest "<<smallest<<std::endl;
    386385    choleskyCondition_=largest/smallest;
    387     // Should save adjustments in ..R_
     386    /* Should save adjustments in ..R_*/
    388387    int n1=0,n2=0;
    389     double * primalR = model_->primalR();
    390     double * dualR = model_->dualR();
     388    CoinWorkDouble * primalR = model_->primalR();
     389    CoinWorkDouble * dualR = model_->dualR();
    391390    for (iRow=0;iRow<numberTotal;iRow++) {
    392391      if (rowsDropped2[iRow]) {
    393392        n1++;
    394         //printf("row region1 %d dropped\n",iRow);
    395         //rowsDropped_[iRow]=1;
     393        /*printf("row region1 %d dropped\n",iRow);*/
     394        /*rowsDropped_[iRow]=1;*/
    396395        rowsDropped_[iRow]=0;
    397396        primalR[iRow]=doubleParameters_[20];
     
    404403      if (rowsDropped2[iRow]) {
    405404        n2++;
    406         //printf("row region2 %d dropped\n",iRow);
    407         //rowsDropped_[iRow]=1;
     405        /*printf("row region2 %d dropped\n",iRow);*/
     406        /*rowsDropped_[iRow]=1;*/
    408407        rowsDropped_[iRow]=0;
    409408        dualR[iRow-numberTotal]=doubleParameters_[34];
     
    425424  diagonal_ = sparseFactor_ + 40000;
    426425  sparseFactor_=diagonal_ + numberRows_;
    427   //memcpy(sparseFactor_,xx,sizeFactor_*sizeof(double));
     426  /*memcpy(sparseFactor_,xx,sizeFactor_*sizeof(double));*/
    428427  CoinMemcpyN(xx,40000,sparseFactor_);
    429428  CoinMemcpyN(yy,numberRows_,diagonal_);
    430429  int numberDropped=0;
    431   double largest=0.0;
    432   double smallest=COIN_DBL_MAX;
     430  CoinWorkDouble largest=0.0;
     431  CoinWorkDouble smallest=COIN_DBL_MAX;
    433432  double dropValue = doubleParameters_[10];
    434433  int firstPositive=integerParameters_[34];
    435434  longDouble * work = sparseFactor_;
    436   // Allow for triangular
     435  /* Allow for triangular*/
    437436  int addOffset=numberRows_-1;
    438437  work--;
     
    441440    int addOffsetNow = numberRows_-1;;
    442441    longDouble * workNow=sparseFactor_-1+iColumn;
    443     double diagonalValue = diagonal_[iColumn];
     442    CoinWorkDouble diagonalValue = diagonal_[iColumn];
    444443    for (iRow=0;iRow<iColumn;iRow++) {
    445444      double aj = *workNow;
    446445      addOffsetNow--;
    447446      workNow += addOffsetNow;
    448       double multiplier = workDouble_[iRow];
    449       diagonalValue -= aj*aj*multiplier;
     447      diagonalValue -= aj*aj*workDouble_[iRow];
    450448    }
    451449    bool dropColumn=false;
    452450    if (iColumn<firstPositive) {
    453       // must be negative
     451      /* must be negative*/
    454452      if (diagonalValue<=-dropValue) {
    455453        smallest = CoinMin(smallest,-diagonalValue);
     
    464462      }
    465463    } else {
    466       // must be positive
     464      /* must be positive*/
    467465      if (diagonalValue>=dropValue) {
    468466        smallest = CoinMin(smallest,diagonalValue);
     
    494492      }
    495493    } else {
    496       // drop column
     494      /* drop column*/
    497495      rowsDropped[iColumn]=2;
    498496      numberDropped++;
     
    511509  diagonal_=yy;
    512510}
    513 //#define POS_DEBUG
     511/*#define POS_DEBUG*/
    514512#ifdef POS_DEBUG
    515513static int counter=0;
     
    522520  iCol=0;
    523521  int kk=k>>BLOCKSQSHIFT;
    524   //printf("%d %d %d %d\n",k,kk,BLOCKSQ,BLOCKSQSHIFT);
     522  /*printf("%d %d %d %d\n",k,kk,BLOCKSQ,BLOCKSQSHIFT);*/
    525523  k=kk;
    526524  while (k>=numberBlocks) {
     
    541539{
    542540  int iColumn;
    543   //longDouble * xx = sparseFactor_;
    544   //longDouble * yy = diagonal_;
    545   //diagonal_ = sparseFactor_ + 40000;
    546   //sparseFactor_=diagonal_ + numberRows_;
    547   //memcpy(sparseFactor_,xx,sizeFactor_*sizeof(double));
    548   //memcpy(sparseFactor_,xx,40000*sizeof(double));
    549   //memcpy(diagonal_,yy,numberRows_*sizeof(double));
     541  /*longDouble * xx = sparseFactor_;*/
     542  /*longDouble * yy = diagonal_;*/
     543  /*diagonal_ = sparseFactor_ + 40000;*/
     544  /*sparseFactor_=diagonal_ + numberRows_;*/
     545  /*memcpy(sparseFactor_,xx,sizeFactor_*sizeof(double));*/
     546  /*memcpy(sparseFactor_,xx,40000*sizeof(double));*/
     547  /*memcpy(diagonal_,yy,numberRows_*sizeof(double));*/
    550548  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
    551   // later align on boundary
     549  /* later align on boundary*/
    552550  longDouble * a = sparseFactor_+BLOCKSQ*numberBlocks;
    553551  int n=numberRows_;
    554552  int nRound = numberRows_&(~(BLOCK-1));
    555   // adjust if exact
     553  /* adjust if exact*/
    556554  if (nRound==n)
    557555    nRound -= BLOCK;
    558556  int sizeLastBlock=n-nRound;
    559   int get = n*(n-1)/2; // ? as no diagonal
     557  int get = n*(n-1)/2; /* ? as no diagonal*/
    560558  int block = numberBlocks*(numberBlocks+1)/2;
    561559  int ifOdd;
     
    566564    ifOdd=1;
    567565    int put = BLOCKSQ;
    568     // do last separately
     566    /* do last separately*/
    569567    put -= (BLOCK-sizeLastBlock)*(BLOCK+1);
    570568    for (iColumn=numberRows_-1;iColumn>=nRound;iColumn--) {
     
    575573        assert (aa+put2>=sparseFactor_+get);
    576574      }
    577       // save diagonal as well
     575      /* save diagonal as well*/
    578576      aa[--put2]=diagonal_[iColumn];
    579577    }
     
    581579    block--;
    582580  } else {
    583     // exact fit
     581    /* exact fit*/
    584582    rowLast=numberRows_-1;
    585583    ifOdd=0;
    586584  }
    587   // Now main loop
     585  /* Now main loop*/
    588586  int nBlock=0;
    589587  for (;n>0;n-=BLOCK) {
     
    592590    int put = BLOCKSQ;
    593591    int putLast=0;
    594     // see if we have small block
     592    /* see if we have small block*/
    595593    if (ifOdd) {
    596594      aaLast = &a[(block-1)*BLOCKSQ];
     
    600598    for (iColumn=n-1;iColumn>=n-BLOCK;iColumn--) {
    601599      if (aaLast) {
    602         // last bit
     600        /* last bit*/
    603601        for (int iRow=numberRows_-1;iRow>rowLast;iRow--) {
    604602          aaLast[--putLast] = sparseFactor_[--get];
     
    617615        }
    618616        if (j-BLOCK<iColumn) {
    619           // save diagonal as well
     617          /* save diagonal as well*/
    620618          aPut[--put2]=diagonal_[iColumn];
    621619        }
     
    628626    block -= nBlock+ifOdd;
    629627  }
    630   factor(a,numberRows_,numberBlocks,
     628  ClpCholeskyDenseC  info;
     629  info.diagonal_=diagonal_;
     630  info.doubleParameters_[0]=doubleParameters_[10];
     631  info.integerParameters_[0]=integerParameters_[34];
     632#ifndef CLP_CILK
     633  ClpCholeskyCfactor(&info,a,numberRows_,numberBlocks,
    631634         diagonal_,workDouble_,rowsDropped);
    632   //sparseFactor_=xx;
    633   //diagonal_=yy;
    634 }
    635 // Non leaf recursive factor
    636 void
    637 ClpCholeskyDense::factor(longDouble * a, int n, int numberBlocks,
     635#else
     636  info.a=a;
     637  info.n=numberRows_;
     638  info.numberBlocks=numberBlocks;
     639  info.work=workDouble_;
     640  info.rowsDropped=rowsDropped;
     641  info.integerParameters_[1]=model_->numberThreads();
     642  ClpCholeskySpawn(&info);
     643#endif
     644  double largest=0.0;
     645  double smallest=COIN_DBL_MAX;
     646  int numberDropped=0;
     647  for (int i=0;i<numberRows_;i++) {
     648    if (diagonal_[i]) {
     649      largest = CoinMax(largest,CoinAbs(diagonal_[i]));
     650      smallest = CoinMin(smallest,CoinAbs(diagonal_[i]));
     651    } else {
     652      numberDropped++;
     653    }
     654  }
     655  doubleParameters_[3]=CoinMax(doubleParameters_[3],1.0/smallest);
     656  doubleParameters_[4]=CoinMin(doubleParameters_[4],1.0/largest);
     657  integerParameters_[20]+= numberDropped;
     658}
     659  /* Non leaf recursive factor*/
     660void 
     661ClpCholeskyCfactor(ClpCholeskyDenseC * thisStruct, longDouble * a, int n, int numberBlocks,
    638662            longDouble * diagonal, longDouble * work, int * rowsDropped)
    639663{
    640664  if (n<=BLOCK) {
    641     factorLeaf(a,n,diagonal,work,rowsDropped);
     665    ClpCholeskyCfactorLeaf(thisStruct, a,n,diagonal,work,rowsDropped);
    642666  } else {
    643667    int nb=number_blocks((n+1)>>1);
     
    647671    int nintri=(nb*(nb+1))>>1;
    648672    int nbelow=(numberBlocks-nb)*nb;
    649     factor(a,nThis,numberBlocks,diagonal,work,rowsDropped);
    650     triRec(a,nThis,a+number_entries(nb),diagonal,work,nLeft,nb,0,numberBlocks);
     673    ClpCholeskyCfactor(thisStruct, a,nThis,numberBlocks,diagonal,work,rowsDropped);
     674    ClpCholeskyCtriRec(thisStruct, a,nThis,a+number_entries(nb),diagonal,work,nLeft,nb,0,numberBlocks);
    651675    aother=a+number_entries(nintri+nbelow);
    652     recTri(a+number_entries(nb),nLeft,nThis,nb,0,aother,diagonal,work,numberBlocks);
    653     factor(aother,nLeft,
     676    ClpCholeskyCrecTri(thisStruct, a+number_entries(nb),nLeft,nThis,nb,0,aother,diagonal,work,numberBlocks);
     677    ClpCholeskyCfactor(thisStruct, aother,nLeft,
    654678        numberBlocks-nb,diagonal+nThis,work+nThis,rowsDropped);
    655679  }
    656680}
    657 // Non leaf recursive triangle rectangle update
     681  /* Non leaf recursive triangle rectangle update*/
    658682void
    659 ClpCholeskyDense::triRec(longDouble * aTri, int nThis, longDouble * aUnder,
     683ClpCholeskyCtriRec(ClpCholeskyDenseC * thisStruct, longDouble * aTri, int nThis, longDouble * aUnder,
    660684                         longDouble * diagonal, longDouble * work,
    661685                         int nLeft, int iBlock, int jBlock,
     
    663687{
    664688  if (nThis<=BLOCK&&nLeft<=BLOCK) {
    665     triRecLeaf(aTri,aUnder,diagonal,work,nLeft);
     689    ClpCholeskyCtriRecLeaf(/*thisStruct,*/ aTri,aUnder,diagonal,work,nLeft);
    666690  } else if (nThis<nLeft) {
    667691    int nb=number_blocks((nLeft+1)>>1);
    668692    int nLeft2=number_rows(nb);
    669     triRec(aTri,nThis,aUnder,diagonal,work,nLeft2,iBlock,jBlock,numberBlocks);
    670     triRec(aTri,nThis,aUnder+number_entries(nb),diagonal,work,nLeft-nLeft2,
     693    ClpCholeskyCtriRec(thisStruct, aTri,nThis,aUnder,diagonal,work,nLeft2,iBlock,jBlock,numberBlocks);
     694    ClpCholeskyCtriRec(thisStruct, aTri,nThis,aUnder+number_entries(nb),diagonal,work,nLeft-nLeft2,
    671695          iBlock+nb,jBlock,numberBlocks);
    672696  } else {
     
    678702    int nintri=(nb*(nb+1))>>1;
    679703    int nbelow=(numberBlocks-nb)*nb;
    680     triRec(aTri,nThis2,aUnder,diagonal,work,nLeft,iBlock,jBlock,numberBlocks);
     704    ClpCholeskyCtriRec(thisStruct, aTri,nThis2,aUnder,diagonal,work,nLeft,iBlock,jBlock,numberBlocks);
    681705    /* and rectangular update */
    682706    i=((numberBlocks-jBlock)*(numberBlocks-jBlock-1)-
    683707      (numberBlocks-jBlock-nb)*(numberBlocks-jBlock-nb-1))>>1;
    684708    aother=aUnder+number_entries(i);
    685     recRec(aTri+number_entries(nb),nThis-nThis2,nLeft,nThis2,aUnder,aother,
     709    ClpCholeskyCrecRec(thisStruct, aTri+number_entries(nb),nThis-nThis2,nLeft,nThis2,aUnder,aother,
    686710          work,kBlock,jBlock,numberBlocks);
    687     triRec(aTri+number_entries(nintri+nbelow),nThis-nThis2,aother,diagonal+nThis2,
     711    ClpCholeskyCtriRec(thisStruct, aTri+number_entries(nintri+nbelow),nThis-nThis2,aother,diagonal+nThis2,
    688712           work+nThis2,nLeft,
    689713      iBlock-nb,kBlock-nb,numberBlocks-nb);
    690714  }
    691715}
    692 // Non leaf recursive rectangle triangle update
     716  /* Non leaf recursive rectangle triangle update*/
    693717void
    694 ClpCholeskyDense::recTri(longDouble * aUnder, int nTri, int nDo,
     718ClpCholeskyCrecTri(ClpCholeskyDenseC * thisStruct, longDouble * aUnder, int nTri, int nDo,
    695719                         int iBlock, int jBlock,longDouble * aTri,
    696720                         longDouble * diagonal, longDouble * work,
     
    698722{
    699723  if (nTri<=BLOCK&&nDo<=BLOCK) {
    700     recTriLeaf(aUnder,aTri,diagonal,work,nTri);
     724    ClpCholeskyCrecTriLeaf(/*thisStruct,*/ aUnder,aTri,/*diagonal,*/work,nTri);
    701725  } else if (nTri<nDo) {
    702726    int nb=number_blocks((nDo+1)>>1);
     
    704728    longDouble * aother;
    705729    int i;
    706     recTri(aUnder,nTri,nDo2,iBlock,jBlock,aTri,diagonal,work,numberBlocks);
     730    ClpCholeskyCrecTri(thisStruct, aUnder,nTri,nDo2,iBlock,jBlock,aTri,diagonal,work,numberBlocks);
    707731    i=((numberBlocks-jBlock)*(numberBlocks-jBlock-1)-
    708732      (numberBlocks-jBlock-nb)*(numberBlocks-jBlock-nb-1))>>1;
    709733    aother=aUnder+number_entries(i);
    710     recTri(aother,nTri,nDo-nDo2,iBlock-nb,jBlock,aTri,diagonal+nDo2,
     734    ClpCholeskyCrecTri(thisStruct, aother,nTri,nDo-nDo2,iBlock-nb,jBlock,aTri,diagonal+nDo2,
    711735           work+nDo2,numberBlocks-nb);
    712736  } else {
     
    715739    longDouble * aother;
    716740    int i;
    717     recTri(aUnder,nTri2,nDo,iBlock,jBlock,aTri,diagonal,work,numberBlocks);
     741    ClpCholeskyCrecTri(thisStruct, aUnder,nTri2,nDo,iBlock,jBlock,aTri,diagonal,work,numberBlocks);
    718742    /* and rectangular update */
    719743    i=((numberBlocks-iBlock)*(numberBlocks-iBlock+1)-
    720744      (numberBlocks-iBlock-nb)*(numberBlocks-iBlock-nb+1))>>1;
    721745    aother=aTri+number_entries(nb);
    722     recRec(aUnder,nTri2,nTri-nTri2,nDo,aUnder+number_entries(nb),aother,
     746    ClpCholeskyCrecRec(thisStruct, aUnder,nTri2,nTri-nTri2,nDo,aUnder+number_entries(nb),aother,
    723747           work,iBlock,jBlock,numberBlocks);
    724     recTri(aUnder+number_entries(nb),nTri-nTri2,nDo,iBlock+nb,jBlock,
     748    ClpCholeskyCrecTri(thisStruct, aUnder+number_entries(nb),nTri-nTri2,nDo,iBlock+nb,jBlock,
    725749         aTri+number_entries(i),diagonal,work,numberBlocks);
    726750  }
     
    731755*/
    732756void
    733 ClpCholeskyDense::recRec(longDouble * above, int nUnder, int nUnderK,
     757ClpCholeskyCrecRec(ClpCholeskyDenseC * thisStruct, longDouble * above, int nUnder, int nUnderK,
    734758                         int nDo, longDouble * aUnder, longDouble *aOther,
    735759                         longDouble * work,
     
    739763  if (nDo<=BLOCK&&nUnder<=BLOCK&&nUnderK<=BLOCK) {
    740764    assert (nDo==BLOCK&&nUnder==BLOCK);
    741     recRecLeaf(above , aUnder ,  aOther, work, nUnderK);
     765    ClpCholeskyCrecRecLeaf(/*thisStruct,*/ above , aUnder ,  aOther, work, nUnderK);
    742766  } else if (nDo<=nUnderK&&nUnder<=nUnderK) {
    743767    int nb=number_blocks((nUnderK+1)>>1);
    744768    int nUnder2=number_rows(nb);
    745     recRec(above,nUnder,nUnder2,nDo,aUnder,aOther,work,
     769    ClpCholeskyCrecRec(thisStruct, above,nUnder,nUnder2,nDo,aUnder,aOther,work,
    746770               iBlock,jBlock,numberBlocks);
    747     recRec(above,nUnder,nUnderK-nUnder2,nDo,aUnder+number_entries(nb),
     771    ClpCholeskyCrecRec(thisStruct, above,nUnder,nUnderK-nUnder2,nDo,aUnder+number_entries(nb),
    748772        aOther+number_entries(nb),work,iBlock,jBlock,numberBlocks);
    749773  } else if (nUnderK<=nDo&&nUnder<=nDo) {
     
    751775    int nDo2=number_rows(nb);
    752776    int i;
    753     recRec(above,nUnder,nUnderK,nDo2,aUnder,aOther,work,
     777    ClpCholeskyCrecRec(thisStruct, above,nUnder,nUnderK,nDo2,aUnder,aOther,work,
    754778         iBlock,jBlock,numberBlocks);
    755779    i=((numberBlocks-jBlock)*(numberBlocks-jBlock-1)-
    756780      (numberBlocks-jBlock-nb)*(numberBlocks-jBlock-nb-1))>>1;
    757     recRec(above+number_entries(i),nUnder,nUnderK,nDo-nDo2,
     781    ClpCholeskyCrecRec(thisStruct, above+number_entries(i),nUnder,nUnderK,nDo-nDo2,
    758782           aUnder+number_entries(i),
    759783           aOther,work+nDo2,
     
    763787    int nUnder2=number_rows(nb);
    764788    int i;
    765     recRec(above,nUnder2,nUnderK,nDo,aUnder,aOther,work,
     789    ClpCholeskyCrecRec(thisStruct, above,nUnder2,nUnderK,nDo,aUnder,aOther,work,
    766790               iBlock,jBlock,numberBlocks);
    767791    i=((numberBlocks-iBlock)*(numberBlocks-iBlock-1)-
    768792      (numberBlocks-iBlock-nb)*(numberBlocks-iBlock-nb-1))>>1;
    769     recRec(above+number_entries(nb),nUnder-nUnder2,nUnderK,nDo,aUnder,
     793    ClpCholeskyCrecRec(thisStruct, above+number_entries(nb),nUnder-nUnder2,nUnderK,nDo,aUnder,
    770794        aOther+number_entries(i),work,iBlock+nb,jBlock,numberBlocks);
    771795  }
    772796}
    773 // Leaf recursive factor
     797  /* Leaf recursive factor*/
    774798void
    775 ClpCholeskyDense::factorLeaf(longDouble * a, int n,
     799ClpCholeskyCfactorLeaf(ClpCholeskyDenseC * thisStruct, longDouble * a, int n,
    776800                longDouble * diagonal, longDouble * work, int * rowsDropped)
    777801{
    778   longDouble largest=doubleParameters_[3];
    779   longDouble smallest=doubleParameters_[4];
    780   double dropValue = doubleParameters_[10];
    781   int firstPositive=integerParameters_[34];
    782   size_t rowOffset=diagonal-diagonal_;
    783   int numberDropped=0;
     802  double dropValue = thisStruct->doubleParameters_[0];
     803  int firstPositive=thisStruct->integerParameters_[0];
     804  int rowOffset=diagonal-thisStruct->diagonal_;
    784805  int i, j, k;
    785   longDouble t00,temp1;
     806  CoinWorkDouble t00,temp1;
    786807  longDouble * aa;
    787808  aa = a-BLOCK;
    788809  for (j = 0; j < n; j ++) {
     810    bool dropColumn;
     811    CoinWorkDouble useT00;
    789812    aa+=BLOCK;
    790813    t00 = aa[j];
    791814    for (k = 0; k < j; ++k) {
    792       longDouble multiplier = work[k];
     815      CoinWorkDouble multiplier = work[k];
    793816      t00 -= a[j + k * BLOCK] * a[j + k * BLOCK] * multiplier;
    794817    }
    795     bool dropColumn=false;
    796     longDouble useT00=t00;
    797     if ((int)(j+rowOffset)<firstPositive) {
    798       // must be negative
     818    dropColumn=false;
     819    useT00=t00;
     820    if (j+rowOffset<firstPositive) {
     821      /* must be negative*/
    799822      if (t00<=-dropValue) {
    800         smallest = CoinMin(smallest,-t00);
    801         largest = CoinMax(largest,-t00);
    802         //aa[j]=t00;
     823        /*aa[j]=t00;*/
    803824        t00 = 1.0/t00;
    804825      } else {
    805826        dropColumn=true;
    806         //aa[j]=-1.0e100;
     827        /*aa[j]=-1.0e100;*/
    807828        useT00=-1.0e-100;
    808829        t00=0.0;
    809         integerParameters_[20]++;
    810830      }
    811831    } else {
    812       // must be positive
     832      /* must be positive*/
    813833      if (t00>=dropValue) {
    814         smallest = CoinMin(smallest,t00);
    815         largest = CoinMax(largest,t00);
    816         //aa[j]=t00;
     834        /*aa[j]=t00;*/
    817835        t00 = 1.0/t00;
    818836      } else {
    819837        dropColumn=true;
    820         //aa[j]=1.0e100;
     838        /*aa[j]=1.0e100;*/
    821839        useT00=1.0e-100;
    822840        t00=0.0;
    823         integerParameters_[20]++;
    824841      }
    825842    }
     
    831848        t00=aa[i];
    832849        for (k = 0; k < j; ++k) {
    833           longDouble multiplier = work[k];
     850          CoinWorkDouble multiplier = work[k];
    834851          t00 -= a[i + k * BLOCK] * a[j + k * BLOCK] * multiplier;
    835852        }
     
    837854      }
    838855    } else {
    839       // drop column
     856      /* drop column*/
    840857      rowsDropped[j+rowOffset]=2;
    841       numberDropped++;
    842858      diagonal[j]=0.0;
    843       //aa[j]=1.0e100;
     859      /*aa[j]=1.0e100;*/
    844860      work[j]=1.0e100;
    845861      for (i=j+1;i<n;i++) {
     
    848864    }
    849865  }
    850   doubleParameters_[3]=largest;
    851   doubleParameters_[4]=smallest;
    852   integerParameters_[20]+=numberDropped;
    853 }
    854 // Leaf recursive triangle rectangle update
     866}
     867  /* Leaf recursive triangle rectangle update*/
    855868void
    856 ClpCholeskyDense::triRecLeaf(longDouble * aTri, longDouble * aUnder, longDouble * diagonal, longDouble * work,
     869ClpCholeskyCtriRecLeaf(/*ClpCholeskyDenseC * thisStruct,*/ longDouble * aTri, longDouble * aUnder, longDouble * diagonal, longDouble * work,
    857870                int nUnder)
    858871{
     
    862875  int irt,ict;
    863876  int it=bNumber(aTri,irt,ict);
    864   //printf("%d %d\n",iu,it);
     877  /*printf("%d %d\n",iu,it);*/
    865878  printf("trirecleaf  under (%d,%d), tri (%d,%d)\n",
    866879         iru,icu,irt,ict);
    867   assert (diagonal==diagonal_+ict*BLOCK);
     880  assert (diagonal==thisStruct->diagonal_+ict*BLOCK);
    868881#endif
    869882  int j;
     
    873886    aa = aTri-2*BLOCK;
    874887    for (j = 0; j < BLOCK; j +=2) {
     888      int i;
     889      CoinWorkDouble temp0 = diagonal[j];
     890      CoinWorkDouble temp1 = diagonal[j+1];
    875891      aa+=2*BLOCK;
    876       longDouble temp0 = diagonal[j];
    877       longDouble temp1 = diagonal[j+1];
    878       for (int i=0;i<BLOCK;i+=2) {
    879         longDouble t00=aUnder[i+j*BLOCK];
    880         longDouble t10=aUnder[i+ BLOCK + j*BLOCK];
    881         longDouble t01=aUnder[i+1+j*BLOCK];
    882         longDouble t11=aUnder[i+1+ BLOCK + j*BLOCK];
    883         for (int k = 0; k < j; ++k) {
    884           longDouble multiplier=work[k];
    885           longDouble au0 = aUnder[i + k * BLOCK]*multiplier;
    886           longDouble au1 = aUnder[i + 1 + k * BLOCK]*multiplier;
    887           longDouble at0 = aTri[j + k * BLOCK];
    888           longDouble at1 = aTri[j + 1 + k * BLOCK];
     892      for ( i=0;i<BLOCK;i+=2) {
     893        CoinWorkDouble at1;
     894        CoinWorkDouble t00=aUnder[i+j*BLOCK];
     895        CoinWorkDouble t10=aUnder[i+ BLOCK + j*BLOCK];
     896        CoinWorkDouble t01=aUnder[i+1+j*BLOCK];
     897        CoinWorkDouble t11=aUnder[i+1+ BLOCK + j*BLOCK];
     898        int k;
     899        for (k = 0; k < j; ++k) {
     900          CoinWorkDouble multiplier=work[k];
     901          CoinWorkDouble au0 = aUnder[i + k * BLOCK]*multiplier;
     902          CoinWorkDouble au1 = aUnder[i + 1 + k * BLOCK]*multiplier;
     903          CoinWorkDouble at0 = aTri[j + k * BLOCK];
     904          CoinWorkDouble at1 = aTri[j + 1 + k * BLOCK];
    889905          t00 -= au0 * at0;
    890906          t10 -= au0 * at1;
     
    893909        }
    894910        t00 *= temp0;
    895         longDouble at1 = aTri[j + 1 + j*BLOCK]*work[j];
     911        at1 = aTri[j + 1 + j*BLOCK]*work[j];
    896912        t10 -= t00 * at1;
    897913        t01 *= temp0;
     
    907923    aa = aTri-BLOCK;
    908924    for (j = 0; j < BLOCK; j ++) {
     925      int i;
     926      CoinWorkDouble temp1 = diagonal[j];
    909927      aa+=BLOCK;
    910       longDouble temp1 = diagonal[j];
    911       for (int i=0;i<nUnder;i++) {
    912         longDouble t00=aUnder[i+j*BLOCK];
    913         for (int k = 0; k < j; ++k) {
    914           longDouble multiplier=work[k];
     928      for (i=0;i<nUnder;i++) {
     929        int k;
     930        CoinWorkDouble t00=aUnder[i+j*BLOCK];
     931        for ( k = 0; k < j; ++k) {
     932          CoinWorkDouble multiplier=work[k];
    915933          t00 -= aUnder[i + k * BLOCK] * aTri[j + k * BLOCK] * multiplier;
    916934        }
     
    922940#endif
    923941}
    924 // Leaf recursive rectangle triangle update
    925 void ClpCholeskyDense::recTriLeaf(longDouble * aUnder, longDouble * aTri,
    926                                   longDouble * diagonal, longDouble * work, int nUnder)
     942  /* Leaf recursive rectangle triangle update*/
     943void ClpCholeskyCrecTriLeaf(/*ClpCholeskyDenseC * thisStruct,*/ longDouble * aUnder, longDouble * aTri,
     944                            /*longDouble * diagonal,*/ longDouble * work, int nUnder)
    927945{
    928946#ifdef POS_DEBUG
     
    931949  int irt,ict;
    932950  int it=bNumber(aTri,irt,ict);
    933   //printf("%d %d\n",iu,it);
     951  /*printf("%d %d\n",iu,it);*/
    934952  printf("rectrileaf  under (%d,%d), tri (%d,%d)\n",
    935953         iru,icu,irt,ict);
    936   assert (diagonal==diagonal_+icu*BLOCK);
     954  assert (diagonal==thisStruct->diagonal_+icu*BLOCK);
    937955#endif
    938956  int i, j, k;
    939   longDouble t00;
     957  CoinWorkDouble t00;
    940958  longDouble * aa;
    941959#ifdef BLOCKUNROLL
     
    944962    aa = aTri-2*BLOCK;
    945963    for (j = 0; j < BLOCK; j +=2) {
    946       longDouble t00,t01,t10,t11;
     964      CoinWorkDouble t00,t01,t10,t11;
    947965      aa+=2*BLOCK;
    948966      aUnder2+=2;
     
    951969      t10=aa[j+1+BLOCK];
    952970      for (k = 0; k < BLOCK; ++k) {
    953         longDouble multiplier = work[k];
    954         longDouble a0 = aUnder2[k * BLOCK];
    955         longDouble a1 = aUnder2[1 + k * BLOCK];
    956         longDouble x0 = a0 * multiplier;
    957         longDouble x1 = a1 * multiplier;
     971        CoinWorkDouble multiplier = work[k];
     972        CoinWorkDouble a0 = aUnder2[k * BLOCK];
     973        CoinWorkDouble a1 = aUnder2[1 + k * BLOCK];
     974        CoinWorkDouble x0 = a0 * multiplier;
     975        CoinWorkDouble x1 = a1 * multiplier;
    958976        t00 -= a0 * x0;
    959977        t01 -= a1 * x0;
     
    969987        t11=aa[i+1+BLOCK];
    970988        for (k = 0; k < BLOCK; ++k) {
    971           longDouble multiplier = work[k];
    972           longDouble a0 = aUnder2[k * BLOCK]*multiplier;
    973           longDouble a1 = aUnder2[1 + k * BLOCK]*multiplier;
     989          CoinWorkDouble multiplier = work[k];
     990          CoinWorkDouble a0 = aUnder2[k * BLOCK]*multiplier;
     991          CoinWorkDouble a1 = aUnder2[1 + k * BLOCK]*multiplier;
    974992          t00 -= aUnder[i + k * BLOCK] * a0;
    975993          t01 -= aUnder[i + k * BLOCK] * a1;
     
    9911009        t00=aa[i];
    9921010        for (k = 0; k < BLOCK; ++k) {
    993           longDouble multiplier = work[k];
     1011          CoinWorkDouble multiplier = work[k];
    9941012          t00 -= aUnder[i + k * BLOCK] * aUnder[j + k * BLOCK]*multiplier;
    9951013        }
     
    10061024*/
    10071025void
    1008 ClpCholeskyDense::recRecLeaf(const longDouble * COIN_RESTRICT above,
     1026ClpCholeskyCrecRecLeaf(/*ClpCholeskyDenseC * thisStruct,*/
     1027                                         const longDouble * COIN_RESTRICT above,
    10091028                             const longDouble * COIN_RESTRICT aUnder,
    10101029                             longDouble * COIN_RESTRICT aOther,
     
    10191038  int iro,ico;
    10201039  int io=bNumber(aOther,iro,ico);
    1021   //printf("%d %d %d\n",ia,iu,io);
     1040  /*printf("%d %d %d\n",ia,iu,io);*/
    10221041  printf("recrecleaf above (%d,%d), under (%d,%d), other (%d,%d)\n",
    10231042         ira,ica,iru,icu,iro,ico);
     
    10281047  aa = aOther-4*BLOCK;
    10291048  if (nUnder==BLOCK) {
    1030     //#define INTEL
     1049    /*#define INTEL*/
    10311050#ifdef INTEL
    10321051    aa+=2*BLOCK;
     
    10341053      aa+=2*BLOCK;
    10351054      for (i=0;i< BLOCK;i+=2) {
    1036         longDouble t00=aa[i+0*BLOCK];
    1037         longDouble t10=aa[i+1*BLOCK];
    1038         longDouble t01=aa[i+1+0*BLOCK];
    1039         longDouble t11=aa[i+1+1*BLOCK];
     1055        CoinWorkDouble t00=aa[i+0*BLOCK];
     1056        CoinWorkDouble t10=aa[i+1*BLOCK];
     1057        CoinWorkDouble t01=aa[i+1+0*BLOCK];
     1058        CoinWorkDouble t11=aa[i+1+1*BLOCK];
    10401059        for (k=0;k<BLOCK;k++) {
    1041           longDouble multiplier = work[k];
    1042           longDouble a00=aUnder[i+k*BLOCK]*multiplier;
    1043           longDouble a01=aUnder[i+1+k*BLOCK]*multiplier;
     1060          CoinWorkDouble multiplier = work[k];
     1061          CoinWorkDouble a00=aUnder[i+k*BLOCK]*multiplier;
     1062          CoinWorkDouble a01=aUnder[i+1+k*BLOCK]*multiplier;
    10441063          t00 -= a00 * above[j + 0 + k * BLOCK];
    10451064          t10 -= a00 * above[j + 1 + k * BLOCK];
     
    10571076      aa+=4*BLOCK;
    10581077      for (i=0;i< BLOCK;i+=4) {
    1059         longDouble t00=aa[i+0+0*BLOCK];
    1060         longDouble t10=aa[i+0+1*BLOCK];
    1061         longDouble t20=aa[i+0+2*BLOCK];
    1062         longDouble t30=aa[i+0+3*BLOCK];
    1063         longDouble t01=aa[i+1+0*BLOCK];
    1064         longDouble t11=aa[i+1+1*BLOCK];
    1065         longDouble t21=aa[i+1+2*BLOCK];
    1066         longDouble t31=aa[i+1+3*BLOCK];
    1067         longDouble t02=aa[i+2+0*BLOCK];
    1068         longDouble t12=aa[i+2+1*BLOCK];
    1069         longDouble t22=aa[i+2+2*BLOCK];
    1070         longDouble t32=aa[i+2+3*BLOCK];
    1071         longDouble t03=aa[i+3+0*BLOCK];
    1072         longDouble t13=aa[i+3+1*BLOCK];
    1073         longDouble t23=aa[i+3+2*BLOCK];
    1074         longDouble t33=aa[i+3+3*BLOCK];
     1078        CoinWorkDouble t00=aa[i+0+0*BLOCK];
     1079        CoinWorkDouble t10=aa[i+0+1*BLOCK];
     1080        CoinWorkDouble t20=aa[i+0+2*BLOCK];
     1081        CoinWorkDouble t30=aa[i+0+3*BLOCK];
     1082        CoinWorkDouble t01=aa[i+1+0*BLOCK];
     1083        CoinWorkDouble t11=aa[i+1+1*BLOCK];
     1084        CoinWorkDouble t21=aa[i+1+2*BLOCK];
     1085        CoinWorkDouble t31=aa[i+1+3*BLOCK];
     1086        CoinWorkDouble t02=aa[i+2+0*BLOCK];
     1087        CoinWorkDouble t12=aa[i+2+1*BLOCK];
     1088        CoinWorkDouble t22=aa[i+2+2*BLOCK];
     1089        CoinWorkDouble t32=aa[i+2+3*BLOCK];
     1090        CoinWorkDouble t03=aa[i+3+0*BLOCK];
     1091        CoinWorkDouble t13=aa[i+3+1*BLOCK];
     1092        CoinWorkDouble t23=aa[i+3+2*BLOCK];
     1093        CoinWorkDouble t33=aa[i+3+3*BLOCK];
    10751094        const longDouble * COIN_RESTRICT aUnderNow = aUnder+i;
    10761095        const longDouble * COIN_RESTRICT aboveNow = above+j;
    10771096        for (k=0;k<BLOCK;k++) {
    1078           longDouble multiplier = work[k];
    1079           longDouble a00=aUnderNow[0]*multiplier;
    1080           longDouble a01=aUnderNow[1]*multiplier;
    1081           longDouble a02=aUnderNow[2]*multiplier;
    1082           longDouble a03=aUnderNow[3]*multiplier;
     1097          CoinWorkDouble multiplier = work[k];
     1098          CoinWorkDouble a00=aUnderNow[0]*multiplier;
     1099          CoinWorkDouble a01=aUnderNow[1]*multiplier;
     1100          CoinWorkDouble a02=aUnderNow[2]*multiplier;
     1101          CoinWorkDouble a03=aUnderNow[3]*multiplier;
    10831102          t00 -= a00 * aboveNow[0];
    10841103          t10 -= a00 * aboveNow[1];
     
    11251144      aa+=4*BLOCK;
    11261145      for (i=0;i< n;i+=2) {
    1127         longDouble t00=aa[i+0*BLOCK];
    1128         longDouble t10=aa[i+1*BLOCK];
    1129         longDouble t20=aa[i+2*BLOCK];
    1130         longDouble t30=aa[i+3*BLOCK];
    1131         longDouble t01=aa[i+1+0*BLOCK];
    1132         longDouble t11=aa[i+1+1*BLOCK];
    1133         longDouble t21=aa[i+1+2*BLOCK];
    1134         longDouble t31=aa[i+1+3*BLOCK];
     1146        CoinWorkDouble t00=aa[i+0*BLOCK];
     1147        CoinWorkDouble t10=aa[i+1*BLOCK];
     1148        CoinWorkDouble t20=aa[i+2*BLOCK];
     1149        CoinWorkDouble t30=aa[i+3*BLOCK];
     1150        CoinWorkDouble t01=aa[i+1+0*BLOCK];
     1151        CoinWorkDouble t11=aa[i+1+1*BLOCK];
     1152        CoinWorkDouble t21=aa[i+1+2*BLOCK];
     1153        CoinWorkDouble t31=aa[i+1+3*BLOCK];
    11351154        const longDouble * COIN_RESTRICT aUnderNow = aUnder+i;
    11361155        const longDouble * COIN_RESTRICT aboveNow = above+j;
    11371156        for (k=0;k<BLOCK;k++) {
    1138           longDouble multiplier = work[k];
    1139           longDouble a00=aUnderNow[0]*multiplier;
    1140           longDouble a01=aUnderNow[1]*multiplier;
     1157          CoinWorkDouble multiplier = work[k];
     1158          CoinWorkDouble a00=aUnderNow[0]*multiplier;
     1159          CoinWorkDouble a01=aUnderNow[1]*multiplier;
    11411160          t00 -= a00 * aboveNow[0];
    11421161          t10 -= a00 * aboveNow[1];
     
    11601179      }
    11611180      if (odd) {
    1162         longDouble t0=aa[n+0*BLOCK];
    1163         longDouble t1=aa[n+1*BLOCK];
    1164         longDouble t2=aa[n+2*BLOCK];
    1165         longDouble t3=aa[n+3*BLOCK];
    1166         longDouble a0;
     1181        CoinWorkDouble t0=aa[n+0*BLOCK];
     1182        CoinWorkDouble t1=aa[n+1*BLOCK];
     1183        CoinWorkDouble t2=aa[n+2*BLOCK];
     1184        CoinWorkDouble t3=aa[n+3*BLOCK];
     1185        CoinWorkDouble a0;
    11671186        for (k=0;k<BLOCK;k++) {
    11681187          a0=aUnder[n+k*BLOCK]*work[k];
     
    11841203    aa+=BLOCK;
    11851204    for (i=0;i< nUnder;i++) {
    1186       longDouble t00=aa[i+0*BLOCK];
     1205      CoinWorkDouble t00=aa[i+0*BLOCK];
    11871206      for (k=0;k<BLOCK;k++) {
    1188         longDouble a00=aUnder[i+k*BLOCK]*work[k];
     1207        CoinWorkDouble a00=aUnder[i+k*BLOCK]*work[k];
    11891208        t00 -= a00 * above[j + k * BLOCK];
    11901209      }
     
    11961215/* Uses factorization to solve. */
    11971216void
    1198 ClpCholeskyDense::solve (double * region)
     1217ClpCholeskyDense::solve (CoinWorkDouble * region)
    11991218{
    12001219#ifdef CHOL_COMPARE
     
    12281247  }
    12291248#endif
    1230   //longDouble * xx = sparseFactor_;
    1231   //longDouble * yy = diagonal_;
    1232   //diagonal_ = sparseFactor_ + 40000;
    1233   //sparseFactor_=diagonal_ + numberRows_;
     1249  /*longDouble * xx = sparseFactor_;*/
     1250  /*longDouble * yy = diagonal_;*/
     1251  /*diagonal_ = sparseFactor_ + 40000;*/
     1252  /*sparseFactor_=diagonal_ + numberRows_;*/
    12341253  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
    1235   // later align on boundary
     1254  /* later align on boundary*/
    12361255  longDouble * a = sparseFactor_+BLOCKSQ*numberBlocks;
    12371256  int iBlock;
     
    12611280    aa+=BLOCKSQ;
    12621281  }
    1263   // do diagonal outside
     1282  /* do diagonal outside*/
    12641283  for (iColumn=0;iColumn<numberRows_;iColumn++)
    12651284    region[iColumn] *= diagonal_[iColumn];
     
    12931312  if (numberRows_<200) {
    12941313    for (int i=0;i<numberRows_;i++) {
    1295       assert(fabs(region[i]-region2[i])<1.0e-3);
     1314      assert(CoinAbs(region[i]-region2[i])<1.0e-3);
    12961315    }
    12971316    delete [] region2;
     
    12991318#endif
    13001319}
    1301 // Forward part of solve 1
     1320/* Forward part of solve 1*/
    13021321void
    1303 ClpCholeskyDense::solveF1(longDouble * a,int n,double * region)
     1322ClpCholeskyDense::solveF1(longDouble * a,int n,CoinWorkDouble * region)
    13041323{
    13051324  int j, k;
    1306   longDouble t00;
     1325  CoinWorkDouble t00;
    13071326  for (j = 0; j < n; j ++) {
    13081327    t00 = region[j];
     
    13101329      t00 -= region[k] * a[j + k * BLOCK];
    13111330    }
    1312     //t00*=a[j + j * BLOCK];
     1331    /*t00*=a[j + j * BLOCK];*/
    13131332    region[j] = t00;
    13141333  }
    13151334}
    1316 // Forward part of solve 2
     1335/* Forward part of solve 2*/
    13171336void
    1318 ClpCholeskyDense::solveF2(longDouble * a,int n,double * region, double * region2)
     1337ClpCholeskyDense::solveF2(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2)
    13191338{
    13201339  int j, k;
     
    13221341  if (n==BLOCK) {
    13231342    for (k = 0; k < BLOCK; k+=4) {
    1324       longDouble t0 = region2[0];
    1325       longDouble t1 = region2[1];
    1326       longDouble t2 = region2[2];
    1327       longDouble t3 = region2[3];
     1343      CoinWorkDouble t0 = region2[0];
     1344      CoinWorkDouble t1 = region2[1];
     1345      CoinWorkDouble t2 = region2[2];
     1346      CoinWorkDouble t3 = region2[3];
    13281347      t0 -= region[0] * a[0 + 0 * BLOCK];
    13291348      t1 -= region[0] * a[1 + 0 * BLOCK];
     
    14161435#endif
    14171436    for (k = 0; k < n; ++k) {
    1418       longDouble t00 = region2[k];
     1437      CoinWorkDouble t00 = region2[k];
    14191438      for (j = 0; j < BLOCK; j ++) {
    14201439        t00 -= region[j] * a[k + j * BLOCK];
     
    14261445#endif
    14271446}
    1428 // Backward part of solve 1
     1447/* Backward part of solve 1*/
    14291448void
    1430 ClpCholeskyDense::solveB1(longDouble * a,int n,double * region)
     1449ClpCholeskyDense::solveB1(longDouble * a,int n,CoinWorkDouble * region)
    14311450{
    14321451  int j, k;
    1433   longDouble t00;
     1452  CoinWorkDouble t00;
    14341453  for (j = n-1; j >=0; j --) {
    14351454    t00 = region[j];
     
    14371456      t00 -= region[k] * a[k + j * BLOCK];
    14381457    }
    1439     //t00*=a[j + j * BLOCK];
     1458    /*t00*=a[j + j * BLOCK];*/
    14401459    region[j] = t00;
    14411460  }
    14421461}
    1443 // Backward part of solve 2
     1462/* Backward part of solve 2*/
    14441463void
    1445 ClpCholeskyDense::solveB2(longDouble * a,int n,double * region, double * region2)
     1464ClpCholeskyDense::solveB2(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2)
    14461465{
    14471466  int j, k;
     
    14491468  if (n==BLOCK) {
    14501469    for (j = 0; j < BLOCK; j +=4) {
    1451       longDouble t0 = region[0];
    1452       longDouble t1 = region[1];
    1453       longDouble t2 = region[2];
    1454       longDouble t3 = region[3];
     1470      CoinWorkDouble t0 = region[0];
     1471      CoinWorkDouble t1 = region[1];
     1472      CoinWorkDouble t2 = region[2];
     1473      CoinWorkDouble t3 = region[3];
    14551474      t0 -= region2[0] * a[0 + 0*BLOCK];
    14561475      t1 -= region2[0] * a[0 + 1*BLOCK];
     
    15441563#endif
    15451564    for (j = 0; j < BLOCK; j ++) {
    1546       longDouble t00 = region[j];
     1565      CoinWorkDouble t00 = region[j];
    15471566      for (k = 0; k < n; ++k) {
    15481567        t00 -= region2[k] * a[k + j * BLOCK];
     
    15541573#endif
    15551574}
    1556 /* Uses factorization to solve. */
    1557 void
    1558 ClpCholeskyDense::solveLong (longDouble * region)
    1559 {
    1560   int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
    1561   // later align on boundary
    1562   longDouble * a = sparseFactor_+BLOCKSQ*numberBlocks;
    1563   int iBlock;
    1564   longDouble * aa = a;
    1565   int iColumn;
    1566   for (iBlock=0;iBlock<numberBlocks;iBlock++) {
    1567     int nChunk;
    1568     int jBlock;
    1569     int iDo = iBlock*BLOCK;
    1570     int base=iDo;
    1571     if (iDo+BLOCK>numberRows_) {
    1572       nChunk=numberRows_-iDo;
    1573     } else {
    1574       nChunk=BLOCK;
    1575     }
    1576     solveF1Long(aa,nChunk,region+iDo);
    1577     for (jBlock=iBlock+1;jBlock<numberBlocks;jBlock++) {
    1578       base+=BLOCK;
    1579       aa+=BLOCKSQ;
    1580       if (base+BLOCK>numberRows_) {
    1581         nChunk=numberRows_-base;
    1582       } else {
    1583         nChunk=BLOCK;
    1584       }
    1585       solveF2Long(aa,nChunk,region+iDo,region+base);
    1586     }
    1587     aa+=BLOCKSQ;
    1588   }
    1589   // do diagonal outside
    1590   for (iColumn=0;iColumn<numberRows_;iColumn++)
    1591     region[iColumn] *= diagonal_[iColumn];
    1592   int offset=((numberBlocks*(numberBlocks+1))>>1);
    1593   aa = a+number_entries(offset-1);
    1594   int lBase=(numberBlocks-1)*BLOCK;
    1595   for (iBlock=numberBlocks-1;iBlock>=0;iBlock--) {
    1596     int nChunk;
    1597     int jBlock;
    1598     int triBase=iBlock*BLOCK;
    1599     int iBase=lBase;
    1600     for (jBlock=iBlock+1;jBlock<numberBlocks;jBlock++) {
    1601       if (iBase+BLOCK>numberRows_) {
    1602         nChunk=numberRows_-iBase;
    1603       } else {
    1604         nChunk=BLOCK;
    1605       }
    1606       solveB2Long(aa,nChunk,region+triBase,region+iBase);
    1607       iBase-=BLOCK;
    1608       aa-=BLOCKSQ;
    1609     }
    1610     if (triBase+BLOCK>numberRows_) {
    1611       nChunk=numberRows_-triBase;
    1612     } else {
    1613       nChunk=BLOCK;
    1614     }
    1615     solveB1Long(aa,nChunk,region+triBase);
    1616     aa-=BLOCKSQ;
    1617   }
    1618 }
    1619 // Forward part of solve 1
    1620 void
    1621 ClpCholeskyDense::solveF1Long(longDouble * a,int n,longDouble * region)
    1622 {
    1623   int j, k;
    1624   longDouble t00;
    1625   for (j = 0; j < n; j ++) {
    1626     t00 = region[j];
    1627     for (k = 0; k < j; ++k) {
    1628       t00 -= region[k] * a[j + k * BLOCK];
    1629     }
    1630     //t00*=a[j + j * BLOCK];
    1631     region[j] = t00;
    1632   }
    1633 }
    1634 // Forward part of solve 2
    1635 void
    1636 ClpCholeskyDense::solveF2Long(longDouble * a,int n,longDouble * region, longDouble * region2)
    1637 {
    1638   int j, k;
    1639 #ifdef BLOCKUNROLL
    1640   if (n==BLOCK) {
    1641     for (k = 0; k < BLOCK; k+=4) {
    1642       longDouble t0 = region2[0];
    1643       longDouble t1 = region2[1];
    1644       longDouble t2 = region2[2];
    1645       longDouble t3 = region2[3];
    1646       t0 -= region[0] * a[0 + 0 * BLOCK];
    1647       t1 -= region[0] * a[1 + 0 * BLOCK];
    1648       t2 -= region[0] * a[2 + 0 * BLOCK];
    1649       t3 -= region[0] * a[3 + 0 * BLOCK];
    1650 
    1651       t0 -= region[1] * a[0 + 1 * BLOCK];
    1652       t1 -= region[1] * a[1 + 1 * BLOCK];
    1653       t2 -= region[1] * a[2 + 1 * BLOCK];
    1654       t3 -= region[1] * a[3 + 1 * BLOCK];
    1655 
    1656       t0 -= region[2] * a[0 + 2 * BLOCK];
    1657       t1 -= region[2] * a[1 + 2 * BLOCK];
    1658       t2 -= region[2] * a[2 + 2 * BLOCK];
    1659       t3 -= region[2] * a[3 + 2 * BLOCK];
    1660 
    1661       t0 -= region[3] * a[0 + 3 * BLOCK];
    1662       t1 -= region[3] * a[1 + 3 * BLOCK];
    1663       t2 -= region[3] * a[2 + 3 * BLOCK];
    1664       t3 -= region[3] * a[3 + 3 * BLOCK];
    1665 
    1666       t0 -= region[4] * a[0 + 4 * BLOCK];
    1667       t1 -= region[4] * a[1 + 4 * BLOCK];
    1668       t2 -= region[4] * a[2 + 4 * BLOCK];
    1669       t3 -= region[4] * a[3 + 4 * BLOCK];
    1670 
    1671       t0 -= region[5] * a[0 + 5 * BLOCK];
    1672       t1 -= region[5] * a[1 + 5 * BLOCK];
    1673       t2 -= region[5] * a[2 + 5 * BLOCK];
    1674       t3 -= region[5] * a[3 + 5 * BLOCK];
    1675 
    1676       t0 -= region[6] * a[0 + 6 * BLOCK];
    1677       t1 -= region[6] * a[1 + 6 * BLOCK];
    1678       t2 -= region[6] * a[2 + 6 * BLOCK];
    1679       t3 -= region[6] * a[3 + 6 * BLOCK];
    1680 
    1681       t0 -= region[7] * a[0 + 7 * BLOCK];
    1682       t1 -= region[7] * a[1 + 7 * BLOCK];
    1683       t2 -= region[7] * a[2 + 7 * BLOCK];
    1684       t3 -= region[7] * a[3 + 7 * BLOCK];
    1685 #if BLOCK>8
    1686       t0 -= region[8] * a[0 + 8 * BLOCK];
    1687       t1 -= region[8] * a[1 + 8 * BLOCK];
    1688       t2 -= region[8] * a[2 + 8 * BLOCK];
    1689       t3 -= region[8] * a[3 + 8 * BLOCK];
    1690 
    1691       t0 -= region[9] * a[0 + 9 * BLOCK];
    1692       t1 -= region[9] * a[1 + 9 * BLOCK];
    1693       t2 -= region[9] * a[2 + 9 * BLOCK];
    1694       t3 -= region[9] * a[3 + 9 * BLOCK];
    1695 
    1696       t0 -= region[10] * a[0 + 10 * BLOCK];
    1697       t1 -= region[10] * a[1 + 10 * BLOCK];
    1698       t2 -= region[10] * a[2 + 10 * BLOCK];
    1699       t3 -= region[10] * a[3 + 10 * BLOCK];
    1700 
    1701       t0 -= region[11] * a[0 + 11 * BLOCK];
    1702       t1 -= region[11] * a[1 + 11 * BLOCK];
    1703       t2 -= region[11] * a[2 + 11 * BLOCK];
    1704       t3 -= region[11] * a[3 + 11 * BLOCK];
    1705 
    1706       t0 -= region[12] * a[0 + 12 * BLOCK];
    1707       t1 -= region[12] * a[1 + 12 * BLOCK];
    1708       t2 -= region[12] * a[2 + 12 * BLOCK];
    1709       t3 -= region[12] * a[3 + 12 * BLOCK];
    1710 
    1711       t0 -= region[13] * a[0 + 13 * BLOCK];
    1712       t1 -= region[13] * a[1 + 13 * BLOCK];
    1713       t2 -= region[13] * a[2 + 13 * BLOCK];
    1714       t3 -= region[13] * a[3 + 13 * BLOCK];
    1715 
    1716       t0 -= region[14] * a[0 + 14 * BLOCK];
    1717       t1 -= region[14] * a[1 + 14 * BLOCK];
    1718       t2 -= region[14] * a[2 + 14 * BLOCK];
    1719       t3 -= region[14] * a[3 + 14 * BLOCK];
    1720 
    1721       t0 -= region[15] * a[0 + 15 * BLOCK];
    1722       t1 -= region[15] * a[1 + 15 * BLOCK];
    1723       t2 -= region[15] * a[2 + 15 * BLOCK];
    1724       t3 -= region[15] * a[3 + 15 * BLOCK];
    1725 #endif
    1726       region2[0] = t0;
    1727       region2[1] = t1;
    1728       region2[2] = t2;
    1729       region2[3] = t3;
    1730       region2+=4;
    1731       a+=4;
    1732     }
    1733   } else {
    1734 #endif
    1735     for (k = 0; k < n; ++k) {
    1736       longDouble t00 = region2[k];
    1737       for (j = 0; j < BLOCK; j ++) {
    1738         t00 -= region[j] * a[k + j * BLOCK];
    1739       }
    1740       region2[k] = t00;
    1741     }
    1742 #ifdef BLOCKUNROLL
    1743   }
    1744 #endif
    1745 }
    1746 // Backward part of solve 1
    1747 void
    1748 ClpCholeskyDense::solveB1Long(longDouble * a,int n,longDouble * region)
    1749 {
    1750   int j, k;
    1751   longDouble t00;
    1752   for (j = n-1; j >=0; j --) {
    1753     t00 = region[j];
    1754     for (k = j+1; k < n; ++k) {
    1755       t00 -= region[k] * a[k + j * BLOCK];