Changeset 14 for branches/devel1/ClpSimplexPrimal.cpp
 Timestamp:
 Aug 7, 2002 11:55:35 AM (18 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

branches/devel1/ClpSimplexPrimal.cpp
r13 r14 216 216 assert(rowUpper_); 217 217 218 #ifdef CLP_DEBUG219 int debugIteration=1;220 #endif221 218 222 219 algorithm_ = +1; … … 249 246 perturb(); 250 247 251 double objectiveChange;252 248 // for primal we will change bounds using infeasibilityCost_ 253 249 if (nonLinearCost_==NULL) { … … 343 339 saveNumber = numberIterations_; 344 340 345 // status stays at 1 while iterating, >=0 finished, 2 to invert 346 // status 3 to go to top without an invert 347 while (problemStatus_==1) { 341 // Iterate 342 whileIterating(firstSuperBasic); 343 } 344 345 // if infeasible get real values 346 if (problemStatus_) { 347 infeasibilityCost_=0.0; 348 createRim(7); 349 nonLinearCost_>checkInfeasibilities(true); 350 sumPrimalInfeasibilities_=nonLinearCost_>sumInfeasibilities(); 351 numberPrimalInfeasibilities_= nonLinearCost_>numberInfeasibilities(); 352 // and get good feasible duals 353 computeDuals(); 354 } 355 // at present we are leaving factorization around 356 // maybe we should empty it 357 deleteRim(); 358 handler_>message(CLP_SIMPLEX_FINISHED+problemStatus_,messages_) 359 <<objectiveValue() 360 <<OsiMessageEol; 361 // Restore any saved stuff 362 perturbation_ = savePerturbation; 363 factorization_>sparseThreshold(saveSparse); 364 infeasibilityCost_ = saveInfeasibilityCost; 365 return problemStatus_; 366 } 367 void 368 ClpSimplexPrimal::whileIterating(int firstSuperBasic) 369 { 370 371 // Say if values pass 372 int ifValuesPass=0; 373 if (firstSuperBasic<numberRows_+numberColumns_) 374 ifValuesPass=1; 375 int saveNumber = numberIterations_; 376 // status stays at 1 while iterating, >=0 finished, 2 to invert 377 // status 3 to go to top without an invert 378 while (problemStatus_==1) { 348 379 #ifdef CLP_DEBUG 349 350 351 352 353 354 355 356 357 358 359 380 { 381 int i; 382 // not [1] as has information 383 for (i=0;i<4;i++) { 384 if (i!=1) 385 rowArray_[i]>checkClear(); 386 } 387 for (i=0;i<2;i++) { 388 columnArray_[i]>checkClear(); 389 } 390 } 360 391 #endif 361 392 #if CLP_DEBUG>2 362 // very expensive 363 if (numberIterations_>0&&numberIterations_<2534) { 364 handler_>setLogLevel(63); 365 double saveValue = objectiveValue_; 366 double * saveRow1 = new double[numberRows_]; 367 double * saveRow2 = new double[numberRows_]; 368 memcpy(saveRow1,rowReducedCost_,numberRows_*sizeof(double)); 369 memcpy(saveRow2,rowActivityWork_,numberRows_*sizeof(double)); 370 double * saveColumn1 = new double[numberColumns_]; 371 double * saveColumn2 = new double[numberColumns_]; 372 memcpy(saveColumn1,reducedCostWork_,numberColumns_*sizeof(double)); 373 memcpy(saveColumn2,columnActivityWork_,numberColumns_*sizeof(double)); 374 createRim(7); 375 gutsOfSolution(rowActivityWork_,columnActivityWork_); 376 printf("xxx %d old obj %g, recomputed %g, sum primal inf %g\n", 377 numberIterations_, 378 saveValue,objectiveValue_,sumPrimalInfeasibilities_); 379 memcpy(rowReducedCost_,saveRow1,numberRows_*sizeof(double)); 380 memcpy(rowActivityWork_,saveRow2,numberRows_*sizeof(double)); 381 memcpy(reducedCostWork_,saveColumn1,numberColumns_*sizeof(double)); 382 memcpy(columnActivityWork_,saveColumn2,numberColumns_*sizeof(double)); 383 delete [] saveRow1; 384 delete [] saveRow2; 385 delete [] saveColumn1; 386 delete [] saveColumn2; 387 objectiveValue_=saveValue; 388 } 389 #endif 390 #ifdef CLP_DEBUG 391 if(numberIterations_==debugIteration) { 392 printf("dodgy iteration coming up\n"); 393 } 394 #endif 395 if (!ifValuesPass) { 396 // choose column to come in 397 // can use pivotRow_ to update weights 398 // pass in list of cost changes so can do row updates (rowArray_[1]) 399 // NOTE rowArray_[0] is used by computeDuals which is a 400 // slow way of getting duals but might be used 401 primalColumn(rowArray_[1],rowArray_[2],rowArray_[3], 402 columnArray_[0],columnArray_[1]); 403 } else { 404 // in values pass 405 if (ifValuesPass>0) { 406 nextSuperBasic(firstSuperBasic); 407 if (firstSuperBasic==numberRows_+numberColumns_) 408 ifValuesPass=1; // signal end of values pass after this 409 } else { 410 // end of values pass  initialize weights etc 411 primalColumnPivot_>saveWeights(this,5); 412 ifValuesPass=0; 413 if(saveNumber != numberIterations_) { 414 problemStatus_=2; // factorize now 415 pivotRow_=1; // say no weights update 416 break; 417 } 418 419 // and get variable 420 primalColumn(rowArray_[1],rowArray_[2],rowArray_[3], 421 columnArray_[0],columnArray_[1]); 422 } 423 } 424 pivotRow_=1; 425 sequenceOut_=1; 426 rowArray_[1]>clear(); 427 if (sequenceIn_>=0) { 428 // we found a pivot column 429 #ifdef CLP_DEBUG 430 if ((handler_>logLevel()&32)) { 431 char x = isColumn(sequenceIn_) ? 'C' :'R'; 432 std::cout<<"pivot column "<< 433 x<<sequenceWithin(sequenceIn_)<<std::endl; 434 } 435 #endif 436 // update the incoming column 437 unpack(rowArray_[1]); 438 // save reduced cost 439 double saveDj = dualIn_; 440 factorization_>updateColumn(rowArray_[2],rowArray_[1],true); 441 // do ratio test and recompute dj 442 primalRow(rowArray_[1],rowArray_[3],rowArray_[2],rowArray_[0], 443 ifValuesPass); 444 if (ifValuesPass) { 445 saveDj=dualIn_; 446 if (pivotRow_==1(pivotRow_>=0&&fabs(alpha_)<1.0e5)) { 447 if(fabs(dualIn_)<1.0e2*dualTolerance_) { 448 // try other way 449 directionIn_=directionIn_; 450 primalRow(rowArray_[1],rowArray_[3],rowArray_[2],rowArray_[0], 451 0); 452 } 453 if (pivotRow_==1(pivotRow_>=0&&fabs(alpha_)<1.0e5)) { 454 // reject it 455 char x = isColumn(sequenceIn_) ? 'C' :'R'; 456 handler_>message(CLP_SIMPLEX_FLAG,messages_) 457 <<x<<sequenceWithin(sequenceIn_) 458 <<OsiMessageEol; 459 setFlagged(sequenceIn_); 460 lastBadIteration_ = numberIterations_; // say be more cautious 461 rowArray_[1]>clear(); 462 pivotRow_=1; 463 continue; 464 } 465 } 466 } 467 468 #ifdef CLP_DEBUG 469 if ((handler_>logLevel()&32)) 470 printf("btran dj %g, ftran dj %g\n",saveDj,dualIn_); 471 #endif 472 if (saveDj*dualIn_<1.0e20 473 fabs(saveDjdualIn_)>1.0e5*(1.0+fabs(saveDj))) { 474 handler_>message(CLP_PRIMAL_DJ,messages_) 475 <<saveDj<<dualIn_ 476 <<OsiMessageEol; 477 if(saveNumber != numberIterations_) { 478 problemStatus_=2; // factorize now 479 rowArray_[1]>clear(); 480 pivotRow_=1; // say no weights update 481 break; 482 } else { 483 // take on more relaxed criterion 484 if (saveDj*dualIn_<1.0e20 485 fabs(saveDjdualIn_)>1.0e4*(1.0+fabs(dualIn_))) { 486 // need to reject something 487 char x = isColumn(sequenceIn_) ? 'C' :'R'; 488 handler_>message(CLP_SIMPLEX_FLAG,messages_) 489 <<x<<sequenceWithin(sequenceIn_) 490 <<OsiMessageEol; 491 setFlagged(sequenceIn_); 492 lastBadIteration_ = numberIterations_; // say be more cautious 493 rowArray_[1]>clear(); 494 pivotRow_=1; 495 continue; 496 } 497 } 498 } 499 if (pivotRow_>=0) { 500 // if stable replace in basis 501 int updateStatus = factorization_>replaceColumn(rowArray_[2], 502 pivotRow_, 503 alpha_); 504 if (updateStatus==1) { 505 // slight error 506 if (factorization_>pivots()>5) 507 problemStatus_=2; // factorize now 508 } else if (updateStatus==2) { 509 // major error 510 // later we may need to unwind more e.g. fake bounds 511 if(saveNumber != numberIterations_) { 512 problemStatus_=2; // factorize now 513 break; 514 } else { 515 // need to reject something 516 char x = isColumn(sequenceIn_) ? 'C' :'R'; 517 handler_>message(CLP_SIMPLEX_FLAG,messages_) 518 <<x<<sequenceWithin(sequenceIn_) 519 <<OsiMessageEol; 520 setFlagged(sequenceIn_); 521 lastBadIteration_ = numberIterations_; // say be more cautious 522 rowArray_[1]>clear(); 523 pivotRow_=1; 524 continue; 525 } 526 } else if (updateStatus==3) { 527 // out of memory 528 // increase space if not many iterations 529 if (factorization_>pivots()< 530 0.5*factorization_>maximumPivots()&& 531 factorization_>pivots()<200) 532 factorization_>areaFactor( 533 factorization_>areaFactor() * 1.1); 534 problemStatus_=2; // factorize now 535 } 536 // here do part of steepest  ready for next iteration 537 primalColumnPivot_>updateWeights(rowArray_[1]); 538 } else { 539 if (pivotRow_==1) { 540 // no outgoing row is valid 541 #ifdef CLP_DEBUG 542 if (handler_>logLevel()&32) 543 printf("** no row pivot\n"); 544 #endif 545 if (!factorization_>pivots()) { 546 problemStatus_=5; //say looks unbounded 547 // do ray 548 delete [] ray_; 549 ray_ = new double [numberColumns_]; 550 CoinFillN(ray_,numberColumns_,0.0); 551 int number=rowArray_[1]>getNumElements(); 552 int * index = rowArray_[1]>getIndices(); 553 double * array = rowArray_[1]>denseVector(); 554 double way=directionIn_; 555 int i; 556 double zeroTolerance=1.0e12; 557 if (sequenceIn_<numberColumns_) 558 ray_[sequenceIn_]=directionIn_; 559 for (i=0;i<number;i++) { 560 int iRow=index[i]; 561 int iPivot=pivotVariable_[iRow]; 562 double arrayValue = array[iRow]; 563 if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance) 564 ray_[iPivot] = way* array[iRow]; 565 } 566 } 567 rowArray_[0]>clear(); 568 break; 569 } else { 570 // flipping from bound to bound 571 } 572 } 573 574 // update primal solution 575 576 objectiveChange=0.0; 577 // Cost on pivot row may change  may need to change dualIn 578 double oldCost=0.0; 579 if (pivotRow_>=0) 580 oldCost = cost(pivotVariable_[pivotRow_]); 581 // rowArray_[1] is not empty  used to update djs 582 updatePrimalsInPrimal(rowArray_[1],theta_, objectiveChange); 583 if (pivotRow_>=0) 584 dualIn_ += (oldCostcost(pivotVariable_[pivotRow_])); 585 586 int whatNext=housekeeping(objectiveChange); 587 588 if (whatNext==1) { 589 problemStatus_ =2; // refactorize 590 } else if (whatNext==2) { 591 // maximum iterations or equivalent 592 problemStatus_= 3; 593 break; 594 } 595 } else { 596 // no pivot column 597 #ifdef CLP_DEBUG 598 if (handler_>logLevel()&32) 599 printf("** no column pivot\n"); 600 #endif 601 if (nonLinearCost_>numberInfeasibilities()) 602 problemStatus_=4; // might be infeasible 603 break; 604 } 605 } 606 #if CLP_DEBUG>2 607 if (numberIterations_>620&&numberIterations_<2534) { 393 // very expensive 394 if (numberIterations_>0&&numberIterations_<2534) { 608 395 handler_>setLogLevel(63); 609 396 double saveValue = objectiveValue_; … … 632 419 } 633 420 #endif 634 } 635 636 // if infeasible get real values 637 if (problemStatus_) { 638 infeasibilityCost_=0.0; 639 createRim(7); 640 nonLinearCost_>checkInfeasibilities(true); 641 sumPrimalInfeasibilities_=nonLinearCost_>sumInfeasibilities(); 642 numberPrimalInfeasibilities_= nonLinearCost_>numberInfeasibilities(); 643 } 644 // at present we are leaving factorization around 645 // maybe we should empty it 646 deleteRim(); 647 handler_>message(CLP_SIMPLEX_FINISHED+problemStatus_,messages_) 648 <<objectiveValue() 649 <<OsiMessageEol; 650 // Restore any saved stuff 651 perturbation_ = savePerturbation; 652 factorization_>sparseThreshold(saveSparse); 653 infeasibilityCost_ = saveInfeasibilityCost; 654 return problemStatus_; 421 if (!ifValuesPass) { 422 // choose column to come in 423 // can use pivotRow_ to update weights 424 // pass in list of cost changes so can do row updates (rowArray_[1]) 425 // NOTE rowArray_[0] is used by computeDuals which is a 426 // slow way of getting duals but might be used 427 primalColumn(rowArray_[1],rowArray_[2],rowArray_[3], 428 columnArray_[0],columnArray_[1]); 429 } else { 430 // in values pass 431 if (ifValuesPass>0) { 432 nextSuperBasic(firstSuperBasic); 433 if (firstSuperBasic==numberRows_+numberColumns_) 434 ifValuesPass=1; // signal end of values pass after this 435 } else { 436 // end of values pass  initialize weights etc 437 primalColumnPivot_>saveWeights(this,5); 438 ifValuesPass=0; 439 if(saveNumber != numberIterations_) { 440 problemStatus_=2; // factorize now 441 pivotRow_=1; // say no weights update 442 break; 443 } 444 445 // and get variable 446 primalColumn(rowArray_[1],rowArray_[2],rowArray_[3], 447 columnArray_[0],columnArray_[1]); 448 } 449 } 450 pivotRow_=1; 451 sequenceOut_=1; 452 rowArray_[1]>clear(); 453 if (sequenceIn_>=0) { 454 // we found a pivot column 455 #ifdef CLP_DEBUG 456 if ((handler_>logLevel()&32)) { 457 char x = isColumn(sequenceIn_) ? 'C' :'R'; 458 std::cout<<"pivot column "<< 459 x<<sequenceWithin(sequenceIn_)<<std::endl; 460 } 461 #endif 462 // update the incoming column 463 unpack(rowArray_[1]); 464 // save reduced cost 465 double saveDj = dualIn_; 466 factorization_>updateColumn(rowArray_[2],rowArray_[1],true); 467 // do ratio test and recompute dj 468 primalRow(rowArray_[1],rowArray_[3],rowArray_[2],rowArray_[0], 469 ifValuesPass); 470 if (ifValuesPass) { 471 saveDj=dualIn_; 472 if (pivotRow_==1(pivotRow_>=0&&fabs(alpha_)<1.0e5)) { 473 if(fabs(dualIn_)<1.0e2*dualTolerance_) { 474 // try other way 475 directionIn_=directionIn_; 476 primalRow(rowArray_[1],rowArray_[3],rowArray_[2],rowArray_[0], 477 0); 478 } 479 if (pivotRow_==1(pivotRow_>=0&&fabs(alpha_)<1.0e5)) { 480 // reject it 481 char x = isColumn(sequenceIn_) ? 'C' :'R'; 482 handler_>message(CLP_SIMPLEX_FLAG,messages_) 483 <<x<<sequenceWithin(sequenceIn_) 484 <<OsiMessageEol; 485 setFlagged(sequenceIn_); 486 lastBadIteration_ = numberIterations_; // say be more cautious 487 rowArray_[1]>clear(); 488 pivotRow_=1; 489 continue; 490 } 491 } 492 } 493 494 #ifdef CLP_DEBUG 495 if ((handler_>logLevel()&32)) 496 printf("btran dj %g, ftran dj %g\n",saveDj,dualIn_); 497 #endif 498 if (saveDj*dualIn_<1.0e20 499 fabs(saveDjdualIn_)>1.0e5*(1.0+fabs(saveDj))) { 500 handler_>message(CLP_PRIMAL_DJ,messages_) 501 <<saveDj<<dualIn_ 502 <<OsiMessageEol; 503 if(saveNumber != numberIterations_) { 504 problemStatus_=2; // factorize now 505 rowArray_[1]>clear(); 506 pivotRow_=1; // say no weights update 507 break; 508 } else { 509 // take on more relaxed criterion 510 if (saveDj*dualIn_<1.0e20 511 fabs(saveDjdualIn_)>1.0e4*(1.0+fabs(dualIn_))) { 512 // need to reject something 513 char x = isColumn(sequenceIn_) ? 'C' :'R'; 514 handler_>message(CLP_SIMPLEX_FLAG,messages_) 515 <<x<<sequenceWithin(sequenceIn_) 516 <<OsiMessageEol; 517 setFlagged(sequenceIn_); 518 lastBadIteration_ = numberIterations_; // say be more cautious 519 rowArray_[1]>clear(); 520 pivotRow_=1; 521 continue; 522 } 523 } 524 } 525 if (pivotRow_>=0) { 526 // if stable replace in basis 527 int updateStatus = factorization_>replaceColumn(rowArray_[2], 528 pivotRow_, 529 alpha_); 530 if (updateStatus==1) { 531 // slight error 532 if (factorization_>pivots()>5) 533 problemStatus_=2; // factorize now 534 } else if (updateStatus==2) { 535 // major error 536 // later we may need to unwind more e.g. fake bounds 537 if(saveNumber != numberIterations_) { 538 problemStatus_=2; // factorize now 539 break; 540 } else { 541 // need to reject something 542 char x = isColumn(sequenceIn_) ? 'C' :'R'; 543 handler_>message(CLP_SIMPLEX_FLAG,messages_) 544 <<x<<sequenceWithin(sequenceIn_) 545 <<OsiMessageEol; 546 setFlagged(sequenceIn_); 547 lastBadIteration_ = numberIterations_; // say be more cautious 548 rowArray_[1]>clear(); 549 pivotRow_=1; 550 continue; 551 } 552 } else if (updateStatus==3) { 553 // out of memory 554 // increase space if not many iterations 555 if (factorization_>pivots()< 556 0.5*factorization_>maximumPivots()&& 557 factorization_>pivots()<200) 558 factorization_>areaFactor( 559 factorization_>areaFactor() * 1.1); 560 problemStatus_=2; // factorize now 561 } 562 // here do part of steepest  ready for next iteration 563 primalColumnPivot_>updateWeights(rowArray_[1]); 564 } else { 565 if (pivotRow_==1) { 566 // no outgoing row is valid 567 #ifdef CLP_DEBUG 568 if (handler_>logLevel()&32) 569 printf("** no row pivot\n"); 570 #endif 571 if (!factorization_>pivots()) { 572 problemStatus_=5; //say looks unbounded 573 // do ray 574 delete [] ray_; 575 ray_ = new double [numberColumns_]; 576 CoinFillN(ray_,numberColumns_,0.0); 577 int number=rowArray_[1]>getNumElements(); 578 int * index = rowArray_[1]>getIndices(); 579 double * array = rowArray_[1]>denseVector(); 580 double way=directionIn_; 581 int i; 582 double zeroTolerance=1.0e12; 583 if (sequenceIn_<numberColumns_) 584 ray_[sequenceIn_]=directionIn_; 585 for (i=0;i<number;i++) { 586 int iRow=index[i]; 587 int iPivot=pivotVariable_[iRow]; 588 double arrayValue = array[iRow]; 589 if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance) 590 ray_[iPivot] = way* array[iRow]; 591 } 592 } 593 rowArray_[0]>clear(); 594 break; 595 } else { 596 // flipping from bound to bound 597 } 598 } 599 600 // update primal solution 601 602 double objectiveChange=0.0; 603 // Cost on pivot row may change  may need to change dualIn 604 double oldCost=0.0; 605 if (pivotRow_>=0) 606 oldCost = cost(pivotVariable_[pivotRow_]); 607 // rowArray_[1] is not empty  used to update djs 608 updatePrimalsInPrimal(rowArray_[1],theta_, objectiveChange); 609 if (pivotRow_>=0) 610 dualIn_ += (oldCostcost(pivotVariable_[pivotRow_])); 611 612 int whatNext=housekeeping(objectiveChange); 613 614 if (whatNext==1) { 615 problemStatus_ =2; // refactorize 616 } else if (whatNext==2) { 617 // maximum iterations or equivalent 618 problemStatus_= 3; 619 break; 620 } 621 } else { 622 // no pivot column 623 #ifdef CLP_DEBUG 624 if (handler_>logLevel()&32) 625 printf("** no column pivot\n"); 626 #endif 627 if (nonLinearCost_>numberInfeasibilities()) 628 problemStatus_=4; // might be infeasible 629 break; 630 } 631 } 655 632 } 656 633 /* Checks if finished. Updates status */
Note: See TracChangeset
for help on using the changeset viewer.