Changeset 1311 for trunk/Clp/src/ClpSolve.cpp
 Timestamp:
 Dec 3, 2008 11:53:28 AM (11 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/Clp/src/ClpSolve.cpp
r1307 r1311 3109 3109 return matched; 3110 3110 } 3111 #include "CoinStructuredModel.hpp" 3112 // Solve using structure of model and maybe in parallel 3113 int 3114 ClpSimplex::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==numberColumnBlocksnumberRowBlocks==numberColumnBlocks+1) { 3122 // looks like DantzigWolfe 3123 bool masterColumns = (numberColumnBlocks==numberRowBlocks); 3124 if ((masterColumns&&numberElementBlocks==2*numberRowBlocks1) 3125 (!masterColumns&&numberElementBlocks==2*numberRowBlocks)) { 3126 // make up problems 3127 int numberBlocks=numberRowBlocks1; 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 = numberColumnBlocks1;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<0masterColumnBlock==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 (endstart<=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()>lastObjective1.0e7&&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]>iPass7) { 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.0e15) { 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.0e8smallest<1.0e8) 3396 printf("For subproblem %d smallest  %g, largest %g  dj %g\n", 3397 iBlock,smallest,largest,dj); 3398 if (dj<1.0e6!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.0e15) { 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.0e8smallest<1.0e8) 3428 printf("For subproblem ray %d smallest  %g, largest %g  dj %g\n", 3429 iBlock,smallest,largest,dj); 3430 if (dj<1.0e6) { 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[inumberMasterColumns]==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.0e8&& 3530 fullSolution[iColumn]>fullLower[iColumn]+1.0e8) { 3531 assert(getStatus(iColumn)==ClpSimplex::basic); 3532 } else if (fullSolution[iColumn]>=fullUpper[iColumn]1.0e8) { 3533 // may help to make rest non basic 3534 setStatus(iColumn,ClpSimplex::atUpperBound); 3535 } else if (fullSolution[iColumn]<=fullLower[iColumn]+1.0e8) { 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 */ 3574 int 3575 ClpSimplex::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[numberRowsnumberRowBlocks+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[numberColumnsnumberColumnBlocks+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[numberRowsnumberRowBlocks+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[numberColumnsnumberColumnBlocks+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.