Changeset 700


Ignore:
Timestamp:
Dec 2, 2005 6:46:13 AM (14 years ago)
Author:
forrest
Message:

for quadratic barrier

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/ClpQuadraticObjective.cpp

    r468 r700  
    337337                if (iColumn!=jColumn) {
    338338                  offset += valueI*valueJ*elementValue;
     339                  //if (fabs(valueI*valueJ*elementValue)>1.0e-12)
     340                  //printf("%d %d %g %g %g -> %g\n",
     341                  //       iColumn,jColumn,valueI,valueJ,elementValue,
     342                  //       valueI*valueJ*elementValue);
    339343                  double gradientI = valueJ*elementValue;
    340344                  double gradientJ = valueI*elementValue;
     
    343347                } else {
    344348                  offset += 0.5*valueI*valueI*elementValue;
     349                  //if (fabs(valueI*valueI*elementValue)>1.0e-12)
     350                  //printf("XX %d %g %g -> %g\n",
     351                  //       iColumn,valueI,elementValue,
     352                  //       0.5*valueI*valueI*elementValue);
    345353                  double gradientI = valueI*elementValue;
    346354                  gradient_[iColumn] += gradientI;
     
    369377        }
    370378      }
    371       offset *= model->optimizationDirection()*model->objectiveScale();
     379      if (model)
     380        offset *= model->optimizationDirection()*model->objectiveScale();
    372381      return gradient_;
    373382    }
     
    462471      }
    463472    }
    464     offset *= model->optimizationDirection();
     473    if (model)
     474      offset *= model->optimizationDirection();
    465475    return gradient_;
    466476  }
  • trunk/ClpSimplex.cpp

    r698 r700  
    46714671  ClpInterior barrier;
    46724672  barrier.borrowModel(*model2);
    4673 // Preference is WSSMP, UFL, TAUCS then base
     4673  // See if quadratic objective
     4674  ClpQuadraticObjective * quadraticObj = NULL;
     4675  if (objective_->type()==2)
     4676    quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
     4677  // If Quadratic we need KKT
     4678  bool doKKT = (quadraticObj!=NULL);
     4679  // Preference is WSSMP, UFL, TAUCS then base
    46744680#ifdef WSSMP_BARRIER
    4675   ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp(CoinMax(100,model2->numberRows()/10));
     4681 if (!doKKT) {
     4682   ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp(CoinMax(100,model2->numberRows()/10));
     4683   barrier.setCholesky(cholesky);
     4684 } else {
    46764685  //ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp();
    4677   //ClpCholeskyWssmpKKT * cholesky = new ClpCholeskyWssmpKKT(CoinMax(100,model2->numberRows()/10));
     4686   ClpCholeskyWssmpKKT * cholesky = new ClpCholeskyWssmpKKT(CoinMax(100,model2->numberRows()/10));
     4687   barrier.setCholesky(cholesky);
     4688 }
    46784689#elif UFL_BARRIER
    4679   ClpCholeskyUfl * cholesky = new ClpCholeskyUfl();
     4690 if (!doKKT) {
     4691   ClpCholeskyUfl * cholesky = new ClpCholeskyUfl();
     4692   barrier.setCholesky(cholesky);
     4693 } else {
     4694   ClpCholeskyUfl * cholesky = new ClpCholeskyUfl();
     4695   cholesky->setKKT(true);
     4696   barrier.setCholesky(cholesky);
     4697 }
    46804698#elif TAUCS_BARRIER
    4681   ClpCholeskyTaucs * cholesky = new ClpCholeskyTaucs();
     4699 assert (!doKKT);
     4700 ClpCholeskyTaucs * cholesky = new ClpCholeskyTaucs();
     4701 barrier.setCholesky(cholesky);
    46824702#else
     4703 if (!doKKT) {
    46834704  ClpCholeskyBase * cholesky = new ClpCholeskyBase();
     4705  barrier.setCholesky(cholesky);
     4706 } else {
     4707   ClpCholeskyBase * cholesky = new ClpCholeskyBase();
     4708   cholesky->setKKT(true);
     4709   barrier.setCholesky(cholesky);
     4710 }
    46844711#endif
    4685   barrier.setCholesky(cholesky);
     4712  barrier.setDiagonalPerturbation(1.0e-14);
    46864713  int numberRows = model2->numberRows();
    46874714  int numberColumns = model2->numberColumns();
     
    63896416    if (updateStatus) {
    63906417      updateSolution=false;
    6391       returnCode=-1;
     6418      returnCode=1;
    63926419    }
    63936420    // if no pivots, bad update but reasonable alpha - take and invert
     
    63986425      // slight error
    63996426      if (factorization_->pivots()>5||updateStatus==4) {
    6400         returnCode=-1;
     6427        returnCode=1;
    64016428      }
    64026429    } else if (updateStatus==2) {
     
    64266453        roundAgain=true;
    64276454      } else {
    6428         returnCode=1;
     6455        returnCode=-1;
    64296456      }
    64306457    } else if (updateStatus==3) {
     
    64366463        factorization_->areaFactor(
    64376464                                   factorization_->areaFactor() * 1.1);
    6438       returnCode =-1; // factorize now
     6465      returnCode =1; // factorize now
    64396466    }
    64406467    rowArray_[1]->clear();
     
    64476474    }
    64486475  }
    6449   if (returnCode == -1) {
     6476  if (returnCode == 1) {
    64506477    // refactorize here
    64516478    int factorStatus = internalFactorize(1);
  • trunk/Test/unitTest.cpp

    r687 r700  
    1919#include "ClpSimplex.hpp"
    2020#include "ClpSimplexOther.hpp"
     21#include "ClpSimplexNonlinear.hpp"
    2122#include "ClpInterior.hpp"
    2223#include "ClpLinearObjective.hpp"
     
    185186  return solveOptions;
    186187}
     188static void printSol(ClpSimplex & model)
     189{
     190  int numberRows = model.numberRows();
     191  int numberColumns = model.numberColumns();
     192 
     193  double * rowPrimal = model.primalRowSolution();
     194  double * rowDual = model.dualRowSolution();
     195  double * rowLower = model.rowLower();
     196  double * rowUpper = model.rowUpper();
     197 
     198  int iRow;
     199  double objValue = model.getObjValue();
     200  printf("Objvalue %g Rows (%d)\n",objValue,numberRows);
     201  for (iRow=0;iRow<numberRows;iRow++) {
     202    printf("%d primal %g dual %g low %g up %g\n",
     203           iRow,rowPrimal[iRow],rowDual[iRow],
     204           rowLower[iRow],rowUpper[iRow]);
     205  }
     206  double * columnPrimal = model.primalColumnSolution();
     207  double * columnDual = model.dualColumnSolution();
     208  double * columnLower = model.columnLower();
     209  double * columnUpper = model.columnUpper();
     210  double offset;
     211  //const double * gradient = model.objectiveAsObject()->gradient(&model,
     212  //                                                       columnPrimal,offset,true,1);
     213  const double * gradient = model.objective(columnPrimal,offset);
     214  int iColumn;
     215  objValue = -offset-model.objectiveOffset();
     216  printf("offset %g (%g)\n",offset,model.objectiveOffset());
     217  printf("Columns (%d)\n",numberColumns);
     218  for (iColumn=0;iColumn<numberColumns;iColumn++) {
     219    printf("%d primal %g dual %g low %g up %g\n",
     220           iColumn,columnPrimal[iColumn],columnDual[iColumn],
     221           columnLower[iColumn],columnUpper[iColumn]);
     222    objValue += columnPrimal[iColumn]*gradient[iColumn];
     223    if (fabs(columnPrimal[iColumn]*gradient[iColumn])>1.0e-8)
     224      printf("obj -> %g gradient %g\n",objValue,gradient[iColumn]);
     225  }
     226  printf("Computed objective %g\n",objValue);
     227}
    187228//----------------------------------------------------------------
    188229// unitTest [-mpsDir=V1] [-netlibDir=V2] [-test]
     
    16161657    delete [] element;
    16171658    //solution.quadraticSLP(50,1.0e-4);
    1618     double objValue = solution.getObjValue();
    16191659    CoinRelFltEq eq(1.0e-4);
    16201660    //assert(eq(objValue,820.0));
    16211661    //solution.setLogLevel(63);
    16221662    solution.primal();
    1623     int numberRows = solution.numberRows();
    1624 
    1625     double * rowPrimal = solution.primalRowSolution();
    1626     double * rowDual = solution.dualRowSolution();
    1627     double * rowLower = solution.rowLower();
    1628     double * rowUpper = solution.rowUpper();
    1629    
    1630     int iRow;
    1631     printf("Rows\n");
    1632     for (iRow=0;iRow<numberRows;iRow++) {
    1633       printf("%d primal %g dual %g low %g up %g\n",
    1634              iRow,rowPrimal[iRow],rowDual[iRow],
    1635              rowLower[iRow],rowUpper[iRow]);
    1636     }
    1637     double * columnPrimal = solution.primalColumnSolution();
    1638     double * columnDual = solution.dualColumnSolution();
    1639     double * columnLower = solution.columnLower();
    1640     double * columnUpper = solution.columnUpper();
    1641     objValue = solution.getObjValue();
    1642     int iColumn;
    1643     printf("Columns\n");
    1644     for (iColumn=0;iColumn<numberColumns;iColumn++) {
    1645       printf("%d primal %g dual %g low %g up %g\n",
    1646              iColumn,columnPrimal[iColumn],columnDual[iColumn],
    1647              columnLower[iColumn],columnUpper[iColumn]);
    1648     }
     1663    printSol(solution);
    16491664    //assert(eq(objValue,3.2368421));
    16501665    //exit(77);
     
    16961711    delete [] element;
    16971712    int numberColumns=model.numberColumns();
     1713    model.scaling(0);
    16981714#if 0
    16991715    model.nonlinearSLP(50,1.0e-4);
     
    17111727    model.setFactorizationFrequency(10);
    17121728    model.primal();
     1729    printSol(model);
    17131730    double objValue = model.getObjValue();
    17141731    CoinRelFltEq eq(1.0e-4);
    17151732    assert(eq(objValue,-400.92));
     1733    // and again for barrier
     1734    model.barrier(false);
     1735    //printSol(model);
     1736    model.allSlackBasis();
     1737    model.primal();
     1738    //printSol(model);
    17161739  }
    17171740  if (0) {   
Note: See TracChangeset for help on using the changeset viewer.