Changeset 1311


Ignore:
Timestamp:
Dec 3, 2008 11:53:28 AM (11 years ago)
Author:
forrest
Message:

To allow solving a structured model
(This is just proof of concept initial coding)

Location:
trunk/Clp
Files:
1 added
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/Clp/src/ClpSimplex.hpp

    r1288 r1311  
    2222class ClpNonLinearCost;
    2323class ClpNodeStuff;
    24 class CoinModel;
     24class CoinStructuredModel;
    2525class OsiClpSolverInterface;
    2626class CoinWarmStartBasis;
     
    293293      =1 use solution */
    294294  int reducedGradient(int phase=0);
     295  /// Solve using structure of model and maybe in parallel
     296  int solve(CoinStructuredModel * model);
     297  /** This loads a model from a CoinStructuredModel object - returns number of errors.
     298      If originalOrder then keep to order stored in blocks,
     299      otherwise first column/rows correspond to first block - etc.
     300      If keepSolution true and size is same as current then
     301      keeps current status and solution
     302  */
     303  int loadProblem (  CoinStructuredModel & modelObject,
     304                     bool originalOrder=true,bool keepSolution=false);
    295305  /**
    296306     When scaling is on it is possible that the scaled problem
  • trunk/Clp/src/ClpSolve.cpp

    r1307 r1311  
    31093109  return matched;
    31103110}
     3111#include "CoinStructuredModel.hpp"
     3112// Solve using structure of model and maybe in parallel
     3113int
     3114ClpSimplex::solve(CoinStructuredModel * model)
     3115{
     3116  // analyze structure
     3117  int numberColumns = model->numberColumns();
     3118  int numberRowBlocks=model->numberRowBlocks();
     3119  int numberColumnBlocks = model->numberColumnBlocks();
     3120  int numberElementBlocks = model->numberElementBlocks();
     3121  if (numberRowBlocks==numberColumnBlocks||numberRowBlocks==numberColumnBlocks+1) {
     3122    // looks like Dantzig-Wolfe
     3123    bool masterColumns = (numberColumnBlocks==numberRowBlocks);
     3124    if ((masterColumns&&numberElementBlocks==2*numberRowBlocks-1)
     3125        ||(!masterColumns&&numberElementBlocks==2*numberRowBlocks)) {
     3126      // make up problems
     3127      int numberBlocks=numberRowBlocks-1;
     3128      // Find master rows and columns
     3129      int * rowCounts = new int [numberRowBlocks];
     3130      CoinZeroN(rowCounts,numberRowBlocks);
     3131      int * columnCounts = new int [numberColumnBlocks+1];
     3132      CoinZeroN(columnCounts,numberColumnBlocks);
     3133      int iBlock;
     3134      for (iBlock = 0;iBlock<numberElementBlocks;iBlock++) {
     3135        int iRowBlock = model->blockType(iBlock).rowBlock;
     3136        rowCounts[iRowBlock]++;
     3137        int iColumnBlock =model->blockType(iBlock).columnBlock;
     3138        columnCounts[iColumnBlock]++;
     3139      }
     3140      int * whichBlock = new int [numberElementBlocks];
     3141      int masterRowBlock=-1;
     3142      for (iBlock = 0;iBlock<numberRowBlocks;iBlock++) {
     3143        if (rowCounts[iBlock]>1) {
     3144          if (masterRowBlock==-1) {
     3145            masterRowBlock=iBlock;
     3146          } else {
     3147            // Can't decode
     3148            masterRowBlock=-2;
     3149            break;
     3150          }
     3151        }
     3152      }
     3153      int masterColumnBlock=-1;
     3154      int kBlock=0;
     3155      for (iBlock = 0;iBlock<numberColumnBlocks;iBlock++) {
     3156        int count=columnCounts[iBlock];
     3157        columnCounts[iBlock]=kBlock;
     3158        kBlock += count;
     3159      }
     3160      for (iBlock = 0;iBlock<numberElementBlocks;iBlock++) {
     3161        int iColumnBlock = model->blockType(iBlock).columnBlock;
     3162        whichBlock[columnCounts[iColumnBlock]]=iBlock;
     3163        columnCounts[iColumnBlock]++;
     3164      }
     3165      for (iBlock = numberColumnBlocks-1;iBlock>=0;iBlock--)
     3166        columnCounts[iBlock+1]=columnCounts[iBlock];
     3167      columnCounts[0]=0;
     3168      for (iBlock = 0;iBlock<numberColumnBlocks;iBlock++) {
     3169        int count=columnCounts[iBlock+1]-columnCounts[iBlock];
     3170        if (count==1) {
     3171          int kBlock = whichBlock[columnCounts[iBlock]];
     3172          int iRowBlock = model->blockType(kBlock).rowBlock;
     3173          if (iRowBlock==masterRowBlock) {
     3174            if (masterColumnBlock==-1) {
     3175              masterColumnBlock=iBlock;
     3176            } else {
     3177              // Can't decode
     3178              masterColumnBlock=-2;
     3179              break;
     3180            }
     3181          }
     3182        }
     3183      }
     3184      if (masterRowBlock<0||masterColumnBlock==-2) {
     3185        // What now
     3186        abort();
     3187      }
     3188      delete [] rowCounts;
     3189      // create all data
     3190      const CoinPackedMatrix ** top = new const CoinPackedMatrix * [numberColumnBlocks];
     3191      ClpSimplex * sub = new ClpSimplex [numberBlocks];
     3192      ClpSimplex master;
     3193      kBlock=0;
     3194      int masterBlock=-1;
     3195      for (iBlock = 0;iBlock<numberColumnBlocks;iBlock++) {
     3196        top[kBlock]=NULL;
     3197        int start=columnCounts[iBlock];
     3198        int end=columnCounts[iBlock+1];
     3199        assert (end-start<=2);
     3200        for (int j=start;j<end;j++) {
     3201          int jBlock = whichBlock[j];
     3202          int iRowBlock = model->blockType(jBlock).rowBlock;
     3203          int iColumnBlock =model->blockType(jBlock).columnBlock;
     3204          assert (iColumnBlock==iBlock);
     3205          if (iColumnBlock!=masterColumnBlock&&iRowBlock==masterRowBlock) {
     3206            // top matrix
     3207            top[kBlock]=model->block(jBlock)->packedMatrix();
     3208          } else {
     3209            const CoinPackedMatrix * matrix
     3210              = model->block(jBlock)->packedMatrix();
     3211            // Get pointers to arrays
     3212            const double * rowLower;
     3213            const double * rowUpper;
     3214            const double * columnLower;
     3215            const double * columnUpper;
     3216            const double * objective;
     3217            model->block(iRowBlock,iColumnBlock,rowLower,rowUpper,
     3218                         columnLower,columnUpper,objective);
     3219            if (iColumnBlock!=masterColumnBlock) {
     3220              // diagonal block
     3221              sub[kBlock].loadProblem(*matrix,columnLower,columnUpper,
     3222                                objective,rowLower,rowUpper);
     3223              // Set columnCounts to be diagonal block index for cleanup
     3224              columnCounts[kBlock]=jBlock;
     3225            } else {
     3226              // master
     3227              masterBlock=jBlock;
     3228              master.loadProblem(*matrix,columnLower,columnUpper,
     3229                                objective,rowLower,rowUpper);
     3230            }
     3231          }
     3232        }
     3233        if (iBlock!=masterColumnBlock)
     3234          kBlock++;
     3235      }
     3236      delete [] whichBlock;
     3237      // For now master must have been defined (does not have to have columns)
     3238      assert (master.numberRows());
     3239      assert (masterBlock>=0);
     3240      // Overkill in terms of space
     3241      int numberMasterRows = master.numberRows();
     3242      int * columnAdd = new int[numberBlocks+1];
     3243      int * rowAdd = new int[numberBlocks*(numberMasterRows+1)];
     3244      double * elementAdd = new double[numberBlocks*(numberMasterRows+1)];
     3245      double * objective = new double[numberBlocks];
     3246      int maxPass=500;
     3247      int iPass;
     3248      double lastObjective=1.0e31;
     3249      // Create convexity rows for proposals
     3250      int numberMasterColumns = master.numberColumns();
     3251      master.resize(numberMasterRows+numberBlocks,numberMasterColumns);
     3252      // Arrays to say which block and when created
     3253      int maximumColumns = 2*numberMasterRows+10*numberBlocks;
     3254      whichBlock = new int[maximumColumns];
     3255      int * when = new int[maximumColumns];
     3256      int numberColumnsGenerated=numberBlocks;
     3257      // fill in rhs and add in artificials
     3258      {
     3259        double * rowLower = master.rowLower();
     3260        double * rowUpper = master.rowUpper();
     3261        int iBlock;
     3262        columnAdd[0] = 0;
     3263        for (iBlock=0;iBlock<numberBlocks;iBlock++) {
     3264          int iRow = iBlock + numberMasterRows;;
     3265          rowLower[iRow]=1.0;
     3266          rowUpper[iRow]=1.0;
     3267          rowAdd[iBlock] = iRow;
     3268          elementAdd[iBlock] = 1.0;
     3269          objective[iBlock] = 1.0e9;
     3270          columnAdd[iBlock+1] = iBlock+1;
     3271          when[iBlock]=-1;
     3272          whichBlock[iBlock] = iBlock;
     3273        }
     3274        master.addColumns(numberBlocks,NULL,NULL,objective,
     3275                          columnAdd, rowAdd, elementAdd);
     3276      }
     3277      // and resize matrix to double check clp will be happy
     3278      //master.matrix()->setDimensions(numberMasterRows+numberBlocks,
     3279      //                         numberMasterColumns+numberBlocks);
     3280     
     3281      for (iPass=0;iPass<maxPass;iPass++) {
     3282        printf("Start of pass %d\n",iPass);
     3283        // Solve master - may be infeasible
     3284        master.scaling(false);
     3285        if (0) {
     3286          master.writeMps("yy.mps");
     3287        }
     3288       
     3289        master.primal(1);
     3290        int problemStatus = master.status(); // do here as can change (delcols)
     3291        if (master.numberIterations()==0&&iPass)
     3292          break; // finished
     3293        if (master.objectiveValue()>lastObjective-1.0e-7&&iPass>555)
     3294          break; // finished
     3295        lastObjective = master.objectiveValue();
     3296        // mark basic ones and delete if necessary
     3297        int iColumn;
     3298        numberColumnsGenerated=master.numberColumns()-numberMasterColumns;
     3299        for (iColumn=0;iColumn<numberColumnsGenerated;iColumn++) {
     3300          if (master.getStatus(iColumn+numberMasterColumns)==ClpSimplex::basic)
     3301            when[iColumn]=iPass;
     3302        }
     3303        if (numberColumnsGenerated+numberBlocks>maximumColumns) {
     3304          // delete
     3305          int numberKeep=0;
     3306          int numberDelete=0;
     3307          int * whichDelete = new int[numberColumnsGenerated];
     3308          for (iColumn=0;iColumn<numberColumnsGenerated;iColumn++) {
     3309            if (when[iColumn]>iPass-7) {
     3310              // keep
     3311              when[numberKeep]=when[iColumn];
     3312              whichBlock[numberKeep++]=whichBlock[iColumn];
     3313            } else {
     3314              // delete
     3315              whichDelete[numberDelete++]=iColumn+numberMasterColumns;
     3316            }
     3317          }
     3318          numberColumnsGenerated -= numberDelete;
     3319          master.deleteColumns(numberDelete,whichDelete);
     3320          delete [] whichDelete;
     3321        }
     3322        const double * dual=NULL;
     3323        bool deleteDual=false;
     3324        if (problemStatus==0) {
     3325          dual = master.dualRowSolution();
     3326        } else if (problemStatus==1) {
     3327          // could do composite objective
     3328          dual = master.infeasibilityRay();
     3329          deleteDual = true;
     3330          printf("The sum of infeasibilities is %g\n",
     3331                 master.sumPrimalInfeasibilities());
     3332        } else if (!master.numberColumns()) {
     3333          assert(!iPass);
     3334          dual = master.dualRowSolution();
     3335          memset(master.dualRowSolution(),
     3336                 0,(numberMasterRows+numberBlocks)*sizeof(double));
     3337        } else {
     3338          abort();
     3339        }
     3340        // Create objective for sub problems and solve
     3341        columnAdd[0]=0;
     3342        int numberProposals=0;
     3343        for (iBlock=0;iBlock<numberBlocks;iBlock++) {
     3344          int numberColumns2 = sub[iBlock].numberColumns();
     3345          double * saveObj = new double [numberColumns2];
     3346          double * objective2 = sub[iBlock].objective();
     3347          memcpy(saveObj,objective2,numberColumns2*sizeof(double));
     3348          // new objective
     3349          top[iBlock]->transposeTimes(dual,objective2);
     3350          int i;
     3351          if (problemStatus==0) {
     3352            for (i=0;i<numberColumns2;i++)
     3353              objective2[i] = saveObj[i]-objective2[i];
     3354          } else {
     3355            for (i=0;i<numberColumns2;i++)
     3356              objective2[i] = -objective2[i];
     3357          }
     3358          if (iPass)
     3359            sub[iBlock].primal();
     3360          else
     3361            sub[iBlock].dual();
     3362          memcpy(objective2,saveObj,numberColumns2*sizeof(double));
     3363          // get proposal
     3364          if (sub[iBlock].numberIterations()||!iPass) {
     3365            double objValue =0.0;
     3366            int start = columnAdd[numberProposals];
     3367            // proposal
     3368            if (sub[iBlock].isProvenOptimal()) {
     3369              const double * solution = sub[iBlock].primalColumnSolution();
     3370              top[iBlock]->times(solution,elementAdd+start);
     3371              for (i=0;i<numberColumns2;i++)
     3372                objValue += solution[i]*saveObj[i];
     3373              // See if good dj and pack down
     3374              int number=start;
     3375              double dj = objValue;
     3376              if (problemStatus)
     3377                dj=0.0;
     3378              double smallest=1.0e100;
     3379              double largest=0.0;
     3380              for (i=0;i<numberMasterRows;i++) {
     3381                double value = elementAdd[start+i];
     3382                if (fabs(value)>1.0e-15) {
     3383                  dj -= dual[i]*value;
     3384                  smallest = CoinMin(smallest,fabs(value));
     3385                  largest = CoinMax(largest,fabs(value));
     3386                  rowAdd[number]=i;
     3387                  elementAdd[number++]=value;
     3388                }
     3389              }
     3390              // and convexity
     3391              dj -= dual[numberMasterRows+iBlock];
     3392              rowAdd[number]=numberMasterRows+iBlock;
     3393              elementAdd[number++]=1.0;
     3394              // if elements large then scale?
     3395              //if (largest>1.0e8||smallest<1.0e-8)
     3396              printf("For subproblem %d smallest - %g, largest %g - dj %g\n",
     3397                     iBlock,smallest,largest,dj);
     3398              if (dj<-1.0e-6||!iPass) {
     3399                // take
     3400                objective[numberProposals]=objValue;
     3401                columnAdd[++numberProposals]=number;
     3402                when[numberColumnsGenerated]=iPass;
     3403                whichBlock[numberColumnsGenerated++]=iBlock;
     3404              }
     3405            } else if (sub[iBlock].isProvenDualInfeasible()) {
     3406              // use ray
     3407              const double * solution = sub[iBlock].unboundedRay();
     3408              top[iBlock]->times(solution,elementAdd+start);
     3409              for (i=0;i<numberColumns2;i++)
     3410                objValue += solution[i]*saveObj[i];
     3411              // See if good dj and pack down
     3412              int number=start;
     3413              double dj = objValue;
     3414              double smallest=1.0e100;
     3415              double largest=0.0;
     3416              for (i=0;i<numberMasterRows;i++) {
     3417                double value = elementAdd[start+i];
     3418                if (fabs(value)>1.0e-15) {
     3419                  dj -= dual[i]*value;
     3420                  smallest = CoinMin(smallest,fabs(value));
     3421                  largest = CoinMax(largest,fabs(value));
     3422                  rowAdd[number]=i;
     3423                  elementAdd[number++]=value;
     3424                }
     3425              }
     3426              // if elements large or small then scale?
     3427              //if (largest>1.0e8||smallest<1.0e-8)
     3428              printf("For subproblem ray %d smallest - %g, largest %g - dj %g\n",
     3429                     iBlock,smallest,largest,dj);
     3430              if (dj<-1.0e-6) {
     3431                // take
     3432                objective[numberProposals]=objValue;
     3433                columnAdd[++numberProposals]=number;
     3434                when[numberColumnsGenerated]=iPass;
     3435                whichBlock[numberColumnsGenerated++]=iBlock;
     3436              }
     3437            } else {
     3438              abort();
     3439            }
     3440          }
     3441          delete [] saveObj;
     3442        }
     3443        if (deleteDual)
     3444          delete [] dual;
     3445        if (numberProposals)
     3446          master.addColumns(numberProposals,NULL,NULL,objective,
     3447                            columnAdd,rowAdd,elementAdd);
     3448      }
     3449      //master.scaling(false);
     3450      //master.primal(1);
     3451      loadProblem(*model);
     3452      // now put back a good solution
     3453      double * lower = new double[numberMasterRows+numberBlocks];
     3454      double * upper = new double[numberMasterRows+numberBlocks];
     3455      numberColumnsGenerated  += numberMasterColumns;
     3456      double * sol = new double[numberColumnsGenerated];
     3457      const double * solution = master.primalColumnSolution();
     3458      const double * masterLower = master.rowLower();
     3459      const double * masterUpper = master.rowUpper();
     3460      double * fullSolution = primalColumnSolution();
     3461      const double * fullLower = columnLower();
     3462      const double * fullUpper = columnUpper();
     3463      const double * rowSolution = master.primalRowSolution();
     3464      double * fullRowSolution = primalRowSolution();
     3465      const int * rowBack = model->block(masterBlock)->originalRows();
     3466      int numberRows2 = model->block(masterBlock)->numberRows();
     3467      const int * columnBack = model->block(masterBlock)->originalColumns();
     3468      int numberColumns2 = model->block(masterBlock)->numberColumns();
     3469      for (int iRow=0;iRow<numberRows2;iRow++) {
     3470        int kRow = rowBack[iRow];
     3471        setRowStatus(kRow,master.getRowStatus(iRow));
     3472        fullRowSolution[kRow]=rowSolution[iRow];
     3473      }
     3474      for (int iColumn=0;iColumn<numberColumns2;iColumn++) {
     3475        int kColumn = columnBack[iColumn];
     3476        setStatus(kColumn,master.getStatus(iColumn));
     3477        fullSolution[kColumn]=solution[iColumn];
     3478      }
     3479      for (iBlock=0;iBlock<numberBlocks;iBlock++) {
     3480        // convert top bit to by rows
     3481        CoinPackedMatrix topMatrix = *top[iBlock];
     3482        topMatrix.reverseOrdering();
     3483        // zero solution
     3484        memset(sol,0,numberColumnsGenerated*sizeof(double));
     3485       
     3486        for (int i=numberMasterColumns;i<numberColumnsGenerated;i++) {
     3487          if (whichBlock[i-numberMasterColumns]==iBlock)
     3488            sol[i] = solution[i];
     3489        }
     3490        memset(lower,0,(numberMasterRows+numberBlocks)*sizeof(double));
     3491        master.times(1.0,sol,lower);
     3492        for (int iRow=0;iRow<numberMasterRows;iRow++) {
     3493          double value = lower[iRow];
     3494          if (masterUpper[iRow]<1.0e20)
     3495            upper[iRow] = value;
     3496          else
     3497            upper[iRow]=COIN_DBL_MAX;
     3498          if (masterLower[iRow]>-1.0e20)
     3499            lower[iRow] = value;
     3500          else
     3501            lower[iRow]=-COIN_DBL_MAX;
     3502        }
     3503        sub[iBlock].addRows(numberMasterRows,lower,upper,
     3504                            topMatrix.getVectorStarts(),
     3505                            topMatrix.getVectorLengths(),
     3506                            topMatrix.getIndices(),
     3507                            topMatrix.getElements());
     3508        sub[iBlock].primal(1);
     3509        const double * subSolution = sub[iBlock].primalColumnSolution();
     3510        const double * subRowSolution = sub[iBlock].primalRowSolution();
     3511        // move solution
     3512        int kBlock = columnCounts[iBlock];
     3513        const int * rowBack = model->block(kBlock)->originalRows();
     3514        int numberRows2 = model->block(kBlock)->numberRows();
     3515        const int * columnBack = model->block(kBlock)->originalColumns();
     3516        int numberColumns2 = model->block(kBlock)->numberColumns();
     3517        for (int iRow=0;iRow<numberRows2;iRow++) {
     3518          int kRow = rowBack[iRow];
     3519          setRowStatus(kRow,sub[iBlock].getRowStatus(iRow));
     3520          fullRowSolution[kRow]=subRowSolution[iRow];
     3521        }
     3522        for (int iColumn=0;iColumn<numberColumns2;iColumn++) {
     3523          int kColumn = columnBack[iColumn];
     3524          setStatus(kColumn,sub[iBlock].getStatus(iColumn));
     3525          fullSolution[kColumn]=subSolution[iColumn];
     3526        }
     3527      }
     3528      for (int iColumn=0;iColumn<numberColumns;iColumn++) {
     3529        if (fullSolution[iColumn]<fullUpper[iColumn]-1.0e-8&&
     3530            fullSolution[iColumn]>fullLower[iColumn]+1.0e-8) {
     3531          assert(getStatus(iColumn)==ClpSimplex::basic);
     3532        } else if (fullSolution[iColumn]>=fullUpper[iColumn]-1.0e-8) {
     3533          // may help to make rest non basic
     3534          setStatus(iColumn,ClpSimplex::atUpperBound);
     3535        } else if (fullSolution[iColumn]<=fullLower[iColumn]+1.0e-8) {
     3536          // may help to make rest non basic
     3537          setStatus(iColumn,ClpSimplex::atLowerBound);
     3538        }
     3539      }
     3540      //int numberRows=model->numberRows();
     3541      //for (int iRow=0;iRow<numberRows;iRow++)
     3542      //setRowStatus(iRow,ClpSimplex::superBasic);
     3543      primal(1);
     3544      delete [] columnCounts;
     3545      delete [] sol;
     3546      delete [] lower;
     3547      delete [] upper;
     3548      delete [] whichBlock;
     3549      delete [] when;
     3550      delete [] columnAdd;
     3551      delete [] rowAdd;
     3552      delete [] elementAdd;
     3553      delete [] objective;
     3554      delete [] top;
     3555      delete [] sub;
     3556    } else {
     3557      printf("What structure - look further?\n");
     3558      loadProblem(*model);
     3559      abort();
     3560    }
     3561  } else {
     3562    printf("Does not look decomposable????\n");
     3563    loadProblem(*model);
     3564  }
     3565  primal(1);
     3566  return 0;
     3567}
     3568/* This loads a model from a CoinStructuredModel object - returns number of errors.
     3569   If originalOrder then keep to order stored in blocks,
     3570   otherwise first column/rows correspond to first block - etc.
     3571   If keepSolution true and size is same as current then
     3572   keeps current status and solution
     3573*/
     3574int
     3575ClpSimplex::loadProblem (  CoinStructuredModel & coinModel,
     3576                           bool originalOrder,
     3577                           bool keepSolution)
     3578{
     3579  unsigned char * status = NULL;
     3580  double * psol = NULL;
     3581  double * dsol = NULL;
     3582  int numberRows=coinModel.numberRows();
     3583  int numberColumns = coinModel.numberColumns();
     3584  int numberRowBlocks=coinModel.numberRowBlocks();
     3585  int numberColumnBlocks = coinModel.numberColumnBlocks();
     3586  int numberElementBlocks = coinModel.numberElementBlocks();
     3587  if (status_&&numberRows_&&numberRows_==numberRows&&
     3588      numberColumns_==numberColumns&&keepSolution) {
     3589    status = new unsigned char [numberRows_+numberColumns_];
     3590    CoinMemcpyN(status_,numberRows_+numberColumns_,status);
     3591    psol = new double [numberRows_+numberColumns_];
     3592    CoinMemcpyN(columnActivity_,numberColumns_,psol);
     3593    CoinMemcpyN(rowActivity_,numberRows_,psol+numberColumns_);
     3594    dsol = new double [numberRows_+numberColumns_];
     3595    CoinMemcpyN(reducedCost_,numberColumns_,dsol);
     3596    CoinMemcpyN(dual_,numberRows_,dsol+numberColumns_);
     3597  }
     3598  int returnCode=0;
     3599  double * rowLower = new double [numberRows];
     3600  double * rowUpper = new double [numberRows];
     3601  double * columnLower = new double [numberColumns];
     3602  double * columnUpper = new double [numberColumns];
     3603  double * objective = new double [numberColumns];
     3604  int * integerType = new int [numberColumns];
     3605  CoinBigIndex numberElements=0;
     3606  // Bases for blocks
     3607  int * rowBase = new int[numberRowBlocks];
     3608  CoinFillN(rowBase,numberRowBlocks,-1);
     3609  // And row to put it
     3610  int * whichRow = new int [numberRows];
     3611  int * columnBase = new int[numberColumnBlocks];
     3612  CoinFillN(columnBase,numberColumnBlocks,-1);
     3613  // And column to put it
     3614  int * whichColumn = new int [numberColumns];
     3615  for (int iBlock=0;iBlock<numberElementBlocks;iBlock++) {
     3616    CoinModel * block = coinModel.block(iBlock);
     3617    numberElements += block->numberElements();
     3618    //and set up elements etc
     3619    double * associated = block->associatedArray();
     3620    // If strings then do copies
     3621    if (block->stringsExist())
     3622      returnCode += block->createArrays(rowLower, rowUpper, columnLower, columnUpper,
     3623                                          objective, integerType,associated);
     3624    const CoinModelBlockInfo & info = coinModel.blockType(iBlock);
     3625    int iRowBlock = info.rowBlock;
     3626    int iColumnBlock = info.columnBlock;
     3627    if (rowBase[iRowBlock]<0) {
     3628      rowBase[iRowBlock]=block->numberRows();
     3629      // Save block number
     3630      whichRow[numberRows-numberRowBlocks+iRowBlock]= iBlock;
     3631    } else {
     3632      assert(rowBase[iRowBlock]==block->numberRows());
     3633    }
     3634    if (columnBase[iColumnBlock]<0) {
     3635      columnBase[iColumnBlock]=block->numberColumns();
     3636      // Save block number
     3637      whichColumn[numberColumns-numberColumnBlocks+iColumnBlock]= iBlock;
     3638    } else {
     3639      assert(columnBase[iColumnBlock]==block->numberColumns());
     3640    }
     3641  }
     3642  // Fill arrays with defaults
     3643  CoinFillN(rowLower,numberRows,-COIN_DBL_MAX);
     3644  CoinFillN(rowUpper,numberRows,COIN_DBL_MAX);
     3645  CoinFillN(columnLower,numberColumns,0.0);
     3646  CoinFillN(columnUpper,numberColumns,COIN_DBL_MAX);
     3647  CoinFillN(objective,numberColumns,0.0);
     3648  CoinFillN(integerType,numberColumns,0);
     3649  int n=0;
     3650  for (int iBlock=0;iBlock<numberRowBlocks;iBlock++) {
     3651    int k = rowBase[iBlock];
     3652    rowBase[iBlock]=n;
     3653    assert (k>=0);
     3654    // block number
     3655    int jBlock = whichRow[numberRows-numberRowBlocks+iBlock];
     3656    if (originalOrder) {
     3657      memcpy(whichRow+n,coinModel.block(jBlock)->originalRows(),k*sizeof(int));
     3658    } else {
     3659      CoinIotaN(whichRow+n,k,n);
     3660    }
     3661    n+=k;
     3662  }
     3663  assert (n==numberRows);
     3664  n=0;
     3665  for (int iBlock=0;iBlock<numberColumnBlocks;iBlock++) {
     3666    int k = columnBase[iBlock];
     3667    columnBase[iBlock]=n;
     3668    assert (k>=0);
     3669    // block number
     3670    int jBlock = whichColumn[numberColumns-numberColumnBlocks+iBlock];
     3671    if (originalOrder) {
     3672      memcpy(whichColumn+n,coinModel.block(jBlock)->originalColumns(),
     3673             k*sizeof(int));
     3674    } else {
     3675      CoinIotaN(whichColumn+n,k,n);
     3676    }
     3677    n+=k;
     3678  }
     3679  assert (n==numberColumns);
     3680  bool gotIntegers=false;
     3681  for (int iBlock=0;iBlock<numberElementBlocks;iBlock++) {
     3682    CoinModel * block = coinModel.block(iBlock);
     3683    const CoinModelBlockInfo & info = coinModel.blockType(iBlock);
     3684    int iRowBlock = info.rowBlock;
     3685    int iRowBase = rowBase[iRowBlock];
     3686    int iColumnBlock = info.columnBlock;
     3687    int iColumnBase = columnBase[iColumnBlock];
     3688    if (info.rhs) {
     3689      int nRows = block->numberRows();
     3690      const double * lower = block->rowLowerArray();
     3691      const double * upper = block->rowUpperArray();
     3692      for (int i=0;i<nRows;i++) {
     3693        int put = whichRow[i+iRowBase];
     3694        rowLower[put] = lower[i];
     3695        rowUpper[put] = upper[i];
     3696      }
     3697    }
     3698    if (info.bounds) {
     3699      int nColumns = block->numberColumns();
     3700      const double * lower = block->columnLowerArray();
     3701      const double * upper = block->columnUpperArray();
     3702      const double * obj = block->objectiveArray();
     3703      for (int i=0;i<nColumns;i++) {
     3704        int put = whichColumn[i+iColumnBase];
     3705        columnLower[put] = lower[i];
     3706        columnUpper[put] = upper[i];
     3707        objective[put] = obj[i];
     3708      }
     3709    }
     3710    if (info.integer) {
     3711      gotIntegers=true;
     3712      int nColumns = block->numberColumns();
     3713      const int * type = block->integerTypeArray();
     3714      for (int i=0;i<nColumns;i++) {
     3715        int put = whichColumn[i+iColumnBase];
     3716        integerType[put] = type[i];
     3717      }
     3718    }
     3719  }
     3720  gutsOfLoadModel(numberRows, numberColumns,
     3721                  columnLower, columnUpper, objective, rowLower, rowUpper, NULL);
     3722  delete [] rowLower;
     3723  delete [] rowUpper;
     3724  delete [] columnLower;
     3725  delete [] columnUpper;
     3726  delete [] objective;
     3727  // Do integers if wanted
     3728  if (gotIntegers) {
     3729    for (int iColumn=0;iColumn<numberColumns;iColumn++) {
     3730      if (integerType[iColumn])
     3731        setInteger(iColumn);
     3732    }
     3733  }
     3734  delete [] integerType;
     3735  setObjectiveOffset(coinModel.objectiveOffset());
     3736  // Space for elements
     3737  int * row = new int[numberElements];
     3738  int * column = new int[numberElements];
     3739  double * element = new double[numberElements];
     3740  numberElements=0;
     3741  for (int iBlock=0;iBlock<numberElementBlocks;iBlock++) {
     3742    CoinModel * block = coinModel.block(iBlock);
     3743    const CoinModelBlockInfo & info = coinModel.blockType(iBlock);
     3744    int iRowBlock = info.rowBlock;
     3745    int iRowBase = rowBase[iRowBlock];
     3746    int iColumnBlock = info.columnBlock;
     3747    int iColumnBase = columnBase[iColumnBlock];
     3748    if (info.rowName) {
     3749      int numberItems = block->rowNames()->numberItems();
     3750      assert( block->numberRows()>=numberItems);
     3751      if (numberItems) {
     3752        const char *const * rowNames=block->rowNames()->names();
     3753        for (int i=0;i<numberItems;i++) {
     3754          int put = whichRow[i+iRowBase];
     3755          std::string name = rowNames[i];
     3756          setRowName(put,name);
     3757        }
     3758      }
     3759    }
     3760    if (info.columnName) {
     3761      int numberItems = block->columnNames()->numberItems();
     3762      assert( block->numberColumns()>=numberItems);
     3763      if (numberItems) {
     3764        const char *const * columnNames=block->columnNames()->names();
     3765        for (int i=0;i<numberItems;i++) {
     3766          int put = whichColumn[i+iColumnBase];
     3767          std::string name = columnNames[i];
     3768          setColumnName(put,name);
     3769        }
     3770      }
     3771    }
     3772    if (info.matrix) {
     3773      CoinPackedMatrix matrix2;
     3774      const CoinPackedMatrix * matrix = block->packedMatrix();
     3775      if (!matrix) {
     3776        double * associated = block->associatedArray();
     3777        block->createPackedMatrix(matrix2,associated);
     3778        matrix = &matrix2;
     3779      }
     3780      // get matrix data pointers
     3781      const int * row2 = matrix->getIndices();
     3782      const CoinBigIndex * columnStart = matrix->getVectorStarts();
     3783      const double * elementByColumn = matrix->getElements();
     3784      const int * columnLength = matrix->getVectorLengths();
     3785      int n = matrix->getNumCols();
     3786      assert (matrix->isColOrdered());
     3787      for (int iColumn=0;iColumn<n;iColumn++) {
     3788        CoinBigIndex j;
     3789        int jColumn = whichColumn[iColumn+iColumnBase];
     3790        for (j=columnStart[iColumn];
     3791             j<columnStart[iColumn]+columnLength[iColumn];j++) {
     3792          row[numberElements]=whichRow[row2[j]+iRowBase];
     3793          column[numberElements]=jColumn;
     3794          element[numberElements++]=elementByColumn[j];
     3795        }
     3796      }
     3797    }
     3798  }
     3799  delete [] whichRow;
     3800  delete [] whichColumn;
     3801  delete [] rowBase;
     3802  delete [] columnBase;
     3803  CoinPackedMatrix * matrix =
     3804    new CoinPackedMatrix (true,row,column,element,numberElements);
     3805  matrix_ = new ClpPackedMatrix(matrix);
     3806  matrix_->setDimensions(numberRows,numberColumns);
     3807  delete [] row;
     3808  delete [] column;
     3809  delete [] element;
     3810  createStatus();
     3811  if (status) {
     3812    // copy back
     3813    CoinMemcpyN(status,numberRows_+numberColumns_,status_);
     3814    CoinMemcpyN(psol,numberColumns_,columnActivity_);
     3815    CoinMemcpyN(psol+numberColumns_,numberRows_,rowActivity_);
     3816    CoinMemcpyN(dsol,numberColumns_,reducedCost_);
     3817    CoinMemcpyN(dsol+numberColumns_,numberRows_,dual_);
     3818    delete [] status;
     3819    delete [] psol;
     3820    delete [] dsol;
     3821  }
     3822  optimizationDirection_ = coinModel.optimizationDirection(); 
     3823  return returnCode;
     3824}
Note: See TracChangeset for help on using the changeset viewer.