- Timestamp:
- Nov 5, 2009 7:34:07 AM (10 years ago)
- Location:
- stable/1.11
- Files:
-
- 112 edited
- 1 copied
Legend:
- Unmodified
- Added
- Removed
-
stable/1.11/Clp/src/CbcOrClpParam.cpp
r1404 r1458 1 /* $Id$ */ 1 2 // Copyright (C) 2002, International Business Machines 2 3 // Corporation and others. All Rights Reserved. … … 33 34 #else 34 35 static char coin_prompt[]="Clp:"; 36 #endif 37 #ifdef CLP_CILK 38 #ifndef CBC_THREAD 39 #define CBC_THREAD 40 #endif 35 41 #endif 36 42 #if WSSMP_BARRIER>=2 … … 67 73 action_(INVALID), 68 74 currentKeyWord_(-1), 69 display_( false),75 display_(0), 70 76 intValue_(-1), 71 77 doubleValue_(-1.0), … … 77 83 CbcOrClpParam::CbcOrClpParam (std::string name, std::string help, 78 84 double lower, double upper, CbcOrClpParameterType type, 79 booldisplay)85 int display) 80 86 : type_(type), 81 87 lowerIntValue_(0), … … 99 105 CbcOrClpParam::CbcOrClpParam (std::string name, std::string help, 100 106 int lower, int upper, CbcOrClpParameterType type, 101 booldisplay)107 int display) 102 108 : type_(type), 103 109 lowerDoubleValue_(0.0), … … 123 129 std::string firstValue, 124 130 CbcOrClpParameterType type,int whereUsed, 125 booldisplay)131 int display) 126 132 : type_(type), 127 133 lowerDoubleValue_(0.0), … … 147 153 CbcOrClpParam::CbcOrClpParam (std::string name, std::string help, 148 154 CbcOrClpParameterType type,int whereUsed, 149 booldisplay)155 int display) 150 156 : type_(type), 151 157 lowerDoubleValue_(0.0), … … 232 238 { 233 239 std::string::size_type shriekPos = name_.find('!'); 234 lengthName_ = (unsigned int)name_.length();240 lengthName_ = name_.length(); 235 241 if ( shriekPos==std::string::npos ) { 236 242 //does not contain '!' 237 243 lengthMatch_= lengthName_; 238 244 } else { 239 lengthMatch_= (unsigned int)shriekPos;245 lengthMatch_=shriekPos; 240 246 name_ = name_.substr(0,shriekPos)+name_.substr(shriekPos+1); 241 247 lengthName_--; 242 248 } 249 } 250 // Returns length of name for printing 251 int 252 CbcOrClpParam::lengthMatchName ( ) const 253 { 254 if (lengthName_==lengthMatch_) 255 return lengthName_; 256 else 257 return lengthName_+2; 243 258 } 244 259 // Insert string (only valid for keywords) … … 285 300 CbcOrClpParam::parameterOption ( std::string check ) const 286 301 { 287 size_t numberItems = definedKeyWords_.size();302 int numberItems = definedKeyWords_.size(); 288 303 if (!numberItems) { 289 304 return -1; 290 305 } else { 291 size_t whichItem=0;306 int whichItem=0; 292 307 unsigned int it; 293 308 for (it=0;it<definedKeyWords_.size();it++) { 294 309 std::string thisOne = definedKeyWords_[it]; 295 std::string::size_type shriekPos = thisOne.find('!');296 std::string::size_typelength1 = thisOne.length();297 std::string::size_typelength2 = length1;310 std::string::size_type shriekPos = thisOne.find('!'); 311 unsigned int length1 = thisOne.length(); 312 unsigned int length2 = length1; 298 313 if ( shriekPos!=std::string::npos ) { 299 314 //contains '!' … … 304 319 } 305 320 if (check.length()<=length1&&length2<=check.length()) { 306 std::string::size_typei;321 unsigned int i; 307 322 for (i=0;i<check.length();i++) { 308 323 if (tolower(thisOne[i])!=tolower(check[i])) … … 319 334 } 320 335 if (whichItem<numberItems) 321 return (int)whichItem;336 return whichItem; 322 337 else 323 338 return -1; … … 364 379 void CoinReadPrintit(const char * input) 365 380 { 366 size_t length =strlen(input);381 int length =strlen(input); 367 382 char temp[101]; 368 size_t i;383 int i; 369 384 int n=0; 370 385 for (i=0;i<length;i++) { … … 644 659 #ifdef COIN_HAS_CBC 645 660 double 646 CbcOrClpParam::doubleParameter (OsiSolverInterface * model) const 661 CbcOrClpParam::doubleParameter (OsiSolverInterface * 662 #ifndef NDEBUG 663 model 664 #endif 665 ) const 647 666 { 648 667 double value=0.0; … … 858 877 model.setNumberAnalyzeIterations(value); 859 878 break; 860 case NUMBERMINI:861 oldValue=model.sizeMiniTree();862 model.setSizeMiniTree(value);863 break;864 879 case CUTPASSINTREE: 865 880 oldValue=model.getMaximumCutPasses(); … … 912 927 case NUMBERANALYZE: 913 928 value=model.numberAnalyzeIterations(); 914 break;915 case NUMBERMINI:916 value=model.sizeMiniTree();917 929 break; 918 930 case CUTPASSINTREE: … … 999 1011 static char * where=NULL; 1000 1012 extern int CbcOrClpRead_mode; 1013 int CbcOrClpEnvironmentIndex=-1; 1014 static 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 } 1001 1048 extern FILE * CbcOrClpReadCommand; 1002 1049 // Simple read stuff … … 1072 1119 while (field=="EOL") { 1073 1120 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 } 1076 1133 if (field=="-") { 1077 1134 std::cout<<"Switching to line mode"<<std::endl; … … 1082 1139 // now allow std::cout<<"skipping non-command "<<field<<std::endl; 1083 1140 // field="EOL"; // skip 1084 } else {1141 } else if (CbcOrClpEnvironmentIndex<0) { 1085 1142 // special dispensation - taken as -import name 1086 1143 CbcOrClpRead_mode--; … … 1119 1176 if (afterEquals=="") { 1120 1177 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 } 1129 1191 } 1130 1192 } else { … … 1145 1207 if (afterEquals=="") { 1146 1208 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 } 1150 1217 } 1151 1218 } else { … … 1156 1223 afterEquals = ""; 1157 1224 } 1158 size_t value=0;1225 int value=0; 1159 1226 //std::cout<<field<<std::endl; 1160 1227 if (field!="EOL") { … … 1172 1239 *valid=2; 1173 1240 } 1174 return (int)value;1241 return value; 1175 1242 } 1176 1243 double … … 1180 1247 if (afterEquals=="") { 1181 1248 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 } 1185 1257 } 1186 1258 } else { … … 1218 1290 numberParameters=0; 1219 1291 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); 1223 1295 parameters[numberParameters++]= 1224 1296 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 \ 1306 are printed out on ?." 1307 ); 1226 1308 #ifdef COIN_HAS_CBC 1227 1309 parameters[numberParameters++]= … … 1245 1327 basis anyway." 1246 1328 ); 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 1247 1341 parameters[numberParameters++]= 1248 1342 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); 1250 1344 parameters[numberParameters-1].append("on"); 1251 1345 parameters[numberParameters-1].setLonghelp … … 1253 1347 "If you think you may get odd objective values or large equality rows etc then\ 1254 1348 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." 1256 1350 ); 1257 1351 parameters[numberParameters++]= … … 1261 1355 ( 1262 1356 "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 \1357 corrector algorithm. You may want to link in an alternative \ 1358 ordering and factorization. It will also solve models \ 1265 1359 with quadratic objectives." 1266 1360 … … 1274 1368 directory given by 'directory'. A name of '$' will use the previous value for the name. This\ 1275 1369 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." 1277 1371 ); 1278 1372 parameters[numberParameters++]= … … 1287 1381 parameters[numberParameters++]= 1288 1382 CbcOrClpParam("biasLU","Whether factorization biased towards U", 1289 "UU",BIASLU,2, false);1383 "UU",BIASLU,2,0); 1290 1384 parameters[numberParameters-1].append("UX"); 1291 1385 parameters[numberParameters-1].append("LX"); … … 1312 1406 were obtained fairly early in the search so the important point is to select the best variable to branch on. \ 1313 1407 See 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"); 1408 and trustPseudoCosts parameters. If cuts did a good job, then you may wish to \ 1409 have 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"); 1321 1420 parameters[numberParameters++]= 1322 1421 CbcOrClpParam("chol!esky","Which cholesky algorithm", … … 1327 1426 parameters[numberParameters-1].append("fudge!Long"); 1328 1427 parameters[numberParameters-1].append("wssmp"); 1329 #define REAL_BARRIER1330 1428 #else 1331 1429 parameters[numberParameters-1].append("fudge!Long_dummy"); … … 1334 1432 #ifdef UFL_BARRIER 1335 1433 parameters[numberParameters-1].append("Uni!versityOfFlorida"); 1336 #define REAL_BARRIER1337 1434 #else 1338 1435 parameters[numberParameters-1].append("Uni!versityOfFlorida_dummy"); … … 1340 1437 #ifdef TAUCS_BARRIER 1341 1438 parameters[numberParameters-1].append("Taucs"); 1342 #define REAL_BARRIER1343 1439 #else 1344 1440 parameters[numberParameters-1].append("Taucs_dummy"); … … 1346 1442 #ifdef MUMPS_BARRIER 1347 1443 parameters[numberParameters-1].append("Mumps"); 1348 #define REAL_BARRIER1349 1444 #else 1350 1445 parameters[numberParameters-1].append("Mumps_dummy"); … … 1353 1448 ( 1354 1449 "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 \ 1450 The native ordering and factorization is not state of the art, although acceptable. \ 1451 You may want to link in one from another source. See Makefile.locations for some \ 1356 1452 possibilities." 1357 1453 ); … … 1366 1462 parameters[numberParameters-1].append("forceOn"); 1367 1463 parameters[numberParameters-1].append("onglobal"); 1368 parameters[numberParameters-1].append("rootglobal");1369 1464 parameters[numberParameters-1].setLonghelp 1370 1465 ( … … 1382 1477 "This switches on a heuristic which does branch and cut on the problem given by just \ 1383 1478 using variables which have appeared in one or more solutions. \ 1479 It obviously only tries after two or more solutions. \ 1480 See 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 \ 1491 fixing variables which have same value in two or more solutions. \ 1384 1492 It obviously only tries after two or more solutions. \ 1385 1493 See Rounding for meaning of on,both,before" … … 1405 1513 parameters[numberParameters-1].setLonghelp 1406 1514 ( 1407 " If the user has Cplex, but wants to use some of Cbc's heuristics \1515 " If the user has Cplex, but wants to use some of Cbc's heuristics \ 1408 1516 then you can! If this is on, then Cbc will get to the root node and then \ 1409 1517 hand over to Cplex. If heuristics find a solution this can be significantly \ … … 1416 1524 parameters[numberParameters++]= 1417 1525 CbcOrClpParam("cpp!Generate","Generates C++ code", 1418 -1,50000,CPP );1526 -1,50000,CPP,1); 1419 1527 parameters[numberParameters-1].setLonghelp 1420 1528 ( … … 1423 1531 0 gives simplest driver, 1 generates saves and restores, 2 \ 1424 1532 generates saves and restores even for variables at default value. \ 1425 4 bit in cbc generates size dependent code rather than computed values." 1533 4 bit in cbc generates size dependent code rather than computed values. \ 1534 This is now deprecated as you can call stand-alone solver - see \ 1535 Cbc/examples/driver4.cpp." 1426 1536 ); 1427 1537 #ifdef COIN_HAS_CLP … … 1438 1548 "If crash is set on and there is an all slack basis then Clp will flip or put structural\ 1439 1549 variables into basis with the aim of getting dual feasible. On the whole dual seems to be\ 1440 better without it and there a lernative 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'. \ 1441 1551 I have also added a variant due to Solow and Halim which is as on but just flip."); 1442 1552 parameters[numberParameters++]= … … 1458 1568 parameters[numberParameters++]= 1459 1569 CbcOrClpParam("csv!Statistics","Create one line of statistics", 1460 CSVSTATISTICS,2, false);1570 CSVSTATISTICS,2,1); 1461 1571 parameters[numberParameters-1].setLonghelp 1462 1572 ( … … 1474 1584 every so often. The original method was every so many nodes but it is more logical \ 1475 1585 to do it whenever depth in tree is a multiple of K. This option does that and defaults \ 1476 to -1 (off)." 1586 to -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. \ 1595 Any value >0 says that all cuts <= this length can be generated both at \ 1596 root node and in tree. 0 says to use some dynamic lengths. If value >=10,000,000 \ 1597 then the length in tree is value%10000000 - so 10000100 means unlimited length \ 1598 at root and 100 in tree." 1477 1599 ); 1478 1600 parameters[numberParameters-1].setIntValue(-1); … … 1502 1624 parameters[numberParameters++]= 1503 1625 CbcOrClpParam("debug!In","read valid solution from file", 1504 DEBUG,7, false);1626 DEBUG,7,1); 1505 1627 parameters[numberParameters-1].setLonghelp 1506 1628 ( … … 1518 1640 parameters[numberParameters++]= 1519 1641 CbcOrClpParam("dense!Threshold","Whether to use dense factorization", 1520 -1,10000,DENSE, false);1642 -1,10000,DENSE,1); 1521 1643 parameters[numberParameters-1].setLonghelp 1522 1644 ( … … 1528 1650 #ifdef COIN_HAS_CBC 1529 1651 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 \ 1658 means 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); 1532 1663 parameters[numberParameters-1].setDoubleValue(0.0); 1533 1664 parameters[numberParameters++]= 1534 CbcOrClpParam("dextra 2","Extra double parameter 2",1535 -COIN_DBL_MAX,COIN_DBL_MAX,DEXTRA 2,false);1665 CbcOrClpParam("dextra4","Extra double parameter 4", 1666 -COIN_DBL_MAX,COIN_DBL_MAX,DEXTRA4,0); 1536 1667 parameters[numberParameters-1].setDoubleValue(0.0); 1537 1668 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++]=1546 1669 CbcOrClpParam("dextra5","Extra double parameter 5", 1547 -COIN_DBL_MAX,COIN_DBL_MAX,DEXTRA5, false);1670 -COIN_DBL_MAX,COIN_DBL_MAX,DEXTRA5,0); 1548 1671 parameters[numberParameters-1].setDoubleValue(0.0); 1549 1672 parameters[numberParameters++]= … … 1580 1703 parameters[numberParameters++]= 1581 1704 CbcOrClpParam("dirSample","Set directory where the COIN-OR sample problems are.", 1582 DIRSAMPLE );1705 DIRSAMPLE,7,1); 1583 1706 parameters[numberParameters-1].setLonghelp 1584 1707 ( … … 1590 1713 parameters[numberParameters++]= 1591 1714 CbcOrClpParam("dirNetlib","Set directory where the netlib problems are.", 1592 DIRNETLIB );1715 DIRNETLIB,7,1); 1593 1716 parameters[numberParameters-1].setLonghelp 1594 1717 ( … … 1602 1725 parameters[numberParameters++]= 1603 1726 CbcOrClpParam("dirMiplib","Set directory where the miplib 2003 problems are.", 1604 DIRMIPLIB );1727 DIRMIPLIB,7,1); 1605 1728 parameters[numberParameters-1].setLonghelp 1606 1729 ( … … 1615 1738 parameters[numberParameters++]= 1616 1739 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." 1621 1749 ); 1622 1750 parameters[numberParameters-1].setIntValue(-1); … … 1676 1804 ( 1677 1805 "Normally heuristics are done in branch and bound. It may be useful to do them outside. \ 1678 Doing this may also set cutoff." 1806 Only those heuristics with 'both' or 'before' set will run. \ 1807 Doing this may also set cutoff, which can help with preprocessing." 1679 1808 ); 1680 1809 #endif … … 1698 1827 parameters[numberParameters++]= 1699 1828 CbcOrClpParam("dualize","Solves dual reformulation", 1700 0,3,DUALIZE, false);1829 0,3,DUALIZE,1); 1701 1830 parameters[numberParameters-1].setLonghelp 1702 1831 ( … … 1705 1834 parameters[numberParameters++]= 1706 1835 CbcOrClpParam("dualP!ivot","Dual pivot choice algorithm", 1707 "auto!matic",DUALPIVOT );1836 "auto!matic",DUALPIVOT,7,1); 1708 1837 parameters[numberParameters-1].append("dant!zig"); 1709 1838 parameters[numberParameters-1].append("partial"); … … 1757 1886 ); 1758 1887 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++]= 1759 1895 CbcOrClpParam("error!sAllowed","Whether to allow import errors", 1760 1896 "off",ERRORSALLOWED,3); … … 1777 1913 parameters[numberParameters++]= 1778 1914 CbcOrClpParam("exp!eriment","Whether to use testing features", 1779 -1,200,EXPERIMENT, false);1915 -1,200,EXPERIMENT,0); 1780 1916 parameters[numberParameters-1].setLonghelp 1781 1917 ( … … 1792 1928 "This will write an MPS format file to the given file name. It will use the default\ 1793 1929 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'. \ 1931 It 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." 1795 1932 ); 1796 1933 #ifdef COIN_HAS_CBC 1797 1934 parameters[numberParameters++]= 1798 1935 CbcOrClpParam("extra1","Extra integer parameter 1", 1799 -COIN_INT_MAX,COIN_INT_MAX,EXTRA1, false);1936 -COIN_INT_MAX,COIN_INT_MAX,EXTRA1,0); 1800 1937 parameters[numberParameters-1].setIntValue(-1); 1801 1938 parameters[numberParameters++]= 1802 1939 CbcOrClpParam("extra2","Extra integer parameter 2", 1803 -100,COIN_INT_MAX,EXTRA2, false);1940 -100,COIN_INT_MAX,EXTRA2,0); 1804 1941 parameters[numberParameters-1].setIntValue(-1); 1805 1942 parameters[numberParameters++]= 1806 1943 CbcOrClpParam("extra3","Extra integer parameter 3", 1807 -1,COIN_INT_MAX,EXTRA3, false);1944 -1,COIN_INT_MAX,EXTRA3,0); 1808 1945 parameters[numberParameters-1].setIntValue(-1); 1809 1946 parameters[numberParameters++]= 1810 1947 CbcOrClpParam("extra4","Extra integer parameter 4", 1811 - COIN_INT_MAX,COIN_INT_MAX,EXTRA4,false);1948 -1,COIN_INT_MAX,EXTRA4,0); 1812 1949 parameters[numberParameters-1].setIntValue(-1); 1950 parameters[numberParameters-1].setLonghelp 1951 ( 1952 "This switches on yet more special options!! \ 1953 The bottom digit is a strategy when to used shadow price stuff e.g. 3 \ 1954 means use until a solution is found. The next two digits say what sort \ 1955 of 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 ); 1813 1961 #endif 1814 1962 #ifdef COIN_HAS_CLP 1815 1963 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 \ 1972 other choices are a dense one, osl's or one designed for small problems." 1973 ); 1974 parameters[numberParameters++]= 1816 1975 CbcOrClpParam("fakeB!ound","All bounds <= this value - DEBUG", 1817 1.0,1.0e15,FAKEBOUND, false);1976 1.0,1.0e15,FAKEBOUND,0); 1818 1977 #ifdef COIN_HAS_CBC 1819 1978 parameters[numberParameters++]= … … 1833 1992 CbcOrClpParam("fix!OnDj","Try heuristic based on fixing variables with \ 1834 1993 reduced costs greater than this", 1835 -1.0e20,1.0e20,DJFIX, false);1994 -1.0e20,1.0e20,DJFIX,1); 1836 1995 parameters[numberParameters-1].setLonghelp 1837 1996 ( … … 1846 2005 parameters[numberParameters-1].append("ifmove"); 1847 2006 parameters[numberParameters-1].append("forceOn"); 1848 parameters[numberParameters-1].append("onglobal"); 1849 parameters[numberParameters-1].append("rootglobal"); 2007 parameters[numberParameters-1].append("onglobal"); 1850 2008 parameters[numberParameters-1].setLonghelp 1851 2009 ( … … 1861 2019 "-1 off. If 1 then tries to branch to solution given by AMPL or priorities file. \ 1862 2020 If 0 then just tries to set as best solution \ 1863 If 1 then also does that many nodes on fixed problem." 2021 If >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 \ 2030 about are fixed and if the preprocessed model is small enough a few nodes \ 2031 of branch and bound are done on reduced problem. Small problem has to be less than this fraction of original." 1864 2032 ); 1865 2033 #endif 1866 2034 parameters[numberParameters++]= 1867 2035 CbcOrClpParam("gamma!(Delta)","Whether to regularize barrier", 1868 "off",GAMMA,7, false);2036 "off",GAMMA,7,1); 1869 2037 parameters[numberParameters-1].append("on"); 1870 2038 parameters[numberParameters-1].append("gamma"); … … 1883 2051 parameters[numberParameters-1].append("forceOn"); 1884 2052 parameters[numberParameters-1].append("onglobal"); 1885 parameters[numberParameters-1].append(" rootglobal");2053 parameters[numberParameters-1].append("forceandglobal"); 1886 2054 parameters[numberParameters-1].append("forceLongOn"); 1887 2055 parameters[numberParameters-1].append("long"); … … 1929 2097 parameters[numberParameters++]= 1930 2098 CbcOrClpParam("hOp!tions","Heuristic options", 1931 -9999999,9999999,HOPTIONS );2099 -9999999,9999999,HOPTIONS,1); 1932 2100 parameters[numberParameters-1].setLonghelp 1933 2101 ( … … 1941 2109 parameters[numberParameters++]= 1942 2110 CbcOrClpParam("hot!StartMaxIts","Maximum iterations on hot start", 1943 0,COIN_INT_MAX,MAXHOTITS ,false);2111 0,COIN_INT_MAX,MAXHOTITS); 1944 2112 #endif 1945 2113 #ifdef COIN_HAS_CLP 1946 2114 parameters[numberParameters++]= 1947 2115 CbcOrClpParam("idiot!Crash","Whether to try idiot crash", 1948 -1,999999 ,IDIOT);2116 -1,99999999,IDIOT); 1949 2117 parameters[numberParameters-1].setLonghelp 1950 2118 ( … … 1963 2131 directory given by 'directory'. A name of '$' will use the previous value for the name. This\ 1964 2132 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'. \ 2134 If 'keepnames' is off, then names are dropped -> Rnnnnnnn and Cnnnnnnn." 1966 2135 ); 1967 2136 #ifdef COIN_HAS_CBC … … 1980 2149 CbcOrClpParam("inf!easibilityWeight","Each integer infeasibility is expected \ 1981 2150 to cost this much", 1982 0.0,1.0e20,INFEASIBILITYWEIGHT );2151 0.0,1.0e20,INFEASIBILITYWEIGHT,1); 1983 2152 parameters[numberParameters-1].setLonghelp 1984 2153 ( … … 2014 2183 parameters[numberParameters++]= 2015 2184 CbcOrClpParam("KKT","Whether to use KKT factorization", 2016 "off",KKT,7, false);2185 "off",KKT,7,1); 2017 2186 parameters[numberParameters-1].append("on"); 2018 2187 #endif … … 2026 2195 parameters[numberParameters-1].append("forceOn"); 2027 2196 parameters[numberParameters-1].append("onglobal"); 2028 parameters[numberParameters-1].append(" rootglobal");2197 parameters[numberParameters-1].append("forceandglobal"); 2029 2198 parameters[numberParameters-1].setLonghelp 2030 2199 ( … … 2098 2267 "This can be used for testing purposes. The corresponding library call\n\ 2099 2268 \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." 2101 2270 ); 2102 2271 #endif … … 2115 2284 parameters[numberParameters-1].setLonghelp 2116 2285 ( 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. \ 2287 This is checked every node in tree, so it is possible to get more solutions from heuristics." 2118 2288 ); 2119 2289 #endif … … 2130 2300 parameters[numberParameters++]= 2131 2301 CbcOrClpParam("mipO!ptions","Dubious options for mip", 2132 0,COIN_INT_MAX,MIPOPTIONS, false);2302 0,COIN_INT_MAX,MIPOPTIONS,0); 2133 2303 parameters[numberParameters++]= 2134 2304 CbcOrClpParam("more!MipOptions","More dubious options for mip", 2135 -1,COIN_INT_MAX,MOREMIPOPTIONS, false);2305 -1,COIN_INT_MAX,MOREMIPOPTIONS,0); 2136 2306 parameters[numberParameters++]= 2137 2307 CbcOrClpParam("mixed!IntegerRoundingCuts","Whether to use Mixed Integer Rounding cuts", … … 2142 2312 parameters[numberParameters-1].append("forceOn"); 2143 2313 parameters[numberParameters-1].append("onglobal"); 2144 parameters[numberParameters-1].append("rootglobal");2145 2314 parameters[numberParameters-1].setLonghelp 2146 2315 ( … … 2158 2327 but this program turns this off to make it look more friendly. It can be useful\ 2159 2328 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." 2161 2330 ); 2162 2331 #ifdef COIN_HAS_CBC 2163 2332 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); 2171 2347 parameters[numberParameters++]= 2172 2348 CbcOrClpParam("miplib","Do some of miplib test set", 2173 MIPLIB,3 );2349 MIPLIB,3,1); 2174 2350 parameters[numberParameters++]= 2175 2351 CbcOrClpParam("naive!Heuristics","Whether to try some stupid heuristic", 2176 "off",NAIVE );2352 "off",NAIVE,7,1); 2177 2353 parameters[numberParameters-1].append("on"); 2178 2354 parameters[numberParameters-1].append("both"); … … 2181 2357 ( 2182 2358 "Really silly stuff e.g. fix all integers with costs to zero!. \ 2183 Do optionsdoes heuristic before preprocessing" );2359 Doh option does heuristic before preprocessing" ); 2184 2360 #endif 2185 2361 #ifdef COIN_HAS_CLP 2186 2362 parameters[numberParameters++]= 2187 2363 CbcOrClpParam("netlib","Solve entire netlib test set", 2188 NETLIB_EITHER,3 );2364 NETLIB_EITHER,3,1); 2189 2365 parameters[numberParameters-1].setLonghelp 2190 2366 ( … … 2192 2368 The user can set options before e.g. clp -presolve off -netlib" 2193 2369 ); 2194 #ifdef REAL_BARRIER2195 2370 parameters[numberParameters++]= 2196 2371 CbcOrClpParam("netlibB!arrier","Solve entire netlib test set with barrier", 2197 NETLIB_BARRIER,3 );2372 NETLIB_BARRIER,3,1); 2198 2373 parameters[numberParameters-1].setLonghelp 2199 2374 ( … … 2201 2376 The user can set options before e.g. clp -kkt on -netlib" 2202 2377 ); 2203 #endif2204 2378 parameters[numberParameters++]= 2205 2379 CbcOrClpParam("netlibD!ual","Solve entire netlib test set (dual)", 2206 NETLIB_DUAL,3 );2380 NETLIB_DUAL,3,1); 2207 2381 parameters[numberParameters-1].setLonghelp 2208 2382 ( … … 2212 2386 parameters[numberParameters++]= 2213 2387 CbcOrClpParam("netlibP!rimal","Solve entire netlib test set (primal)", 2214 NETLIB_PRIMAL,3 );2388 NETLIB_PRIMAL,3,1); 2215 2389 parameters[numberParameters-1].setLonghelp 2216 2390 ( … … 2220 2394 parameters[numberParameters++]= 2221 2395 CbcOrClpParam("netlibT!une","Solve entire netlib test set with 'best' algorithm", 2222 NETLIB_TUNE,3 );2396 NETLIB_TUNE,3,1); 2223 2397 parameters[numberParameters-1].setLonghelp 2224 2398 ( … … 2230 2404 parameters[numberParameters++]= 2231 2405 CbcOrClpParam("network","Tries to make network matrix", 2232 NETWORK,7, false);2406 NETWORK,7,0); 2233 2407 parameters[numberParameters-1].setLonghelp 2234 2408 ( … … 2257 2431 parameters[numberParameters++]= 2258 2432 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); 2260 2434 parameters[numberParameters-1].setLonghelp 2261 2435 ( … … 2266 2440 parameters[numberParameters++]= 2267 2441 CbcOrClpParam("objective!Scale","Scale factor to apply to objective", 2268 -1.0e20,1.0e20,OBJSCALE, false);2442 -1.0e20,1.0e20,OBJSCALE,1); 2269 2443 parameters[numberParameters-1].setLonghelp 2270 2444 ( 2271 2445 "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." 2273 2447 ); 2274 2448 parameters[numberParameters-1].setDoubleValue(1.0); … … 2277 2451 parameters[numberParameters++]= 2278 2452 CbcOrClpParam("outDup!licates","takes duplicate rows etc out of integer model", 2279 OUTDUPROWS,7, false);2453 OUTDUPROWS,7,0); 2280 2454 #endif 2281 2455 parameters[numberParameters++]= … … 2312 2486 parameters[numberParameters++]= 2313 2487 CbcOrClpParam("passP!resolve","How many passes in presolve", 2314 -200,100,PRESOLVEPASS, false);2488 -200,100,PRESOLVEPASS,1); 2315 2489 parameters[numberParameters-1].setLonghelp 2316 2490 ( … … 2332 2506 parameters[numberParameters++]= 2333 2507 CbcOrClpParam("pertV!alue","Method of perturbation", 2334 -5000,102,PERTVALUE, false);2508 -5000,102,PERTVALUE,1); 2335 2509 parameters[numberParameters++]= 2336 2510 CbcOrClpParam("perturb!ation","Whether to perturb problem", … … 2346 2520 parameters[numberParameters++]= 2347 2521 CbcOrClpParam("PFI","Whether to use Product Form of Inverse in simplex", 2348 "off",PFI,7, false);2522 "off",PFI,7,0); 2349 2523 parameters[numberParameters-1].append("on"); 2350 2524 parameters[numberParameters-1].setLonghelp … … 2355 2529 #ifdef COIN_HAS_CBC 2356 2530 parameters[numberParameters++]= 2357 CbcOrClpParam("pivot !AndFix","Whether to try Pivot and Fixheuristic",2358 "off",PIVOTAND FIX);2531 CbcOrClpParam("pivotAndC!omplement","Whether to try Pivot and Complement heuristic", 2532 "off",PIVOTANDCOMPLEMENT); 2359 2533 parameters[numberParameters-1].append("on"); 2360 2534 parameters[numberParameters-1].append("both"); … … 2363 2537 ( 2364 2538 "stuff needed. \ 2365 Do options does heuristic before preprocessing" ); 2539 Doh 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. \ 2549 Doh option does heuristic before preprocessing" ); 2366 2550 #endif 2367 2551 #ifdef COIN_HAS_CLP 2368 2552 parameters[numberParameters++]= 2369 2553 CbcOrClpParam("plus!Minus","Tries to make +- 1 matrix", 2370 PLUSMINUS,7, false);2554 PLUSMINUS,7,0); 2371 2555 parameters[numberParameters-1].setLonghelp 2372 2556 ( … … 2377 2561 parameters[numberParameters++]= 2378 2562 CbcOrClpParam("pO!ptions","Dubious print options", 2379 0,COIN_INT_MAX,PRINTOPTIONS, false);2563 0,COIN_INT_MAX,PRINTOPTIONS,1); 2380 2564 parameters[numberParameters-1].setIntValue(0); 2381 2565 parameters[numberParameters-1].setLonghelp … … 2385 2569 parameters[numberParameters++]= 2386 2570 CbcOrClpParam("preO!pt","Presolve options", 2387 0,COIN_INT_MAX,PRESOLVEOPTIONS, false);2571 0,COIN_INT_MAX,PRESOLVEOPTIONS,0); 2388 2572 #endif 2389 2573 parameters[numberParameters++]= … … 2436 2620 parameters[numberParameters++]= 2437 2621 CbcOrClpParam("primalP!ivot","Primal pivot choice algorithm", 2438 "auto!matic",PRIMALPIVOT );2622 "auto!matic",PRIMALPIVOT,7,1); 2439 2623 parameters[numberParameters-1].append("exa!ct"); 2440 2624 parameters[numberParameters-1].append("dant!zig"); … … 2536 2720 parameters[numberParameters-1].append("forceOn"); 2537 2721 parameters[numberParameters-1].append("onglobal"); 2538 parameters[numberParameters-1].append(" rootglobal");2722 parameters[numberParameters-1].append("forceonglobal"); 2539 2723 parameters[numberParameters-1].append("forceOnBut"); 2540 2724 parameters[numberParameters-1].append("forceOnStrong"); … … 2548 2732 ); 2549 2733 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 \ 2749 when solution found in feasibility pump" 2750 ); 2751 parameters[numberParameters++]= 2550 2752 CbcOrClpParam("pumpT!une","Dubious ideas for feasibility pump", 2551 2753 0,100000000,FPUMPTUNE); … … 2553 2755 ( 2554 2756 "This fine tunes Feasibility Pump \n\ 2757 \t>=10000000 use as objective weight switch\n\ 2555 2758 \t>=1000000 use as accumulate switch\n\ 2556 2759 \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. \ 2763 If accumulate is on then after a major pass, variables which have not moved \ 2764 are fixed and a small branch and bound is tried." 2560 2765 ); 2561 2766 parameters[numberParameters-1].setIntValue(0); … … 2578 2783 ( 2579 2784 "stuff needed. \ 2580 Do optionsdoes heuristic before preprocessing" );2785 Doh option does heuristic before preprocessing" ); 2581 2786 parameters[numberParameters++]= 2582 2787 CbcOrClpParam("ratio!Gap","Stop when gap between best possible and \ … … 2592 2797 parameters[numberParameters++]= 2593 2798 CbcOrClpParam("readS!tored","Import stored cuts from file", 2594 STOREDFILE,3, false);2799 STOREDFILE,3,0); 2595 2800 #endif 2596 2801 #ifdef COIN_HAS_CLP 2597 2802 parameters[numberParameters++]= 2598 2803 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); 2600 2805 parameters[numberParameters-1].setLonghelp 2601 2806 ( … … 2605 2810 parameters[numberParameters++]= 2606 2811 CbcOrClpParam("reallyS!cale","Scales model in place", 2607 REALLY_SCALE,7, false);2812 REALLY_SCALE,7,0); 2608 2813 #endif 2609 2814 #ifdef COIN_HAS_CBC … … 2636 2841 parameters[numberParameters++]= 2637 2842 CbcOrClpParam("restore!Model","Restore model from binary file", 2638 RESTORE );2843 RESTORE,7,1); 2639 2844 parameters[numberParameters-1].setLonghelp 2640 2845 ( … … 2645 2850 parameters[numberParameters++]= 2646 2851 CbcOrClpParam("reverse","Reverses sign of objective", 2647 REVERSE,7, false);2852 REVERSE,7,0); 2648 2853 parameters[numberParameters-1].setLonghelp 2649 2854 ( … … 2652 2857 parameters[numberParameters++]= 2653 2858 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); 2655 2860 parameters[numberParameters-1].setLonghelp 2656 2861 ( 2657 2862 "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." 2659 2864 ); 2660 2865 parameters[numberParameters-1].setDoubleValue(1.0); … … 2673 2878 ( 2674 2879 "This switches on Relaxation enforced neighborhood Search. \ 2675 on just does feasibility pump\2880 on just does 50 nodes \ 2676 2881 200 or 1000 does that many nodes. \ 2677 Do optionsdoes heuristic before preprocessing" );2882 Doh option does heuristic before preprocessing" ); 2678 2883 parameters[numberParameters++]= 2679 2884 CbcOrClpParam("Rins","Whether to try Relaxed Induced Neighborhood Search", … … 2686 2891 ( 2687 2892 "This switches on Relaxed induced neighborhood Search. \ 2688 Do optionsdoes heuristic before preprocessing" );2893 Doh option does heuristic before preprocessing" ); 2689 2894 parameters[numberParameters++]= 2690 2895 CbcOrClpParam("round!ingHeuristic","Whether to use Rounding heuristic", … … 2704 2909 parameters[numberParameters++]= 2705 2910 CbcOrClpParam("saveM!odel","Save model to binary file", 2706 SAVE );2911 SAVE,7,1); 2707 2912 parameters[numberParameters-1].setLonghelp 2708 2913 ( … … 2730 2935 parameters[numberParameters-1].append("geo!metric"); 2731 2936 parameters[numberParameters-1].append("auto!matic"); 2937 parameters[numberParameters-1].append("dynamic"); 2938 parameters[numberParameters-1].append("rows!only"); 2732 2939 parameters[numberParameters-1].setLonghelp 2733 2940 ( … … 2745 2952 ( 2746 2953 "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)." 2750 2955 ); 2751 2956 #else … … 2760 2965 parameters[numberParameters++]= 2761 2966 CbcOrClpParam("sleep","for debug", 2762 DUMMY,7, false);2967 DUMMY,7,0); 2763 2968 parameters[numberParameters-1].setLonghelp 2764 2969 ( … … 2768 2973 parameters[numberParameters++]= 2769 2974 CbcOrClpParam("slp!Value","Number of slp passes before primal", 2770 -1,50000,SLPVALUE, false);2975 -1,50000,SLPVALUE,1); 2771 2976 parameters[numberParameters-1].setLonghelp 2772 2977 ( … … 2777 2982 parameters[numberParameters++]= 2778 2983 CbcOrClpParam("small!Factorization","Whether to use small factorization", 2779 -1,10000,SMALLFACT, false);2984 -1,10000,SMALLFACT,1); 2780 2985 parameters[numberParameters-1].setLonghelp 2781 2986 ( … … 2815 3020 ); 2816 3021 parameters[numberParameters++]= 2817 CbcOrClpParam("slog!Level","Level of detail in Solver output",3022 CbcOrClpParam("slog!Level","Level of detail in (LP) Solver output", 2818 3023 -1,63,SOLVERLOGLEVEL); 2819 3024 parameters[numberParameters-1].setLonghelp 2820 3025 ( 2821 3026 "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'" 2823 3028 ); 2824 3029 #else … … 2836 3041 parameters[numberParameters++]= 2837 3042 CbcOrClpParam("spars!eFactor","Whether factorization treated as sparse", 2838 "on",SPARSEFACTOR,7, false);3043 "on",SPARSEFACTOR,7,0); 2839 3044 parameters[numberParameters-1].append("off"); 2840 3045 parameters[numberParameters++]= 2841 3046 CbcOrClpParam("special!Options","Dubious options for Simplex - see ClpSimplex.hpp", 2842 0,COIN_INT_MAX,SPECIALOPTIONS, false);3047 0,COIN_INT_MAX,SPECIALOPTIONS,0); 2843 3048 parameters[numberParameters++]= 2844 3049 CbcOrClpParam("sprint!Crash","Whether to try sprint crash", … … 2883 3088 ); 2884 3089 parameters[numberParameters-1].setIntValue(1); 3090 #ifdef CBC_KEEP_DEPRECATED 2885 3091 parameters[numberParameters++]= 2886 3092 CbcOrClpParam("strengthen","Create strengthened problem", … … 2890 3096 "This creates a new problem by applying the root node cuts. All tight constraints \ 2891 3097 will be in resulting problem" 2892 ); 3098 ); 3099 #endif 2893 3100 parameters[numberParameters++]= 2894 3101 CbcOrClpParam("strong!Branching","Number of variables to look at in strong branching", … … 2905 3112 parameters[numberParameters++]= 2906 3113 CbcOrClpParam("subs!titution","How long a column to substitute for in presolve", 2907 0,10000,SUBSTITUTION, false);3114 0,10000,SUBSTITUTION,0); 2908 3115 parameters[numberParameters-1].setLonghelp 2909 3116 ( … … 2916 3123 parameters[numberParameters++]= 2917 3124 CbcOrClpParam("testO!si","Test OsiObject stuff", 2918 -1,COIN_INT_MAX,TESTOSI, false);3125 -1,COIN_INT_MAX,TESTOSI,0); 2919 3126 #endif 2920 3127 #ifdef CBC_THREAD 2921 3128 parameters[numberParameters++]= 2922 3129 CbcOrClpParam("thread!s","Number of threads to try and use", 2923 -100,100000,THREADS, false);3130 -100,100000,THREADS,1); 2924 3131 parameters[numberParameters-1].setLonghelp 2925 3132 ( … … 2934 3141 CbcOrClpParam("tighten!Factor","Tighten bounds using this times largest \ 2935 3142 activity at continuous solution", 2936 1.0e-3,1.0e20,TIGHTENFACTOR, false);3143 1.0e-3,1.0e20,TIGHTENFACTOR,0); 2937 3144 parameters[numberParameters-1].setLonghelp 2938 3145 ( … … 2943 3150 parameters[numberParameters++]= 2944 3151 CbcOrClpParam("tightLP","Poor person's preSolve for now", 2945 TIGHTEN,7, false);3152 TIGHTEN,7,0); 2946 3153 #endif 2947 3154 #ifdef COIN_HAS_CBC … … 2958 3165 parameters[numberParameters++]= 2959 3166 CbcOrClpParam("tune!PreProcess","Dubious tuning parameters", 2960 0,20000000,PROCESSTUNE, false);3167 0,20000000,PROCESSTUNE,1); 2961 3168 parameters[numberParameters-1].setLonghelp 2962 3169 ( … … 2977 3184 parameters[numberParameters-1].append("forceOn"); 2978 3185 parameters[numberParameters-1].append("onglobal"); 2979 parameters[numberParameters-1].append(" rootglobal");3186 parameters[numberParameters-1].append("forceandglobal"); 2980 3187 parameters[numberParameters-1].append("forceLongOn"); 2981 3188 parameters[numberParameters-1].setLonghelp … … 2987 3194 parameters[numberParameters++]= 2988 3195 CbcOrClpParam("unitTest","Do unit test", 2989 UNITTEST,3 );3196 UNITTEST,3,1); 2990 3197 parameters[numberParameters-1].setLonghelp 2991 3198 ( … … 2994 3201 parameters[numberParameters++]= 2995 3202 CbcOrClpParam("userClp","Hand coded Clp stuff", 2996 USERCLP );3203 USERCLP,0,0); 2997 3204 parameters[numberParameters-1].setLonghelp 2998 3205 ( … … 3003 3210 parameters[numberParameters++]= 3004 3211 CbcOrClpParam("userCbc","Hand coded Cbc stuff", 3005 USERCBC );3212 USERCBC,0,0); 3006 3213 parameters[numberParameters-1].setLonghelp 3007 3214 ( 3008 3215 "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 ); 3216 Look for USERCBC in main driver and modify sample code. \ 3217 It 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. \ 3229 Doh option does heuristic before preprocessing" ); 3011 3230 #endif 3012 3231 parameters[numberParameters++]= 3013 3232 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 andClpPackedMatrix 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." 3019 3238 ); 3020 3239 parameters[numberParameters++]= 3021 3240 CbcOrClpParam("verbose","Switches on longer help on single ?", 3022 0, 15,VERBOSE,false);3241 0,31,VERBOSE,0); 3023 3242 parameters[numberParameters-1].setLonghelp 3024 3243 ( … … 3029 3248 parameters[numberParameters++]= 3030 3249 CbcOrClpParam("vub!heuristic","Type of vub heuristic", 3031 -2,20,VUBTRY, false);3250 -2,20,VUBTRY,0); 3032 3251 parameters[numberParameters-1].setLonghelp 3033 3252 ( … … 3035 3254 ); 3036 3255 parameters[numberParameters-1].setIntValue(-1); 3256 #ifdef ZERO_HALF_CUTS 3037 3257 parameters[numberParameters++]= 3038 3258 CbcOrClpParam("zero!HalfCuts","Whether to use zero half cuts", … … 3043 3263 parameters[numberParameters-1].append("forceOn"); 3044 3264 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) \ 3049 3268 See branchAndCut for information on options." 3050 ); 3269 ); 3270 #endif 3051 3271 #endif 3052 3272 assert(numberParameters<CBCMAXPARAMETERS); … … 3078 3298 int numberColumnsFile; 3079 3299 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"); 3083 3310 double * dualRowSolution = lpSolver->dualRowSolution(); 3084 3311 double * primalRowSolution = lpSolver->primalRowSolution(); … … 3103 3330 lpSolver->setObjectiveValue(objectiveValue); 3104 3331 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"); 3109 3344 } else { 3110 3345 std::cout<<"Mismatch on rows and/or columns - truncating"<<std::endl; 3111 3346 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); 3120 3363 delete [] temp; 3121 3364 } … … 3177 3420 int numberColumns=lpSolver->numberColumns(); 3178 3421 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"); 3182 3432 double * dualRowSolution = lpSolver->dualRowSolution(); 3183 3433 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"); 3186 3440 double * dualColumnSolution = lpSolver->dualColumnSolution(); 3187 3441 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"); 3190 3448 fclose(fp); 3191 3449 } else { -
stable/1.11/Clp/src/CbcOrClpParam.hpp
r1344 r1458 1 2 /* $Id$ */ 1 3 // Copyright (C) 2002, International Business Machines 2 4 // Corporation and others. All Rights Reserved. … … 53 55 54 56 DJFIX = 81, TIGHTENFACTOR,PRESOLVETOLERANCE,OBJSCALE2, 55 DEXTRA1, DEXTRA2, DEXTRA3, DEXTRA4, DEXTRA5,57 FAKEINCREMENT, FAKECUTOFF, ARTIFICIALCOST,DEXTRA3, SMALLBAB, DEXTRA4, DEXTRA5, 56 58 57 59 SOLVERLOGLEVEL=101, … … 64 66 65 67 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, 68 70 THREADS,CUTPASS,VUBTRY,DENSE,EXPERIMENT,DIVEOPT,STRATEGY,SMALLFACT, 69 HOPTIONS, 71 HOPTIONS,CUTLENGTH,FPUMPTUNE2, 70 72 #ifdef COIN_HAS_CBC 71 73 LOGLEVEL , … … 75 77 PRIMALPIVOT,PRESOLVE,CRASH,BIASLU,PERTURBATION,MESSAGES,AUTOSCALE, 76 78 CHOLESKY,KKT,BARRIERSCALE,GAMMA,CROSSOVER,PFI,INTPRINT,VECTOR, 79 FACTORIZATION,ALLCOMMANDS, 77 80 78 81 NODESTRATEGY = 251,BRANCHSTRATEGY,CUTSSTRATEGY,HEURISTICSTRATEGY, … … 82 85 LANDPCUTS,RINS,RESIDCUTS,RENS,DIVINGS,DIVINGC,DIVINGF,DIVINGG,DIVINGL, 83 86 DIVINGP,DIVINGV,DINS,PIVOTANDFIX,RANDROUND,NAIVE,ZEROHALFCUTS,CPX, 87 CROSSOVER2,PIVOTANDCOMPLEMENT,VND, 84 88 85 89 DIRECTORY=301,DIRSAMPLE,DIRNETLIB,DIRMIPLIB,IMPORT,EXPORT,RESTORE,SAVE,DUALSIMPLEX,PRIMALSIMPLEX,EITHERSIMPLEX, … … 87 91 TIGHTEN,FAKEBOUND,HELP,PLUSMINUS,NETWORK,ALLSLACK,REVERSE,BARRIER,NETLIB_BARRIER,NETLIB_TUNE, 88 92 REALLY_SCALE,BASISIN,BASISOUT,SOLVECONTINUOUS,CLEARCUTS,VERSION,STATISTICS,DEBUG,DUMMY,PRINTMASK, 89 OUTDUPROWS,USERCLP,MODELIN,CSVSTATISTICS,STOREDFILE, 93 OUTDUPROWS,USERCLP,MODELIN,CSVSTATISTICS,STOREDFILE,ENVIRONMENT, 90 94 91 95 BAB=351,MIPLIB,STRENGTHEN,PRIORITYIN,USERCBC,DOHEURISTIC, … … 107 111 CbcOrClpParam ( ); 108 112 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); 110 114 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); 112 116 // Other strings will be added by insert 113 117 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); 115 119 // Action 116 120 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); 118 122 /// Copy constructor. 119 123 CbcOrClpParam(const CbcOrClpParam &); … … 178 182 /// Returns name which could match 179 183 std::string matchName ( ) const; 184 /// Returns length of name for ptinting 185 int lengthMatchName ( ) const; 180 186 /// Returns parameter option which matches (-1 if none) 181 187 int parameterOption ( std::string check ) const; … … 212 218 { return type_;} 213 219 /// whether to display 214 inline booldisplayThis() const220 inline int displayThis() const 215 221 { return display_;} 216 222 /// Set Long help … … 265 271 int currentKeyWord_; 266 272 /// Display on ? 267 booldisplay_;273 int display_; 268 274 /// Integer parameter - current value 269 275 int intValue_; -
stable/1.11/Clp/src/ClpCholeskyBase.cpp
r1321 r1458 1 /* $Id$ */ 1 2 // Copyright (C) 2002, International Business Machines 2 3 // 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 3 21 4 22 #include "CoinPragma.hpp" 5 23 6 #include <iostream> 24 #include <iostream> 7 25 8 26 #include "ClpCholeskyBase.hpp" … … 86 104 #else 87 105 // actually long double 88 workDouble_ = (double *) ClpCopyOfArray((longWork *) rhs.workDouble_,numberRows_);106 workDouble_ = reinterpret_cast<double *> (ClpCopyOfArray(reinterpret_cast<CoinWorkDouble *> (rhs.workDouble_),numberRows_)); 89 107 #endif 90 108 link_ = ClpCopyOfArray(rhs.link_,numberRows_); … … 171 189 #else 172 190 // actually long double 173 workDouble_ = (double *) ClpCopyOfArray((longWork *) rhs.workDouble_,numberRows_);191 workDouble_ = reinterpret_cast<double *> (ClpCopyOfArray(reinterpret_cast<CoinWorkDouble *> (rhs.workDouble_),numberRows_)); 174 192 #endif 175 193 link_ = ClpCopyOfArray(rhs.link_,numberRows_); … … 195 213 region1 is rows+columns, region2 is rows */ 196 214 void 197 ClpCholeskyBase::solveKKT ( double * region1, double * region2, const double * diagonal,198 double diagonalScaleFactor)215 ClpCholeskyBase::solveKKT (CoinWorkDouble * region1, CoinWorkDouble * region2, const CoinWorkDouble * diagonal, 216 CoinWorkDouble diagonalScaleFactor) 199 217 { 200 218 if (!doKKT_) { … … 202 220 int numberColumns = model_->numberColumns(); 203 221 int numberTotal = numberRows_+numberColumns; 204 double * region1Save = new double[numberTotal];222 CoinWorkDouble * region1Save = new CoinWorkDouble[numberTotal]; 205 223 for (iColumn=0;iColumn<numberTotal;iColumn++) { 206 224 region1[iColumn] *= diagonal[iColumn]; … … 209 227 multiplyAdd(region1+numberColumns,numberRows_,-1.0,region2,1.0); 210 228 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; 214 232 if (maximumRHS>1.0e-30) { 215 233 if (maximumRHS<=0.5) { 216 double factor=2.0;234 CoinWorkDouble factor=2.0; 217 235 while (maximumRHS<=0.5) { 218 236 maximumRHS*=factor; … … 220 238 } /* endwhile */ 221 239 } else if (maximumRHS>=2.0&&maximumRHS<=COIN_DBL_MAX) { 222 double factor=0.5;240 CoinWorkDouble factor=0.5; 223 241 while (maximumRHS>=2.0) { 224 242 maximumRHS*=factor; … … 246 264 int numberColumns = model_->numberColumns(); 247 265 int numberTotal = numberColumns + numberRowsModel; 248 double * array = new double [numberRows_];266 CoinWorkDouble * array = new CoinWorkDouble [numberRows_]; 249 267 CoinMemcpyN(region1,numberTotal,array); 250 268 CoinMemcpyN(region2,numberRowsModel,array+numberTotal); … … 253 271 int iRow; 254 272 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) { 256 274 printf("row region1 %d dropped %g\n",iRow,array[iRow]); 257 275 } 258 276 } 259 277 for (;iRow<numberRows_;iRow++) { 260 if (rowsDropped_[iRow]&& fabs(array[iRow])>1.0e-8) {278 if (rowsDropped_[iRow]&&CoinAbs(array[iRow])>1.0e-8) { 261 279 printf("row region2 %d dropped %g\n",iRow,array[iRow]); 262 280 } … … 621 639 { 622 640 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 623 654 int numberRowsModel = model_->numberRows(); 624 655 int numberColumns = model_->numberColumns(); … … 636 667 rowsDropped_ = new char [numberRows_]; 637 668 numberRowsDropped_=0; 669 memset(rowsDropped_,0,numberRows_); 638 670 rowCopy_ = model->clpMatrix()->reverseOrderedCopy(); 639 671 const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts(); … … 749 781 delete [] count; 750 782 permuteInverse_ = new int [numberRows_]; 751 memset(rowsDropped_,0,numberRows_);752 783 for (iRow=0;iRow<numberRows_;iRow++) { 753 784 //permute_[iRow]=iRow; // force no permute … … 757 788 } 758 789 return 0; 790 } 791 #if BASE_ORDER==1 792 /* Orders rows and saves pointer to matrix.and model */ 793 int 794 ClpCholeskyBase::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 1101 typedef int WSI; 1102 1103 #define NORDTHRESH 7 1104 #define MAXIW 2147000000 1105 1106 //#define WSSMP_DBG 1107 #ifdef WSSMP_DBG 1108 static 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 1116 static void 1117 myamlf(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 { 1121 WSI l, i, j, k; 1122 double x, y; 1123 WSI 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); 1153 L100:; 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 1175 chk(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 1191 chk(2,deg,n-1); 1192 #endif 1193 node = head[deg - 1]; 1194 #ifdef WSSMP_DBG 1195 chk(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 1216 chk(4,l,adjln); 1217 #endif 1218 tn = adjncy[l-1]; 1219 l++; 1220 #ifdef WSSMP_DBG 1221 chk(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 1228 chk(6,i0,adjln); 1229 #endif 1230 i = adjncy[i0 - 1]; 1231 i0++; 1232 #ifdef WSSMP_DBG 1233 chk(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 1249 chk(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 1258 chk(9,j6,adjln); 1259 #endif 1260 j = -adjncy[j6 - 1]; 1261 j6++; 1262 if (j > 0) { 1263 #ifdef WSSMP_DBG 1264 chk(10,j7,adjln); 1265 chk(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 1275 chk(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 1291 if (ipp) chk(13,ipp,n); 1292 if (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 1309 chk(15,j,adjln); 1310 #endif 1311 i = adjncy[j - 1]; 1312 #ifdef WSSMP_DBG 1313 chk(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 1320 chk(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 1327 if (ipp) chk(18,ipp,n); 1328 if (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 1338 chk(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 1348 chk(21,j3,adjln); 1349 #endif 1350 i = adjncy[j3 - 1]; 1351 #ifdef WSSMP_DBG 1352 chk(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 1360 chk(23,l,adjln); 1361 #endif 1362 tn = adjncy[l - 1]; 1363 #ifdef WSSMP_DBG 1364 chk(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 1374 chk(25,j3,adjln); 1375 #endif 1376 i = adjncy[j3 - 1]; 1377 i5 = deg = 0; 1378 #ifdef WSSMP_DBG 1379 chk(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 1384 chk(27,l,adjln); 1385 #endif 1386 tn = adjncy[l - 1]; 1387 #ifdef WSSMP_DBG 1388 chk(70,tn,n); 1389 #endif 1390 j8 = flag[tn - 1]; 1391 if (j8) { 1392 deg += (j8 - fltag); 1393 #ifdef WSSMP_DBG 1394 chk(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 1406 chk(29,l,adjln); 1407 #endif 1408 j = adjncy[l - 1]; 1409 #ifdef WSSMP_DBG 1410 chk(30,j,n); 1411 #endif 1412 numii = varbl[j - 1]; 1413 if (numii > 0) { 1414 deg += numii; 1415 #ifdef WSSMP_DBG 1416 chk(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 1434 chk(32,j2,adjln); 1435 chk(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 1443 chk(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 1457 chk(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 1467 chk(37,j3,adjln); 1468 #endif 1469 i = adjncy[j3 - 1]; 1470 #ifdef WSSMP_DBG 1471 chk(38,i,n); 1472 #endif 1473 if (varbl[i - 1] < 0) { 1474 i5 = perm[i - 1]; 1475 #ifdef WSSMP_DBG 1476 chk(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 1484 chk(40,j,n); 1485 #endif 1486 i = perm[j - 1]; 1487 perm[j - 1] = 0; 1488 } 1489 while (i) { 1490 #ifdef WSSMP_DBG 1491 chk(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 1500 chk(42,l,adjln); 1501 chk(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 1509 chk(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 1516 chk(45,l,adjln); 1517 chk(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 1529 chk(47,j,n); 1530 chk(48,jpp,n); 1531 #endif 1532 j = snxt[j - 1]; 1533 snxt[jpp - 1] = j; 1534 } else { 1535 jpp = j; 1536 #ifdef WSSMP_DBG 1537 chk(49,j,n); 1538 #endif 1539 j = snxt[j - 1]; 1540 } 1541 } else { 1542 jpp = j; 1543 #ifdef WSSMP_DBG 1544 chk(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 1554 chk(51,i,n); 1555 #endif 1556 i = snxt[i - 1]; 1557 } 1558 } 1559 } 1560 } 1561 } 1562 j8 = nm1 - counter; 1563 switch (speed) { 1564 case 1: 1565 for (j = j3 = j4; j3 < j5; j3++) { 1566 #ifdef WSSMP_DBG 1567 chk(52,j3,adjln); 1568 #endif 1569 i = adjncy[j3 - 1]; 1570 #ifdef WSSMP_DBG 1571 chk(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 1579 chk(54,l,adjln); 1580 chk(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 1603 chk(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 1610 chk(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; 1618 case 2: 1619 for (j = j3 = j4; j3 < j5; j3++) { 1620 #ifdef WSSMP_DBG 1621 chk(58,j3,adjln); 1622 #endif 1623 i = adjncy[j3 - 1]; 1624 #ifdef WSSMP_DBG 1625 chk(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 1647 chk(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 1654 chk(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; 1662 default: 1663 for (j = j3 = j4; j3 < j5; j3++) { 1664 #ifdef WSSMP_DBG 1665 chk(62,j3,adjln); 1666 #endif 1667 i = adjncy[j3 - 1]; 1668 #ifdef WSSMP_DBG 1669 chk(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 1679 chk(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 1686 chk(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 1697 chk(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 1710 chk(67,tn,n); 1711 #endif 1712 k = -invp[tn - 1]; 1713 i = l; 1714 #ifdef WSSMP_DBG 1715 chk(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 1729 chk(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--------------------------------------------------------------------------*/ 1739 void 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 { 1743 WSI 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 1819 int 1820 ClpCholeskyBase::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; 759 1918 } 760 1919 /* Does Symbolic factorization given permutation. … … 1237 2396 #else 1238 2397 // actually long double 1239 workDouble_ = (double *) (new longWork[numberRows_]);2398 workDouble_ = reinterpret_cast<double *> (new CoinWorkDouble[numberRows_]); 1240 2399 #endif 1241 2400 diagonal_ = new longDouble[numberRows_]; … … 1338 2497 int k=iRow; 1339 2498 int linked = link_[iRow]; 2499 #ifndef NDEBUG 2500 int count=0; 2501 #endif 1340 2502 while (linked<=kRow) { 1341 2503 k=linked; 1342 2504 linked = link_[k]; 2505 #ifndef NDEBUG 2506 count++; 2507 assert (count<numberRows_); 2508 #endif 1343 2509 } 1344 2510 nz++; … … 1521 2687 /* Factorize - filling in rowsDropped and returning number dropped */ 1522 2688 int 1523 ClpCholeskyBase::factorize(const double * diagonal, int * rowsDropped)2689 ClpCholeskyBase::factorize(const CoinWorkDouble * diagonal, int * rowsDropped) 1524 2690 { 1525 2691 const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts(); … … 1533 2699 int numberColumns=model_->clpMatrix()->getNumCols(); 1534 2700 //perturbation 1535 longDouble perturbation=model_->diagonalPerturbation()*model_->diagonalNorm();2701 CoinWorkDouble perturbation=model_->diagonalPerturbation()*model_->diagonalNorm(); 1536 2702 //perturbation=perturbation*perturbation*100000000.0; 1537 2703 if (perturbation>1.0) { … … 1540 2706 std::cout <<"large perturbation "<<perturbation<<std::endl; 1541 2707 #endif 1542 perturbation= sqrt(perturbation);;2708 perturbation=CoinSqrt(perturbation); 1543 2709 perturbation=1.0; 1544 2710 } … … 1548 2714 CoinZeroN(work,numberRows_); 1549 2715 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; 1552 2718 int numberDense=0; 1553 2719 if (!doKKT_) { 1554 const double * diagonalSlack = diagonal + numberColumns;2720 const CoinWorkDouble * diagonalSlack = diagonal + numberColumns; 1555 2721 if (dense_) 1556 2722 numberDense=dense_->numberRows(); … … 1577 2743 } 1578 2744 } 1579 longDouble delta2 = model_->delta(); // add delta*delta to diagonal2745 CoinWorkDouble delta2 = model_->delta(); // add delta*delta to diagonal 1580 2746 delta2 *= delta2; 1581 2747 // largest in initial matrix 1582 double largest2=1.0e-20;2748 CoinWorkDouble largest2=1.0e-20; 1583 2749 for (iRow=0;iRow<numberRows_;iRow++) { 1584 2750 longDouble * put = sparseFactor_+choleskyStart_[iRow]; … … 1597 2763 CoinBigIndex start=columnStart[iColumn]; 1598 2764 CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn]; 1599 longDouble multiplier = diagonal[iColumn]*elementByRow[k];2765 CoinWorkDouble multiplier = diagonal[iColumn]*elementByRow[k]; 1600 2766 for (CoinBigIndex j=start;j<end;j++) { 1601 2767 int jRow=row[j]; 1602 2768 int jNewRow = permuteInverse_[jRow]; 1603 2769 if (jNewRow>=iRow&&!rowsDropped_[jRow]) { 1604 longDouble value=element[j]*multiplier;2770 CoinWorkDouble value=element[j]*multiplier; 1605 2771 work[jNewRow] += value; 1606 2772 } … … 1609 2775 } 1610 2776 diagonal_[iRow]=work[iRow]; 1611 largest2 = CoinMax(largest2, fabs(work[iRow]));2777 largest2 = CoinMax(largest2,CoinAbs(work[iRow])); 1612 2778 work[iRow]=0.0; 1613 2779 int j; … … 1615 2781 int jRow = which[j]; 1616 2782 put[j]=work[jRow]; 1617 largest2 = CoinMax(largest2, fabs(work[jRow]));2783 largest2 = CoinMax(largest2,CoinAbs(work[jRow])); 1618 2784 work[jRow]=0.0; 1619 2785 } … … 1629 2795 //check sizes 1630 2796 largest2*=1.0e-20; 1631 largest = CoinMin(largest2, 1.0e-11);2797 largest = CoinMin(largest2,CHOL_SMALL_VALUE); 1632 2798 int numberDroppedBefore=0; 1633 2799 for (iRow=0;iRow<numberRows_;iRow++) { … … 1636 2802 rowsDropped[iRow]=dropped; 1637 2803 if (!dropped) { 1638 longDouble diagonal = diagonal_[iRow];2804 CoinWorkDouble diagonal = diagonal_[iRow]; 1639 2805 if (diagonal>largest2) { 1640 2806 diagonal_[iRow]=diagonal+perturbation; … … 1663 2829 for ( i=0;i<numberRows_;i++) { 1664 2830 assert (diagonal_[i]>=0.0); 1665 diagonal_[i]= sqrt(diagonal_[i]);2831 diagonal_[i]=CoinSqrt(diagonal_[i]); 1666 2832 } 1667 2833 // Update dense columns (just L) … … 1673 2839 a[j]=0.0; 1674 2840 } 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 } 1676 2858 } 1677 2859 dense_->resetRowsDropped(); … … 1682 2864 const longDouble * a = denseColumn_+i*numberRows_; 1683 2865 // do diagonal 1684 longDouble value = denseDiagonal[i];2866 CoinWorkDouble value = denseDiagonal[i]; 1685 2867 const longDouble * b = denseColumn_+i*numberRows_; 1686 2868 for (int k=0;k<numberRows_;k++) … … 1688 2870 denseDiagonal[i]=value; 1689 2871 for (int j=i+1;j<numberDense;j++) { 1690 longDouble value = 0.0;2872 CoinWorkDouble value = 0.0; 1691 2873 const longDouble * b = denseColumn_+j*numberRows_; 1692 2874 for (int k=0;k<numberRows_;k++) … … 1778 2960 } 1779 2961 if (permuted) { 1780 longDouble delta2 = model_->delta(); // add delta*delta to bottom2962 CoinWorkDouble delta2 = model_->delta(); // add delta*delta to bottom 1781 2963 delta2 *= delta2; 1782 2964 // Need to permute - ugly … … 1788 2970 if (iOriginalRow<numberColumns) { 1789 2971 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) { 1792 2974 value = 1.0/value; 1793 largest = CoinMax(largest, fabs(value));2975 largest = CoinMax(largest,CoinAbs(value)); 1794 2976 diagonal_[iRow] = -value; 1795 2977 CoinBigIndex start=columnStart[iColumn]; … … 1800 2982 if (kRow>iRow) { 1801 2983 work[kRow]=element[j]; 1802 largest = CoinMax(largest, fabs(element[j]));2984 largest = CoinMax(largest,CoinAbs(element[j])); 1803 2985 } 1804 2986 } … … 1807 2989 } 1808 2990 } 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) { 1811 2993 value = 1.0/value; 1812 largest = CoinMax(largest, fabs(value));2994 largest = CoinMax(largest,CoinAbs(value)); 1813 2995 } else { 1814 2996 value = 1.0e100; … … 1828 3010 if (jNewRow>iRow) { 1829 3011 work[jNewRow]=elementByRow[j]; 1830 largest = CoinMax(largest, fabs(elementByRow[j]));3012 largest = CoinMax(largest,CoinAbs(elementByRow[j])); 1831 3013 } 1832 3014 } … … 1857 3039 CoinBigIndex j; 1858 3040 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) { 1861 3043 value = 1.0/value; 1862 3044 for (j=columnQuadraticStart[iColumn]; … … 1870 3052 } 1871 3053 } 1872 largest = CoinMax(largest, fabs(value));3054 largest = CoinMax(largest,CoinAbs(value)); 1873 3055 diagonal_[iRow] = -value; 1874 3056 CoinBigIndex start=columnStart[iColumn]; … … 1879 3061 if (kRow>iRow) { 1880 3062 work[kRow]=element[j]; 1881 largest = CoinMax(largest, fabs(element[j]));3063 largest = CoinMax(largest,CoinAbs(element[j])); 1882 3064 } 1883 3065 } … … 1886 3068 } 1887 3069 } 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) { 1890 3072 value = 1.0/value; 1891 largest = CoinMax(largest, fabs(value));3073 largest = CoinMax(largest,CoinAbs(value)); 1892 3074 } else { 1893 3075 value = 1.0e100; … … 1907 3089 if (jNewRow>iRow) { 1908 3090 work[jNewRow]=elementByRow[j]; 1909 largest = CoinMax(largest, fabs(elementByRow[j]));3091 largest = CoinMax(largest,CoinAbs(elementByRow[j])); 1910 3092 } 1911 3093 } … … 1931 3113 longDouble * put = sparseFactor_+choleskyStart_[iColumn]; 1932 3114 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) { 1935 3117 value = 1.0/value; 1936 largest = CoinMax(largest, fabs(value));3118 largest = CoinMax(largest,CoinAbs(value)); 1937 3119 diagonal_[iColumn] = -value; 1938 3120 CoinBigIndex start=columnStart[iColumn]; … … 1942 3124 //sparseFactor_[numberElements++]=element[j]; 1943 3125 work[row[j]+numberTotal]=element[j]; 1944 largest = CoinMax(largest, fabs(element[j]));3126 largest = CoinMax(largest,CoinAbs(element[j])); 1945 3127 } 1946 3128 } else { … … 1965 3147 int * which = choleskyRow_+indexStart_[iColumn]; 1966 3148 int number = choleskyStart_[iColumn+1]-choleskyStart_[iColumn]; 1967 longDouble value = diagonal[iColumn];3149 CoinWorkDouble value = diagonal[iColumn]; 1968 3150 CoinBigIndex j; 1969 if ( fabs(value)>1.0e-100) {3151 if (CoinAbs(value)>1.0e-100) { 1970 3152 value = 1.0/value; 1971 3153 for (j=columnQuadraticStart[iColumn]; … … 1978 3160 } 1979 3161 } 1980 largest = CoinMax(largest, fabs(value));3162 largest = CoinMax(largest,CoinAbs(value)); 1981 3163 diagonal_[iColumn] = -value; 1982 3164 CoinBigIndex start=columnStart[iColumn]; … … 1984 3166 for (j=start;j<end;j++) { 1985 3167 work[row[j]+numberTotal]=element[j]; 1986 largest = CoinMax(largest, fabs(element[j]));3168 largest = CoinMax(largest,CoinAbs(element[j])); 1987 3169 } 1988 3170 for (j=0;j<number;j++) { … … 2005 3187 longDouble * put = sparseFactor_+choleskyStart_[iColumn]; 2006 3188 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) { 2009 3191 value = 1.0/value; 2010 largest = CoinMax(largest, fabs(value));3192 largest = CoinMax(largest,CoinAbs(value)); 2011 3193 } else { 2012 3194 value = 1.0e100; … … 2023 3205 } 2024 3206 // Finish diagonal 2025 longDouble delta2 = model_->delta(); // add delta*delta to bottom3207 CoinWorkDouble delta2 = model_->delta(); // add delta*delta to bottom 2026 3208 delta2 *= delta2; 2027 3209 for (iRow=0;iRow<numberRowsModel;iRow++) { … … 2037 3219 //check sizes 2038 3220 largest*=1.0e-20; 2039 largest = CoinMin(largest, 1.0e-11);3221 largest = CoinMin(largest,CHOL_SMALL_VALUE); 2040 3222 doubleParameters_[10]=CoinMax(1.0e-20,largest); 2041 3223 integerParameters_[20]=0; … … 2056 3238 // Should save adjustments in ..R_ 2057 3239 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(); 2060 3242 for (iRow=0;iRow<numberTotal;iRow++) { 2061 3243 if (rowsDropped2[iRow]) { … … 2093 3275 ClpCholeskyBase::factorizePart2(int * rowsDropped) 2094 3276 { 2095 longDouble largest=doubleParameters_[3];2096 longDouble smallest=doubleParameters_[4];3277 CoinWorkDouble largest=doubleParameters_[3]; 3278 CoinWorkDouble smallest=doubleParameters_[4]; 2097 3279 int numberDropped=integerParameters_[20]; 2098 3280 // probably done before … … 2159 3341 for ( jRow=lastRow;jRow<iRow;jRow++) { 2160 3342 int jCount = jRow-lastRow; 2161 longDouble diagonalValue = diagonal_[jRow];3343 CoinWorkDouble diagonalValue = diagonal_[jRow]; 2162 3344 CoinBigIndex start=choleskyStart_[jRow]; 2163 3345 CoinBigIndex end=choleskyStart_[jRow+1]; … … 2165 3347 jCount--; 2166 3348 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; 2169 3351 diagonalValue -= a_jk*value1; 2170 3352 for (CoinBigIndex j=start;j<end;j++) … … 2222 3404 } 2223 3405 // for each column L[*,kRow] that affects L[*,iRow] 2224 longDouble diagonalValue=diagonal_[iRow];3406 CoinWorkDouble diagonalValue=diagonal_[iRow]; 2225 3407 int nextRow = link_[iRow]; 2226 3408 int kRow=0; … … 2234 3416 CoinBigIndex end=choleskyStart_[kRow+1]; 2235 3417 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; 2238 3420 // update first 2239 3421 first[kRow]=k; … … 2259 3441 CoinBigIndex j=first[kkRow]; 2260 3442 //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; 2263 3445 diagonalValue -= a*dValue; 2264 3446 work[kkRow]=dValue; … … 2271 3453 for (int i=0;i<length;i++) { 2272 3454 int lRow = choleskyRow_[currentIndex++]; 2273 longDouble t0 = work[lRow];3455 CoinWorkDouble t0 = work[lRow]; 2274 3456 for (int kkRow=kRow;kkRow<last;kkRow++) { 2275 3457 CoinBigIndex j = first[kkRow]+i; … … 2340 3522 for (CoinBigIndex j=start;j<end;j++) { 2341 3523 int jRow = choleskyRow_[j+offset]; 2342 longDouble value = sparseFactor_[j]-work[jRow];3524 CoinWorkDouble value = sparseFactor_[j]-work[jRow]; 2343 3525 work[jRow]=0.0; 2344 3526 sparseFactor_[j]= diagonalValue*value; … … 2350 3532 // do dense 2351 3533 // update dense part 2352 updateDense(d, work,first);3534 updateDense(d,/*work,*/first); 2353 3535 ClpCholeskyDense dense; 2354 3536 // just borrow space … … 2371 3553 dense.setIntegerParameter(20,0); 2372 3554 dense.setIntegerParameter(34,firstPositive); 3555 dense.setModel(model_); 2373 3556 dense.factorizePart2(dropped); 2374 3557 largest=dense.getDoubleParameter(3); … … 2387 3570 } 2388 3571 // Updates dense part (broken out for profiling) 2389 void ClpCholeskyBase::updateDense(longDouble * d, longDouble * work,int * first)3572 void ClpCholeskyBase::updateDense(longDouble * d, /*longDouble * work,*/ int * first) 2390 3573 { 2391 3574 for (int iRow=0;iRow<firstDense_;iRow++) { … … 2395 3578 CoinBigIndex offset = indexStart_[iRow]-choleskyStart_[iRow]; 2396 3579 if (clique_[iRow]<2) { 2397 longDouble dValue=d[iRow];3580 CoinWorkDouble dValue=d[iRow]; 2398 3581 for (CoinBigIndex k=start;k<end;k++) { 2399 3582 int kRow = choleskyRow_[k+offset]; 2400 3583 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; 2403 3586 diagonal_[kRow] -= value1*a_ik; 2404 3587 CoinBigIndex base = choleskyStart_[kRow]-kRow-1; 2405 3588 for (CoinBigIndex j=k+1;j<end;j++) { 2406 3589 int jRow=choleskyRow_[j+offset]; 2407 longDouble a_jk = sparseFactor_[j];3590 CoinWorkDouble a_jk = sparseFactor_[j]; 2408 3591 sparseFactor_[base+jRow] -= a_jk*value1; 2409 3592 } … … 2411 3594 } else if (clique_[iRow]<3) { 2412 3595 // 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]; 2415 3598 int offset1 = first[iRow+1]-start; 2416 3599 // skip row … … 2419 3602 int kRow = choleskyRow_[k+offset]; 2420 3603 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; 2425 3608 diagonal_[kRow] -= value0*a_ik0 + value1*a_ik1; 2426 3609 CoinBigIndex base = choleskyStart_[kRow]-kRow-1; 2427 3610 for (CoinBigIndex j=k+1;j<end;j++) { 2428 3611 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]; 2431 3614 sparseFactor_[base+jRow] -= a_jk0*value0+a_jk1*value1; 2432 3615 } … … 2441 3624 // maybe later get fancy on big cliques and do transpose copy 2442 3625 // 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]; 2446 3629 // get offsets and skip rows 2447 3630 int offset1 = first[++iRow]-start; … … 2450 3633 int kRow = choleskyRow_[k+offset]; 2451 3634 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; 2459 3642 CoinBigIndex base = choleskyStart_[kRow]-kRow-1; 2460 3643 diagonal_[kRow] = diagonalValue - value0*a_ik0 - value1*a_ik1 - value2*a_ik2; 2461 3644 for (CoinBigIndex j=k+1;j<end;j++) { 2462 3645 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]; 2466 3649 sparseFactor_[base+jRow] -= a_jk0*value0+a_jk1*value1+a_jk2*value2; 2467 3650 } … … 2472 3655 // maybe later get fancy on big cliques and do transpose copy 2473 3656 // 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]; 2478 3661 // get offsets and skip rows 2479 3662 int offset1 = first[++iRow]-start; … … 2483 3666 int kRow = choleskyRow_[k+offset]; 2484 3667 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; 2494 3677 CoinBigIndex base = choleskyStart_[kRow]-kRow-1; 2495 3678 diagonal_[kRow] = diagonalValue - (value0*a_ik0 + value1*a_ik1 + value2*a_ik2+value3*a_ik3); 2496 3679 for (CoinBigIndex j=k+1;j<end;j++) { 2497 3680 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]; 2502 3685 sparseFactor_[base+jRow] -= a_jk0*value0+a_jk1*value1+a_jk2*value2+a_jk3*value3; 2503 3686 } … … 2510 3693 /* Uses factorization to solve. */ 2511 3694 void 2512 ClpCholeskyBase::solve ( double * region)3695 ClpCholeskyBase::solve (CoinWorkDouble * region) 2513 3696 { 2514 3697 if (!whichDense_) { … … 2520 3703 // do change; 2521 3704 int numberDense = dense_->numberRows(); 2522 double * change = new double[numberDense];3705 CoinWorkDouble * change = new CoinWorkDouble[numberDense]; 2523 3706 for (i=0;i<numberDense;i++) { 2524 3707 const longDouble * a = denseColumn_+i*numberRows_; 2525 longDouble value =0.0;3708 CoinWorkDouble value =0.0; 2526 3709 for (int iRow=0;iRow<numberRows_;iRow++) 2527 3710 value += a[iRow]*region[iRow]; … … 2532 3715 for (i=0;i<numberDense;i++) { 2533 3716 const longDouble * a = denseColumn_+i*numberRows_; 2534 longDouble value = change[i];3717 CoinWorkDouble value = change[i]; 2535 3718 for (int iRow=0;iRow<numberRows_;iRow++) 2536 3719 region[iRow] -= value*a[iRow]; … … 2545 3728 */ 2546 3729 void 2547 ClpCholeskyBase::solve( double * region, int type)3730 ClpCholeskyBase::solve(CoinWorkDouble * region, int type) 2548 3731 { 2549 3732 #ifdef CLP_DEBUG 2550 3733 double * regionX=NULL; 2551 if (sizeof( longWork)!=sizeof(double)&&type==3) {3734 if (sizeof(CoinWorkDouble)!=sizeof(double)&&type==3) { 2552 3735 regionX=ClpCopyOfArray(region,numberRows_); 2553 3736 } 2554 3737 #endif 2555 longWork * work = reinterpret_cast<longWork*> (workDouble_);3738 CoinWorkDouble * work = reinterpret_cast<CoinWorkDouble *> (workDouble_); 2556 3739 int i; 2557 3740 CoinBigIndex j; … … 2563 3746 case 1: 2564 3747 for (i=0;i<numberRows_;i++) { 2565 longDouble value=work[i];3748 CoinWorkDouble value=work[i]; 2566 3749 CoinBigIndex offset = indexStart_[i]-choleskyStart_[i]; 2567 3750 for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) { … … 2578 3761 for (i=numberRows_-1;i>=0;i--) { 2579 3762 CoinBigIndex offset = indexStart_[i]-choleskyStart_[i]; 2580 longDouble value=work[i]*diagonal_[i];3763 CoinWorkDouble value=work[i]*diagonal_[i]; 2581 3764 for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) { 2582 3765 int iRow = choleskyRow_[j+offset]; … … 2591 3774 for (i=0;i<firstDense_;i++) { 2592 3775 CoinBigIndex offset = indexStart_[i]-choleskyStart_[i]; 2593 longDouble value=work[i];3776 CoinWorkDouble value=work[i]; 2594 3777 for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) { 2595 3778 int iRow = choleskyRow_[j+offset]; … … 2603 3786 int nDense = numberRows_-firstDense_; 2604 3787 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_); 2610 3789 for (i=numberRows_-1;i>=firstDense_;i--) { 2611 longDouble value=work[i];3790 CoinWorkDouble value=work[i]; 2612 3791 int iRow = permute_[i]; 2613 3792 region[iRow]=value; … … 2616 3795 for (i=firstDense_-1;i>=0;i--) { 2617 3796 CoinBigIndex offset = indexStart_[i]-choleskyStart_[i]; 2618 longDouble value=work[i]*diagonal_[i];3797 CoinWorkDouble value=work[i]*diagonal_[i]; 2619 3798 for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) { 2620 3799 int iRow = choleskyRow_[j+offset]; … … 2634 3813 double largestO=0.0; 2635 3814 for (i=0;i<numberRows_;i++) { 2636 largestO = CoinMax(largestO, fabs(regionX[i]));3815 largestO = CoinMax(largestO,CoinAbs(regionX[i])); 2637 3816 } 2638 3817 for (i=0;i<numberRows_;i++) { … … 2642 3821 for (i=0;i<firstDense_;i++) { 2643 3822 CoinBigIndex offset = indexStart_[i]-choleskyStart_[i]; 2644 longDouble value=work[i];3823 CoinWorkDouble value=work[i]; 2645 3824 for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) { 2646 3825 int iRow = choleskyRow_[j+offset]; … … 2654 3833 int nDense = numberRows_-firstDense_; 2655 3834 dense.reserveSpace(this,nDense); 2656 dense.solve Long(work+firstDense_);3835 dense.solve(work+firstDense_); 2657 3836 for (i=numberRows_-1;i>=firstDense_;i--) { 2658 longDouble value=work[i];3837 CoinWorkDouble value=work[i]; 2659 3838 int iRow = permute_[i]; 2660 3839 regionX[iRow]=value; … … 2663 3842 for (i=firstDense_-1;i>=0;i--) { 2664 3843 CoinBigIndex offset = indexStart_[i]-choleskyStart_[i]; 2665 longDouble value=work[i]*diagonal_[i];3844 CoinWorkDouble value=work[i]*diagonal_[i]; 2666 3845 for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) { 2667 3846 int iRow = choleskyRow_[j+offset]; … … 2675 3854 double largestV=0.0; 2676 3855 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])); 2679 3858 } 2680 3859 printf("largest difference %g, largest %g, largest original %g\n", … … 2684 3863 #endif 2685 3864 } 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. */ 2689 3867 void 2690 ClpCholeskyBase::solve Long(longDouble * region, int type)3868 ClpCholeskyBase::solve (CoinWorkDouble * region) 2691 3869 { 3870 assert (!whichDense_) ; 3871 CoinWorkDouble * work = reinterpret_cast<CoinWorkDouble *> (workDouble_); 2692 3872 int i; 2693 3873 CoinBigIndex j; 2694 3874 for (i=0;i<numberRows_;i++) { 2695 3875 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]; 2722 3895 int iRow = permute_[i]; 2723 3896 region[iRow]=value; 2724 3897 } 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; 2760 3909 } 2761 3910 } 2762 3911 #endif -
stable/1.11/Clp/src/ClpCholeskyBase.hpp
r1055 r1458 1 /* $Id$ */ 1 2 // Copyright (C) 2003, International Business Machines 2 3 // Corporation and others. All Rights Reserved. … … 6 7 #include "CoinPragma.hpp" 7 8 #include "CoinFinite.hpp" 9 //#define CLP_LONG_CHOLESKY 0 10 #ifndef CLP_LONG_CHOLESKY 8 11 #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 9 31 #if CLP_LONG_CHOLESKY>1 10 32 typedef long double longDouble; 11 typedef long double longWork; 33 #define CHOL_SMALL_VALUE 1.0e-15 12 34 #elif CLP_LONG_CHOLESKY==1 13 35 typedef double longDouble; 14 typedef long double longWork; 36 #define CHOL_SMALL_VALUE 1.0e-11 15 37 #else 16 38 typedef double longDouble; 17 typedef double longWork; 39 #define CHOL_SMALL_VALUE 1.0e-11 18 40 #endif 19 41 class ClpInterior; … … 46 68 /** Factorize - filling in rowsDropped and returning number dropped. 47 69 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) ; 49 71 /** Uses factorization to solve. */ 50 virtual void solve ( double * region) ;72 virtual void solve (CoinWorkDouble * region) ; 51 73 /** Uses factorization to solve. - given as if KKT. 52 74 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); 77 private: 78 /// AMD ordering 79 int orderAMD(); 80 public: 55 81 //@} 56 82 … … 141 167 protected: 142 168 /// 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;} 144 173 //@} 145 174 … … 162 191 If 1 and 2 then diagonal has sqrt of inverse otherwise inverse 163 192 */ 164 void solve(double * region, int type); 165 void solveLong(longDouble * region, int type); 193 void solve(CoinWorkDouble * region, int type); 166 194 /// Forms ADAT - returns nonzero if not enough memory 167 195 int preOrder(bool lowerTriangular, bool includeDiagonal, bool doKKT); 168 196 /// 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); 170 198 //@} 171 199 -
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. */ 6 4 #include "CoinPragma.hpp" 7 5 #include "CoinHelperFunctions.hpp" 6 #include "ClpHelperFunctions.hpp" 8 7 9 8 #include "ClpInterior.hpp" … … 12 11 #include "ClpQuadraticObjective.hpp" 13 12 14 / /#############################################################################15 / / Constructors / Destructor / Assignment16 / /#############################################################################17 18 / /-------------------------------------------------------------------19 / / Default Constructor20 / /-------------------------------------------------------------------13 /*#############################################################################*/ 14 /* Constructors / Destructor / Assignment*/ 15 /*#############################################################################*/ 16 17 /*-------------------------------------------------------------------*/ 18 /* Default Constructor */ 19 /*-------------------------------------------------------------------*/ 21 20 ClpCholeskyDense::ClpCholeskyDense () 22 21 : ClpCholeskyBase(), … … 26 25 } 27 26 28 / /-------------------------------------------------------------------29 / / Copy constructor30 / /-------------------------------------------------------------------27 /*-------------------------------------------------------------------*/ 28 /* Copy constructor */ 29 /*-------------------------------------------------------------------*/ 31 30 ClpCholeskyDense::ClpCholeskyDense (const ClpCholeskyDense & rhs) 32 31 : ClpCholeskyBase(rhs), 33 32 borrowSpace_(rhs.borrowSpace_) 34 33 { 35 assert(!rhs.borrowSpace_||!rhs.sizeFactor_); / / can't do if borrowing space36 } 37 38 39 / /-------------------------------------------------------------------40 / / Destructor41 / /-------------------------------------------------------------------34 assert(!rhs.borrowSpace_||!rhs.sizeFactor_); /* can't do if borrowing space*/ 35 } 36 37 38 /*-------------------------------------------------------------------*/ 39 /* Destructor */ 40 /*-------------------------------------------------------------------*/ 42 41 ClpCholeskyDense::~ClpCholeskyDense () 43 42 { 44 43 if (borrowSpace_) { 45 / / set NULL44 /* set NULL*/ 46 45 sparseFactor_=NULL; 47 46 workDouble_=NULL; … … 50 49 } 51 50 52 / /----------------------------------------------------------------53 / / Assignment operator54 / /-------------------------------------------------------------------51 /*----------------------------------------------------------------*/ 52 /* Assignment operator */ 53 /*-------------------------------------------------------------------*/ 55 54 ClpCholeskyDense & 56 55 ClpCholeskyDense::operator=(const ClpCholeskyDense& rhs) 57 56 { 58 57 if (this != &rhs) { 59 assert(!rhs.borrowSpace_||!rhs.sizeFactor_); / / can't do if borrowing space58 assert(!rhs.borrowSpace_||!rhs.sizeFactor_); /* can't do if borrowing space*/ 60 59 ClpCholeskyBase::operator=(rhs); 61 60 borrowSpace_=rhs.borrowSpace_; … … 63 62 return *this; 64 63 } 65 / /-------------------------------------------------------------------66 / / Clone67 / /-------------------------------------------------------------------64 /*-------------------------------------------------------------------*/ 65 /* Clone*/ 66 /*-------------------------------------------------------------------*/ 68 67 ClpCholeskyBase * ClpCholeskyDense::clone() const 69 68 { 70 69 return new ClpCholeskyDense(*this); 71 70 } 72 / / If not power of 2 then need to redo a bit71 /* If not power of 2 then need to redo a bit*/ 73 72 #define BLOCK 16 74 73 #define BLOCKSHIFT 4 75 / / Block unroll if power of 2 and at least 874 /* Block unroll if power of 2 and at least 8*/ 76 75 #define BLOCKUNROLL 77 76 … … 87 86 numberRows_ = numberRows; 88 87 int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT; 89 / / allow one stripe extra88 /* allow one stripe extra*/ 90 89 numberBlocks = numberBlocks + ((numberBlocks*(numberBlocks+1))/2); 91 90 sizeFactor_=numberBlocks*BLOCKSQ; 92 / /#define CHOL_COMPARE91 /*#define CHOL_COMPARE*/ 93 92 #ifdef CHOL_COMPARE 94 93 sizeFactor_ += 95000; … … 115 114 { 116 115 int numberBlocks = (numberRows+BLOCK-1)>>BLOCKSHIFT; 117 // allow one stripe extra116 /* allow one stripe extra*/ 118 117 numberBlocks = numberBlocks + ((numberBlocks*(numberBlocks+1))/2); 119 118 CoinBigIndex sizeFactor=numberBlocks*BLOCKSQ; … … 151 150 /* Factorize - filling in rowsDropped and returning number dropped */ 152 151 int 153 ClpCholeskyDense::factorize(const double * diagonal, int * rowsDropped)152 ClpCholeskyDense::factorize(const CoinWorkDouble * diagonal, int * rowsDropped) 154 153 { 155 154 assert (!borrowSpace_); … … 164 163 int numberColumns=model_->clpMatrix()->getNumCols(); 165 164 CoinZeroN(sparseFactor_,sizeFactor_); 166 / /perturbation167 longDouble perturbation=model_->diagonalPerturbation()*model_->diagonalNorm();165 /*perturbation*/ 166 CoinWorkDouble perturbation=model_->diagonalPerturbation()*model_->diagonalNorm(); 168 167 perturbation=perturbation*perturbation; 169 168 if (perturbation>1.0) { 170 169 #ifdef COIN_DEVELOP 171 / /if (model_->model()->logLevel()&4)170 /*if (model_->model()->logLevel()&4) */ 172 171 std::cout <<"large perturbation "<<perturbation<<std::endl; 173 172 #endif 174 perturbation= sqrt(perturbation);;173 perturbation=CoinSqrt(perturbation);; 175 174 perturbation=1.0; 176 175 } 177 176 int iRow; 178 177 int newDropped=0; 179 double largest=1.0;180 double smallest=COIN_DBL_MAX;181 longDouble delta2 = model_->delta(); // add delta*delta to diagonal178 CoinWorkDouble largest=1.0; 179 CoinWorkDouble smallest=COIN_DBL_MAX; 180 CoinWorkDouble delta2 = model_->delta(); /* add delta*delta to diagonal*/ 182 181 delta2 *= delta2; 183 182 if (!doKKT_) { 184 183 longDouble * work = sparseFactor_; 185 work--; / / skip diagonal184 work--; /* skip diagonal*/ 186 185 int addOffset=numberRows_-1; 187 const double * diagonalSlack = diagonal + numberColumns;188 / / largest in initial matrix189 double largest2=1.0e-20;186 const CoinWorkDouble * diagonalSlack = diagonal + numberColumns; 187 /* largest in initial matrix*/ 188 CoinWorkDouble largest2=1.0e-20; 190 189 for (iRow=0;iRow<numberRows_;iRow++) { 191 190 if (!rowsDropped_[iRow]) { 192 191 CoinBigIndex startRow=rowStart[iRow]; 193 192 CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow]; 194 longDouble diagonalValue = diagonalSlack[iRow]+delta2;193 CoinWorkDouble diagonalValue = diagonalSlack[iRow]+delta2; 195 194 for (CoinBigIndex k=startRow;k<endRow;k++) { 196 195 int iColumn=column[k]; 197 196 CoinBigIndex start=columnStart[iColumn]; 198 197 CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn]; 199 longDouble multiplier = diagonal[iColumn]*elementByRow[k];198 CoinWorkDouble multiplier = diagonal[iColumn]*elementByRow[k]; 200 199 for (CoinBigIndex j=start;j<end;j++) { 201 200 int jRow=row[j]; 202 201 if (!rowsDropped_[jRow]) { 203 202 if (jRow>iRow) { 204 longDouble value=element[j]*multiplier;203 CoinWorkDouble value=element[j]*multiplier; 205 204 work[jRow] += value; 206 205 } else if (jRow==iRow) { 207 longDouble value=element[j]*multiplier;206 CoinWorkDouble value=element[j]*multiplier; 208 207 diagonalValue += value; 209 208 } … … 212 211 } 213 212 for (int j=iRow+1;j<numberRows_;j++) 214 largest2 = CoinMax(largest2, fabs(work[j]));213 largest2 = CoinMax(largest2,CoinAbs(work[j])); 215 214 diagonal_[iRow]=diagonalValue; 216 largest2 = CoinMax(largest2, fabs(diagonalValue));215 largest2 = CoinMax(largest2,CoinAbs(diagonalValue)); 217 216 } else { 218 / / dropped217 /* dropped*/ 219 218 diagonal_[iRow]=1.0; 220 219 } … … 222 221 work += addOffset; 223 222 } 224 / /check sizes223 /*check sizes*/ 225 224 largest2*=1.0e-20; 226 largest = CoinMin(largest2, 1.0e-11);225 largest = CoinMin(largest2,CHOL_SMALL_VALUE); 227 226 int numberDroppedBefore=0; 228 227 for (iRow=0;iRow<numberRows_;iRow++) { 229 228 int dropped=rowsDropped_[iRow]; 230 / / Move to int array229 /* Move to int array*/ 231 230 rowsDropped[iRow]=dropped; 232 231 if (!dropped) { 233 longDouble diagonal = diagonal_[iRow];232 CoinWorkDouble diagonal = diagonal_[iRow]; 234 233 if (diagonal>largest2) { 235 234 diagonal_[iRow]=diagonal+perturbation; … … 245 244 doubleParameters_[3]=0.0; 246 245 doubleParameters_[4]=COIN_DBL_MAX; 247 integerParameters_[34]=0; / / say all must be positive246 integerParameters_[34]=0; /* say all must be positive*/ 248 247 #ifdef CHOL_COMPARE 249 248 if (numberRows_<200) … … 257 256 std::cout<<"Cholesky - largest "<<largest<<" smallest "<<smallest<<std::endl; 258 257 choleskyCondition_=largest/smallest; 259 / /drop fresh makes some formADAT easier258 /*drop fresh makes some formADAT easier*/ 260 259 if (newDropped||numberRowsDropped_) { 261 260 newDropped=0; … … 264 263 rowsDropped_[i]=dropped; 265 264 if (dropped==2) { 266 / /dropped this time265 /*dropped this time*/ 267 266 rowsDropped[newDropped++]=i; 268 267 rowsDropped_[i]=0; … … 273 272 } 274 273 } else { 275 / / KKT274 /* KKT*/ 276 275 CoinPackedMatrix * quadratic = NULL; 277 276 ClpQuadraticObjective * quadraticObj = … … 279 278 if (quadraticObj) 280 279 quadratic = quadraticObj->quadraticObjective(); 281 / / matrix280 /* matrix*/ 282 281 int numberRowsModel = model_->numberRows(); 283 282 int numberColumns = model_->numberColumns(); 284 283 int numberTotal = numberColumns + numberRowsModel; 285 284 longDouble * work = sparseFactor_; 286 work--; / / skip diagonal285 work--; /* skip diagonal*/ 287 286 int addOffset=numberRows_-1; 288 287 int iColumn; 289 288 if (!quadratic) { 290 289 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) { 293 292 value = 1.0/value; 294 largest = CoinMax(largest, fabs(value));293 largest = CoinMax(largest,CoinAbs(value)); 295 294 diagonal_[iColumn] = -value; 296 295 CoinBigIndex start=columnStart[iColumn]; 297 296 CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn]; 298 297 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];*/ 301 300 work[row[j]+numberTotal]=element[j]; 302 largest = CoinMax(largest, fabs(element[j]));301 largest = CoinMax(largest,CoinAbs(element[j])); 303 302 } 304 303 } else { … … 309 308 } 310 309 } else { 311 / / Quadratic310 /* Quadratic*/ 312 311 const int * columnQuadratic = quadratic->getIndices(); 313 312 const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts(); … … 315 314 const double * quadraticElement = quadratic->getElements(); 316 315 for (iColumn=0;iColumn<numberColumns;iColumn++) { 317 longDouble value = diagonal[iColumn];316 CoinWorkDouble value = diagonal[iColumn]; 318 317 CoinBigIndex j; 319 if ( fabs(value)>1.0e-100) {318 if (CoinAbs(value)>1.0e-100) { 320 319 value = 1.0/value; 321 320 for (j=columnQuadraticStart[iColumn]; … … 328 327 } 329 328 } 330 largest = CoinMax(largest, fabs(value));329 largest = CoinMax(largest,CoinAbs(value)); 331 330 diagonal_[iColumn] = -value; 332 331 CoinBigIndex start=columnStart[iColumn]; … … 334 333 for (j=start;j<end;j++) { 335 334 work[row[j]+numberTotal]=element[j]; 336 largest = CoinMax(largest, fabs(element[j]));335 largest = CoinMax(largest,CoinAbs(element[j])); 337 336 } 338 337 } else { … … 344 343 } 345 344 } 346 / / slacks345 /* slacks*/ 347 346 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) { 350 349 value = 1.0/value; 351 largest = CoinMax(largest, fabs(value));350 largest = CoinMax(largest,CoinAbs(value)); 352 351 } else { 353 352 value = 1.0e100; … … 358 357 work += addOffset; 359 358 } 360 / / Finish diagonal359 /* Finish diagonal*/ 361 360 for (iRow=0;iRow<numberRowsModel;iRow++) { 362 361 diagonal_[iRow+numberTotal]=delta2; 363 362 } 364 / /check sizes363 /*check sizes*/ 365 364 largest*=1.0e-20; 366 largest = CoinMin(largest, 1.0e-11);365 largest = CoinMin(largest,CHOL_SMALL_VALUE); 367 366 doubleParameters_[10]=CoinMax(1.0e-20,largest); 368 367 integerParameters_[20]=0; 369 368 doubleParameters_[3]=0.0; 370 369 doubleParameters_[4]=COIN_DBL_MAX; 371 / / Set up LDL cutoff370 /* Set up LDL cutoff*/ 372 371 integerParameters_[34]=numberTotal; 373 / / KKT372 /* KKT*/ 374 373 int * rowsDropped2 = new int[numberRows_]; 375 374 CoinZeroN(rowsDropped2,numberRows_); … … 385 384 std::cout<<"Cholesky - largest "<<largest<<" smallest "<<smallest<<std::endl; 386 385 choleskyCondition_=largest/smallest; 387 / / Should save adjustments in ..R_386 /* Should save adjustments in ..R_*/ 388 387 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(); 391 390 for (iRow=0;iRow<numberTotal;iRow++) { 392 391 if (rowsDropped2[iRow]) { 393 392 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;*/ 396 395 rowsDropped_[iRow]=0; 397 396 primalR[iRow]=doubleParameters_[20]; … … 404 403 if (rowsDropped2[iRow]) { 405 404 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;*/ 408 407 rowsDropped_[iRow]=0; 409 408 dualR[iRow-numberTotal]=doubleParameters_[34]; … … 425 424 diagonal_ = sparseFactor_ + 40000; 426 425 sparseFactor_=diagonal_ + numberRows_; 427 / /memcpy(sparseFactor_,xx,sizeFactor_*sizeof(double));426 /*memcpy(sparseFactor_,xx,sizeFactor_*sizeof(double));*/ 428 427 CoinMemcpyN(xx,40000,sparseFactor_); 429 428 CoinMemcpyN(yy,numberRows_,diagonal_); 430 429 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; 433 432 double dropValue = doubleParameters_[10]; 434 433 int firstPositive=integerParameters_[34]; 435 434 longDouble * work = sparseFactor_; 436 / / Allow for triangular435 /* Allow for triangular*/ 437 436 int addOffset=numberRows_-1; 438 437 work--; … … 441 440 int addOffsetNow = numberRows_-1;; 442 441 longDouble * workNow=sparseFactor_-1+iColumn; 443 double diagonalValue = diagonal_[iColumn];442 CoinWorkDouble diagonalValue = diagonal_[iColumn]; 444 443 for (iRow=0;iRow<iColumn;iRow++) { 445 444 double aj = *workNow; 446 445 addOffsetNow--; 447 446 workNow += addOffsetNow; 448 double multiplier = workDouble_[iRow]; 449 diagonalValue -= aj*aj*multiplier; 447 diagonalValue -= aj*aj*workDouble_[iRow]; 450 448 } 451 449 bool dropColumn=false; 452 450 if (iColumn<firstPositive) { 453 / / must be negative451 /* must be negative*/ 454 452 if (diagonalValue<=-dropValue) { 455 453 smallest = CoinMin(smallest,-diagonalValue); … … 464 462 } 465 463 } else { 466 / / must be positive464 /* must be positive*/ 467 465 if (diagonalValue>=dropValue) { 468 466 smallest = CoinMin(smallest,diagonalValue); … … 494 492 } 495 493 } else { 496 / / drop column494 /* drop column*/ 497 495 rowsDropped[iColumn]=2; 498 496 numberDropped++; … … 511 509 diagonal_=yy; 512 510 } 513 / /#define POS_DEBUG511 /*#define POS_DEBUG*/ 514 512 #ifdef POS_DEBUG 515 513 static int counter=0; … … 522 520 iCol=0; 523 521 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);*/ 525 523 k=kk; 526 524 while (k>=numberBlocks) { … … 541 539 { 542 540 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));*/ 550 548 int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT; 551 / / later align on boundary549 /* later align on boundary*/ 552 550 longDouble * a = sparseFactor_+BLOCKSQ*numberBlocks; 553 551 int n=numberRows_; 554 552 int nRound = numberRows_&(~(BLOCK-1)); 555 / / adjust if exact553 /* adjust if exact*/ 556 554 if (nRound==n) 557 555 nRound -= BLOCK; 558 556 int sizeLastBlock=n-nRound; 559 int get = n*(n-1)/2; / / ? as no diagonal557 int get = n*(n-1)/2; /* ? as no diagonal*/ 560 558 int block = numberBlocks*(numberBlocks+1)/2; 561 559 int ifOdd; … … 566 564 ifOdd=1; 567 565 int put = BLOCKSQ; 568 / / do last separately566 /* do last separately*/ 569 567 put -= (BLOCK-sizeLastBlock)*(BLOCK+1); 570 568 for (iColumn=numberRows_-1;iColumn>=nRound;iColumn--) { … … 575 573 assert (aa+put2>=sparseFactor_+get); 576 574 } 577 / / save diagonal as well575 /* save diagonal as well*/ 578 576 aa[--put2]=diagonal_[iColumn]; 579 577 } … … 581 579 block--; 582 580 } else { 583 / / exact fit581 /* exact fit*/ 584 582 rowLast=numberRows_-1; 585 583 ifOdd=0; 586 584 } 587 / / Now main loop585 /* Now main loop*/ 588 586 int nBlock=0; 589 587 for (;n>0;n-=BLOCK) { … … 592 590 int put = BLOCKSQ; 593 591 int putLast=0; 594 / / see if we have small block592 /* see if we have small block*/ 595 593 if (ifOdd) { 596 594 aaLast = &a[(block-1)*BLOCKSQ]; … … 600 598 for (iColumn=n-1;iColumn>=n-BLOCK;iColumn--) { 601 599 if (aaLast) { 602 / / last bit600 /* last bit*/ 603 601 for (int iRow=numberRows_-1;iRow>rowLast;iRow--) { 604 602 aaLast[--putLast] = sparseFactor_[--get]; … … 617 615 } 618 616 if (j-BLOCK<iColumn) { 619 / / save diagonal as well617 /* save diagonal as well*/ 620 618 aPut[--put2]=diagonal_[iColumn]; 621 619 } … … 628 626 block -= nBlock+ifOdd; 629 627 } 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, 631 634 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*/ 660 void 661 ClpCholeskyCfactor(ClpCholeskyDenseC * thisStruct, longDouble * a, int n, int numberBlocks, 638 662 longDouble * diagonal, longDouble * work, int * rowsDropped) 639 663 { 640 664 if (n<=BLOCK) { 641 factorLeaf(a,n,diagonal,work,rowsDropped);665 ClpCholeskyCfactorLeaf(thisStruct, a,n,diagonal,work,rowsDropped); 642 666 } else { 643 667 int nb=number_blocks((n+1)>>1); … … 647 671 int nintri=(nb*(nb+1))>>1; 648 672 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); 651 675 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, 654 678 numberBlocks-nb,diagonal+nThis,work+nThis,rowsDropped); 655 679 } 656 680 } 657 // Non leaf recursive triangle rectangle update 681 /* Non leaf recursive triangle rectangle update*/ 658 682 void 659 ClpCholesky Dense::triRec(longDouble * aTri, int nThis, longDouble * aUnder,683 ClpCholeskyCtriRec(ClpCholeskyDenseC * thisStruct, longDouble * aTri, int nThis, longDouble * aUnder, 660 684 longDouble * diagonal, longDouble * work, 661 685 int nLeft, int iBlock, int jBlock, … … 663 687 { 664 688 if (nThis<=BLOCK&&nLeft<=BLOCK) { 665 triRecLeaf(aTri,aUnder,diagonal,work,nLeft);689 ClpCholeskyCtriRecLeaf(/*thisStruct,*/ aTri,aUnder,diagonal,work,nLeft); 666 690 } else if (nThis<nLeft) { 667 691 int nb=number_blocks((nLeft+1)>>1); 668 692 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, 671 695 iBlock+nb,jBlock,numberBlocks); 672 696 } else { … … 678 702 int nintri=(nb*(nb+1))>>1; 679 703 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); 681 705 /* and rectangular update */ 682 706 i=((numberBlocks-jBlock)*(numberBlocks-jBlock-1)- 683 707 (numberBlocks-jBlock-nb)*(numberBlocks-jBlock-nb-1))>>1; 684 708 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, 686 710 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, 688 712 work+nThis2,nLeft, 689 713 iBlock-nb,kBlock-nb,numberBlocks-nb); 690 714 } 691 715 } 692 // Non leaf recursive rectangle triangle update 716 /* Non leaf recursive rectangle triangle update*/ 693 717 void 694 ClpCholesky Dense::recTri(longDouble * aUnder, int nTri, int nDo,718 ClpCholeskyCrecTri(ClpCholeskyDenseC * thisStruct, longDouble * aUnder, int nTri, int nDo, 695 719 int iBlock, int jBlock,longDouble * aTri, 696 720 longDouble * diagonal, longDouble * work, … … 698 722 { 699 723 if (nTri<=BLOCK&&nDo<=BLOCK) { 700 recTriLeaf(aUnder,aTri,diagonal,work,nTri);724 ClpCholeskyCrecTriLeaf(/*thisStruct,*/ aUnder,aTri,/*diagonal,*/work,nTri); 701 725 } else if (nTri<nDo) { 702 726 int nb=number_blocks((nDo+1)>>1); … … 704 728 longDouble * aother; 705 729 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); 707 731 i=((numberBlocks-jBlock)*(numberBlocks-jBlock-1)- 708 732 (numberBlocks-jBlock-nb)*(numberBlocks-jBlock-nb-1))>>1; 709 733 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, 711 735 work+nDo2,numberBlocks-nb); 712 736 } else { … … 715 739 longDouble * aother; 716 740 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); 718 742 /* and rectangular update */ 719 743 i=((numberBlocks-iBlock)*(numberBlocks-iBlock+1)- 720 744 (numberBlocks-iBlock-nb)*(numberBlocks-iBlock-nb+1))>>1; 721 745 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, 723 747 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, 725 749 aTri+number_entries(i),diagonal,work,numberBlocks); 726 750 } … … 731 755 */ 732 756 void 733 ClpCholesky Dense::recRec(longDouble * above, int nUnder, int nUnderK,757 ClpCholeskyCrecRec(ClpCholeskyDenseC * thisStruct, longDouble * above, int nUnder, int nUnderK, 734 758 int nDo, longDouble * aUnder, longDouble *aOther, 735 759 longDouble * work, … … 739 763 if (nDo<=BLOCK&&nUnder<=BLOCK&&nUnderK<=BLOCK) { 740 764 assert (nDo==BLOCK&&nUnder==BLOCK); 741 recRecLeaf(above , aUnder , aOther, work, nUnderK);765 ClpCholeskyCrecRecLeaf(/*thisStruct,*/ above , aUnder , aOther, work, nUnderK); 742 766 } else if (nDo<=nUnderK&&nUnder<=nUnderK) { 743 767 int nb=number_blocks((nUnderK+1)>>1); 744 768 int nUnder2=number_rows(nb); 745 recRec(above,nUnder,nUnder2,nDo,aUnder,aOther,work,769 ClpCholeskyCrecRec(thisStruct, above,nUnder,nUnder2,nDo,aUnder,aOther,work, 746 770 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), 748 772 aOther+number_entries(nb),work,iBlock,jBlock,numberBlocks); 749 773 } else if (nUnderK<=nDo&&nUnder<=nDo) { … … 751 775 int nDo2=number_rows(nb); 752 776 int i; 753 recRec(above,nUnder,nUnderK,nDo2,aUnder,aOther,work,777 ClpCholeskyCrecRec(thisStruct, above,nUnder,nUnderK,nDo2,aUnder,aOther,work, 754 778 iBlock,jBlock,numberBlocks); 755 779 i=((numberBlocks-jBlock)*(numberBlocks-jBlock-1)- 756 780 (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, 758 782 aUnder+number_entries(i), 759 783 aOther,work+nDo2, … … 763 787 int nUnder2=number_rows(nb); 764 788 int i; 765 recRec(above,nUnder2,nUnderK,nDo,aUnder,aOther,work,789 ClpCholeskyCrecRec(thisStruct, above,nUnder2,nUnderK,nDo,aUnder,aOther,work, 766 790 iBlock,jBlock,numberBlocks); 767 791 i=((numberBlocks-iBlock)*(numberBlocks-iBlock-1)- 768 792 (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, 770 794 aOther+number_entries(i),work,iBlock+nb,jBlock,numberBlocks); 771 795 } 772 796 } 773 // Leaf recursive factor 797 /* Leaf recursive factor*/ 774 798 void 775 ClpCholesky Dense::factorLeaf(longDouble * a, int n,799 ClpCholeskyCfactorLeaf(ClpCholeskyDenseC * thisStruct, longDouble * a, int n, 776 800 longDouble * diagonal, longDouble * work, int * rowsDropped) 777 801 { 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_; 784 805 int i, j, k; 785 longDouble t00,temp1;806 CoinWorkDouble t00,temp1; 786 807 longDouble * aa; 787 808 aa = a-BLOCK; 788 809 for (j = 0; j < n; j ++) { 810 bool dropColumn; 811 CoinWorkDouble useT00; 789 812 aa+=BLOCK; 790 813 t00 = aa[j]; 791 814 for (k = 0; k < j; ++k) { 792 longDouble multiplier = work[k];815 CoinWorkDouble multiplier = work[k]; 793 816 t00 -= a[j + k * BLOCK] * a[j + k * BLOCK] * multiplier; 794 817 } 795 booldropColumn=false;796 longDoubleuseT00=t00;797 if ( (int)(j+rowOffset)<firstPositive) {798 / / must be negative818 dropColumn=false; 819 useT00=t00; 820 if (j+rowOffset<firstPositive) { 821 /* must be negative*/ 799 822 if (t00<=-dropValue) { 800 smallest = CoinMin(smallest,-t00); 801 largest = CoinMax(largest,-t00); 802 //aa[j]=t00; 823 /*aa[j]=t00;*/ 803 824 t00 = 1.0/t00; 804 825 } else { 805 826 dropColumn=true; 806 / /aa[j]=-1.0e100;827 /*aa[j]=-1.0e100;*/ 807 828 useT00=-1.0e-100; 808 829 t00=0.0; 809 integerParameters_[20]++;810 830 } 811 831 } else { 812 / / must be positive832 /* must be positive*/ 813 833 if (t00>=dropValue) { 814 smallest = CoinMin(smallest,t00); 815 largest = CoinMax(largest,t00); 816 //aa[j]=t00; 834 /*aa[j]=t00;*/ 817 835 t00 = 1.0/t00; 818 836 } else { 819 837 dropColumn=true; 820 / /aa[j]=1.0e100;838 /*aa[j]=1.0e100;*/ 821 839 useT00=1.0e-100; 822 840 t00=0.0; 823 integerParameters_[20]++;824 841 } 825 842 } … … 831 848 t00=aa[i]; 832 849 for (k = 0; k < j; ++k) { 833 longDouble multiplier = work[k];850 CoinWorkDouble multiplier = work[k]; 834 851 t00 -= a[i + k * BLOCK] * a[j + k * BLOCK] * multiplier; 835 852 } … … 837 854 } 838 855 } else { 839 / / drop column856 /* drop column*/ 840 857 rowsDropped[j+rowOffset]=2; 841 numberDropped++;842 858 diagonal[j]=0.0; 843 / /aa[j]=1.0e100;859 /*aa[j]=1.0e100;*/ 844 860 work[j]=1.0e100; 845 861 for (i=j+1;i<n;i++) { … … 848 864 } 849 865 } 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*/ 855 868 void 856 ClpCholesky Dense::triRecLeaf(longDouble * aTri, longDouble * aUnder, longDouble * diagonal, longDouble * work,869 ClpCholeskyCtriRecLeaf(/*ClpCholeskyDenseC * thisStruct,*/ longDouble * aTri, longDouble * aUnder, longDouble * diagonal, longDouble * work, 857 870 int nUnder) 858 871 { … … 862 875 int irt,ict; 863 876 int it=bNumber(aTri,irt,ict); 864 / /printf("%d %d\n",iu,it);877 /*printf("%d %d\n",iu,it);*/ 865 878 printf("trirecleaf under (%d,%d), tri (%d,%d)\n", 866 879 iru,icu,irt,ict); 867 assert (diagonal== diagonal_+ict*BLOCK);880 assert (diagonal==thisStruct->diagonal_+ict*BLOCK); 868 881 #endif 869 882 int j; … … 873 886 aa = aTri-2*BLOCK; 874 887 for (j = 0; j < BLOCK; j +=2) { 888 int i; 889 CoinWorkDouble temp0 = diagonal[j]; 890 CoinWorkDouble temp1 = diagonal[j+1]; 875 891 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 ( intk = 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]; 889 905 t00 -= au0 * at0; 890 906 t10 -= au0 * at1; … … 893 909 } 894 910 t00 *= temp0; 895 longDoubleat1 = aTri[j + 1 + j*BLOCK]*work[j];911 at1 = aTri[j + 1 + j*BLOCK]*work[j]; 896 912 t10 -= t00 * at1; 897 913 t01 *= temp0; … … 907 923 aa = aTri-BLOCK; 908 924 for (j = 0; j < BLOCK; j ++) { 925 int i; 926 CoinWorkDouble temp1 = diagonal[j]; 909 927 aa+=BLOCK; 910 longDouble temp1 = diagonal[j];911 for (int i=0;i<nUnder;i++) { 912 longDouble t00=aUnder[i+j*BLOCK];913 for ( intk = 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]; 915 933 t00 -= aUnder[i + k * BLOCK] * aTri[j + k * BLOCK] * multiplier; 916 934 } … … 922 940 #endif 923 941 } 924 // Leaf recursive rectangle triangle update 925 void ClpCholesky Dense::recTriLeaf(longDouble * aUnder, longDouble * aTri,926 longDouble * diagonal,longDouble * work, int nUnder)942 /* Leaf recursive rectangle triangle update*/ 943 void ClpCholeskyCrecTriLeaf(/*ClpCholeskyDenseC * thisStruct,*/ longDouble * aUnder, longDouble * aTri, 944 /*longDouble * diagonal,*/ longDouble * work, int nUnder) 927 945 { 928 946 #ifdef POS_DEBUG … … 931 949 int irt,ict; 932 950 int it=bNumber(aTri,irt,ict); 933 / /printf("%d %d\n",iu,it);951 /*printf("%d %d\n",iu,it);*/ 934 952 printf("rectrileaf under (%d,%d), tri (%d,%d)\n", 935 953 iru,icu,irt,ict); 936 assert (diagonal== diagonal_+icu*BLOCK);954 assert (diagonal==thisStruct->diagonal_+icu*BLOCK); 937 955 #endif 938 956 int i, j, k; 939 longDouble t00;957 CoinWorkDouble t00; 940 958 longDouble * aa; 941 959 #ifdef BLOCKUNROLL … … 944 962 aa = aTri-2*BLOCK; 945 963 for (j = 0; j < BLOCK; j +=2) { 946 longDouble t00,t01,t10,t11;964 CoinWorkDouble t00,t01,t10,t11; 947 965 aa+=2*BLOCK; 948 966 aUnder2+=2; … … 951 969 t10=aa[j+1+BLOCK]; 952 970 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; 958 976 t00 -= a0 * x0; 959 977 t01 -= a1 * x0; … … 969 987 t11=aa[i+1+BLOCK]; 970 988 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; 974 992 t00 -= aUnder[i + k * BLOCK] * a0; 975 993 t01 -= aUnder[i + k * BLOCK] * a1; … … 991 1009 t00=aa[i]; 992 1010 for (k = 0; k < BLOCK; ++k) { 993 longDouble multiplier = work[k];1011 CoinWorkDouble multiplier = work[k]; 994 1012 t00 -= aUnder[i + k * BLOCK] * aUnder[j + k * BLOCK]*multiplier; 995 1013 } … … 1006 1024 */ 1007 1025 void 1008 ClpCholeskyDense::recRecLeaf(const longDouble * COIN_RESTRICT above, 1026 ClpCholeskyCrecRecLeaf(/*ClpCholeskyDenseC * thisStruct,*/ 1027 const longDouble * COIN_RESTRICT above, 1009 1028 const longDouble * COIN_RESTRICT aUnder, 1010 1029 longDouble * COIN_RESTRICT aOther, … … 1019 1038 int iro,ico; 1020 1039 int io=bNumber(aOther,iro,ico); 1021 / /printf("%d %d %d\n",ia,iu,io);1040 /*printf("%d %d %d\n",ia,iu,io);*/ 1022 1041 printf("recrecleaf above (%d,%d), under (%d,%d), other (%d,%d)\n", 1023 1042 ira,ica,iru,icu,iro,ico); … … 1028 1047 aa = aOther-4*BLOCK; 1029 1048 if (nUnder==BLOCK) { 1030 / /#define INTEL1049 /*#define INTEL*/ 1031 1050 #ifdef INTEL 1032 1051 aa+=2*BLOCK; … … 1034 1053 aa+=2*BLOCK; 1035 1054 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]; 1040 1059 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; 1044 1063 t00 -= a00 * above[j + 0 + k * BLOCK]; 1045 1064 t10 -= a00 * above[j + 1 + k * BLOCK]; … … 1057 1076 aa+=4*BLOCK; 1058 1077 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]; 1075 1094 const longDouble * COIN_RESTRICT aUnderNow = aUnder+i; 1076 1095 const longDouble * COIN_RESTRICT aboveNow = above+j; 1077 1096 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; 1083 1102 t00 -= a00 * aboveNow[0]; 1084 1103 t10 -= a00 * aboveNow[1]; … … 1125 1144 aa+=4*BLOCK; 1126 1145 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]; 1135 1154 const longDouble * COIN_RESTRICT aUnderNow = aUnder+i; 1136 1155 const longDouble * COIN_RESTRICT aboveNow = above+j; 1137 1156 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; 1141 1160 t00 -= a00 * aboveNow[0]; 1142 1161 t10 -= a00 * aboveNow[1]; … … 1160 1179 } 1161 1180 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; 1167 1186 for (k=0;k<BLOCK;k++) { 1168 1187 a0=aUnder[n+k*BLOCK]*work[k]; … … 1184 1203 aa+=BLOCK; 1185 1204 for (i=0;i< nUnder;i++) { 1186 longDouble t00=aa[i+0*BLOCK];1205 CoinWorkDouble t00=aa[i+0*BLOCK]; 1187 1206 for (k=0;k<BLOCK;k++) { 1188 longDouble a00=aUnder[i+k*BLOCK]*work[k];1207 CoinWorkDouble a00=aUnder[i+k*BLOCK]*work[k]; 1189 1208 t00 -= a00 * above[j + k * BLOCK]; 1190 1209 } … … 1196 1215 /* Uses factorization to solve. */ 1197 1216 void 1198 ClpCholeskyDense::solve ( double * region)1217 ClpCholeskyDense::solve (CoinWorkDouble * region) 1199 1218 { 1200 1219 #ifdef CHOL_COMPARE … … 1228 1247 } 1229 1248 #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_;*/ 1234 1253 int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT; 1235 / / later align on boundary1254 /* later align on boundary*/ 1236 1255 longDouble * a = sparseFactor_+BLOCKSQ*numberBlocks; 1237 1256 int iBlock; … … 1261 1280 aa+=BLOCKSQ; 1262 1281 } 1263 / / do diagonal outside1282 /* do diagonal outside*/ 1264 1283 for (iColumn=0;iColumn<numberRows_;iColumn++) 1265 1284 region[iColumn] *= diagonal_[iColumn]; … … 1293 1312 if (numberRows_<200) { 1294 1313 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); 1296 1315 } 1297 1316 delete [] region2; … … 1299 1318 #endif 1300 1319 } 1301 / / Forward part of solve 11320 /* Forward part of solve 1*/ 1302 1321 void 1303 ClpCholeskyDense::solveF1(longDouble * a,int n, double * region)1322 ClpCholeskyDense::solveF1(longDouble * a,int n,CoinWorkDouble * region) 1304 1323 { 1305 1324 int j, k; 1306 longDouble t00;1325 CoinWorkDouble t00; 1307 1326 for (j = 0; j < n; j ++) { 1308 1327 t00 = region[j]; … … 1310 1329 t00 -= region[k] * a[j + k * BLOCK]; 1311 1330 } 1312 / /t00*=a[j + j * BLOCK];1331 /*t00*=a[j + j * BLOCK];*/ 1313 1332 region[j] = t00; 1314 1333 } 1315 1334 } 1316 / / Forward part of solve 21335 /* Forward part of solve 2*/ 1317 1336 void 1318 ClpCholeskyDense::solveF2(longDouble * a,int n, double * region, double * region2)1337 ClpCholeskyDense::solveF2(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2) 1319 1338 { 1320 1339 int j, k; … … 1322 1341 if (n==BLOCK) { 1323 1342 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]; 1328 1347 t0 -= region[0] * a[0 + 0 * BLOCK]; 1329 1348 t1 -= region[0] * a[1 + 0 * BLOCK]; … … 1416 1435 #endif 1417 1436 for (k = 0; k < n; ++k) { 1418 longDouble t00 = region2[k];1437 CoinWorkDouble t00 = region2[k]; 1419 1438 for (j = 0; j < BLOCK; j ++) { 1420 1439 t00 -= region[j] * a[k + j * BLOCK]; … … 1426 1445 #endif 1427 1446 } 1428 / / Backward part of solve 11447 /* Backward part of solve 1*/ 1429 1448 void 1430 ClpCholeskyDense::solveB1(longDouble * a,int n, double * region)1449 ClpCholeskyDense::solveB1(longDouble * a,int n,CoinWorkDouble * region) 1431 1450 { 1432 1451 int j, k; 1433 longDouble t00;1452 CoinWorkDouble t00; 1434 1453 for (j = n-1; j >=0; j --) { 1435 1454 t00 = region[j]; … … 1437 1456 t00 -= region[k] * a[k + j * BLOCK]; 1438 1457 } 1439 / /t00*=a[j + j * BLOCK];1458 /*t00*=a[j + j * BLOCK];*/ 1440 1459 region[j] = t00; 1441 1460 } 1442 1461 } 1443 / / Backward part of solve 21462 /* Backward part of solve 2*/ 1444 1463 void 1445 ClpCholeskyDense::solveB2(longDouble * a,int n, double * region, double * region2)1464 ClpCholeskyDense::solveB2(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2) 1446 1465 { 1447 1466 int j, k; … … 1449 1468 if (n==BLOCK) { 1450 1469 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]; 1455 1474 t0 -= region2[0] * a[0 + 0*BLOCK]; 1456 1475 t1 -= region2[0] * a[0 + 1*BLOCK]; … … 1544 1563 #endif 1545 1564 for (j = 0; j < BLOCK; j ++) { 1546 longDouble t00 = region[j];1565 CoinWorkDouble t00 = region[j]; 1547 1566 for (k = 0; k < n; ++k) { 1548 1567 t00 -= region2[k] * a[k + j * BLOCK]; … … 1554 1573 #endif 1555 1574 } 1556 /* Uses factorization to solve. */1557 void1558 ClpCholeskyDense::solveLong (longDouble * region)1559 {1560 int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;1561 // later align on boundary1562 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 outside1590 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 11620 void1621 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 21635 void1636 ClpCholeskyDense::solveF2Long(longDouble * a,int n,longDouble * region, longDouble * region2)1637 {1638 int j, k;1639 #ifdef BLOCKUNROLL1640 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>81686 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 #endif1726 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 #endif1735 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 BLOCKUNROLL1743 }1744 #endif1745 }1746 // Backward part of solve 11747 void1748 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];