Changeset 636


Ignore:
Timestamp:
Jun 8, 2005 12:22:14 PM (15 years ago)
Author:
forrest
Message:

for faster mip

Location:
trunk
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • trunk/ClpEventHandler.cpp

    r344 r636  
    5656ClpEventHandler::event(Event whichEvent)
    5757{
    58   return -1; // do nothing
     58  if (whichEvent!=theta)
     59    return -1; // do nothing
     60  else
     61    return 0; // say normal exit
    5962}
    6063/* set model. */
  • trunk/ClpSimplex.cpp

    r628 r636  
    8787  rowActivityWork_(NULL),
    8888  columnActivityWork_(NULL),
     89  auxiliaryModel_(NULL),
    8990  numberDualInfeasibilities_(0),
    9091  numberDualInfeasibilitiesWithoutFree_(0),
     
    192193  rowActivityWork_(NULL),
    193194  columnActivityWork_(NULL),
     195  auxiliaryModel_(NULL),
    194196  numberDualInfeasibilities_(0),
    195197  numberDualInfeasibilitiesWithoutFree_(0),
     
    16151617  rowActivityWork_(NULL),
    16161618  columnActivityWork_(NULL),
     1619  auxiliaryModel_(NULL),
    16171620  numberDualInfeasibilities_(0),
    16181621  numberDualInfeasibilitiesWithoutFree_(0),
     
    17141717  rowActivityWork_(NULL),
    17151718  columnActivityWork_(NULL),
     1719  auxiliaryModel_(NULL),
    17161720  numberDualInfeasibilities_(0),
    17171721  numberDualInfeasibilitiesWithoutFree_(0),
     
    17781782  maximumBasic_ = rhs.maximumBasic_;
    17791783  int numberRows2 = numberRows_+numberExtraRows_;
     1784  auxiliaryModel_ = NULL;
    17801785  if ((whatsChanged_&1)!=0) {
    17811786    lower_ = ClpCopyOfArray(rhs.lower_,numberColumns_+numberRows2);
     
    19611966  if (!type) {
    19621967    // delete everything
     1968    delete auxiliaryModel_;
     1969    auxiliaryModel_ = NULL;
    19631970    delete factorization_;
    19641971    factorization_ = NULL;
     
    24862493  }
    24872494  bool oldMatrix = ((startFinishOptions&4)!=0&&(whatsChanged_&1)!=0);
     2495  if (auxiliaryModel_)
     2496    oldMatrix=true;
    24882497  if (what==63) {
    24892498    if (!status_)
     
    25252534  }
    25262535  int numberRows2 = numberRows_+numberExtraRows_;
    2527   //int numberTotal = numberRows_+numberColumns_;
     2536  int numberTotal = numberRows2+numberColumns_;
    25282537  int i;
    25292538  bool doSanityCheck=true;
    2530   if (what==63) {
     2539  if (what==63&&!auxiliaryModel_) {
    25312540    // We may want to switch stuff off for speed
    25322541    if ((specialOptions_&256)!=0)
     
    26822691      delete [] cost_;
    26832692      // extra copy with original costs
    2684       int nTotal = numberRows2+numberColumns_;
    2685       //cost_ = new double[2*nTotal];
    2686       cost_ = new double[nTotal];
     2693      //cost_ = new double[2*numberTotal];
     2694      cost_ = new double[numberTotal];
    26872695      delete [] lower_;
    26882696      delete [] upper_;
     
    26912699      delete [] dj_;
    26922700      dj_ = new double[numberRows2+numberColumns_];
    2693       reducedCostWork_ = dj_;
    2694       rowReducedCost_ = dj_+numberColumns_;
    26952701      delete [] solution_;
    26962702      solution_ = new double[numberRows2+numberColumns_];
    2697       columnActivityWork_ = solution_;
    2698       rowActivityWork_ = solution_+numberColumns_;
    2699     }
     2703    } else if (auxiliaryModel_) {
     2704      assert (!cost_);
     2705      cost_=auxiliaryModel_->cost_;
     2706      assert (!lower_);
     2707      lower_=auxiliaryModel_->lower_;
     2708      assert (!upper_);
     2709      upper_=auxiliaryModel_->upper_;
     2710      assert (!dj_);
     2711      dj_=auxiliaryModel_->dj_;
     2712      assert (!solution_);
     2713      solution_=auxiliaryModel_->solution_;
     2714      assert (!rowScale_);
     2715      assert (!columnScale_);
     2716    }
     2717    reducedCostWork_ = dj_;
     2718    rowReducedCost_ = dj_+numberColumns_;
     2719    columnActivityWork_ = solution_;
     2720    rowActivityWork_ = solution_+numberColumns_;
    27002721    objectiveWork_ = cost_;
    27012722    rowObjectiveWork_ = cost_+numberColumns_;
     
    27062727  }
    27072728  if ((what&4)!=0) {
    2708     double direction = optimizationDirection_*objectiveScale_;
    2709     const double * obj = objective();
    2710     // and also scale by scale factors
    2711     if (rowScale_) {
    2712       if (rowObjective_) {
    2713         for (i=0;i<numberRows_;i++)
    2714           rowObjectiveWork_[i] = rowObjective_[i]*direction/rowScale_[i];
     2729    if (!auxiliaryModel_||(what==63&&(auxiliaryModel_->specialOptions_&4)==0)) {
     2730      double direction = optimizationDirection_*objectiveScale_;
     2731      const double * obj = objective();
     2732      const double * rowScale = auxiliaryModel_ ? auxiliaryModel_->rowScale_ : rowScale_;
     2733      const double * columnScale = auxiliaryModel_ ? auxiliaryModel_->columnScale_ : columnScale_;
     2734      // and also scale by scale factors
     2735      if (rowScale) {
     2736        if (rowObjective_) {
     2737          for (i=0;i<numberRows_;i++)
     2738            rowObjectiveWork_[i] = rowObjective_[i]*direction/rowScale[i];
     2739        } else {
     2740          memset(rowObjectiveWork_,0,numberRows_*sizeof(double));
     2741        }
     2742        // If scaled then do all columns later in one loop
     2743        if (what!=63||auxiliaryModel_) {
     2744          for (i=0;i<numberColumns_;i++) {
     2745            CoinAssert(fabs(obj[i])<1.0e25);
     2746            objectiveWork_[i] = obj[i]*direction*columnScale[i];
     2747          }
     2748        }
    27152749      } else {
    2716         memset(rowObjectiveWork_,0,numberRows_*sizeof(double));
    2717       }
    2718       // If scaled then do all columns later in one loop
    2719       if (what!=63||!rowScale_) {
     2750        if (rowObjective_) {
     2751          for (i=0;i<numberRows_;i++)
     2752            rowObjectiveWork_[i] = rowObjective_[i]*direction;
     2753        } else {
     2754          memset(rowObjectiveWork_,0,numberRows_*sizeof(double));
     2755        }
    27202756        for (i=0;i<numberColumns_;i++) {
    27212757          CoinAssert(fabs(obj[i])<1.0e25);
    2722           objectiveWork_[i] = obj[i]*direction*columnScale_[i];
     2758          objectiveWork_[i] = obj[i]*direction;
    27232759        }
    27242760      }
     2761      if (auxiliaryModel_) {
     2762        // save costs
     2763        memcpy(auxiliaryModel_->cost_+numberTotal,cost_,numberTotal*sizeof(double));
     2764      }
    27252765    } else {
    2726       if (rowObjective_) {
    2727         for (i=0;i<numberRows_;i++)
    2728           rowObjectiveWork_[i] = rowObjective_[i]*direction;
    2729       } else {
    2730         memset(rowObjectiveWork_,0,numberRows_*sizeof(double));
    2731       }
    2732       for (i=0;i<numberColumns_;i++) {
    2733         CoinAssert(fabs(obj[i])<1.0e25);
    2734         objectiveWork_[i] = obj[i]*direction;
    2735       }
     2766      // just copy
     2767      memcpy(cost_,auxiliaryModel_->cost_+numberTotal,numberTotal*sizeof(double));
    27362768    }
    27372769  }
    27382770  if ((what&1)!=0) {
     2771    const double * rowScale = auxiliaryModel_ ? auxiliaryModel_->rowScale_ : rowScale_;
     2772    const double * columnScale = auxiliaryModel_ ? auxiliaryModel_->columnScale_ : columnScale_;
    27392773    // clean up any mismatches on infinity
    27402774    // and fix any variables with tiny gaps
    27412775    double primalTolerance=dblParam_[ClpPrimalTolerance];
    2742     if(rowScale_) {
    2743       // If scaled then do all columns later in one loop
    2744       if (what!=63||!rowScale_) {
     2776    if (!auxiliaryModel_) {
     2777      if(rowScale) {
     2778        // If scaled then do all columns later in one loop
     2779        if (what!=63) {
     2780          for (i=0;i<numberColumns_;i++) {
     2781            double multiplier = rhsScale_/columnScale[i];
     2782            double lowerValue = columnLower_[i];
     2783            double upperValue = columnUpper_[i];
     2784            if (lowerValue>-1.0e20) {
     2785              columnLowerWork_[i]=lowerValue*multiplier;
     2786              if (upperValue>=1.0e20) {
     2787                columnUpperWork_[i]=COIN_DBL_MAX;
     2788              } else {
     2789                columnUpperWork_[i]=upperValue*multiplier;
     2790                if (columnUpperWork_[i]-columnLowerWork_[i]<=primalTolerance) {
     2791                  if (columnLowerWork_[i]>=0.0) {
     2792                    columnUpperWork_[i] = columnLowerWork_[i];
     2793                  } else if (columnUpperWork_[i]<=0.0) {
     2794                    columnLowerWork_[i] = columnUpperWork_[i];
     2795                  } else {
     2796                    columnUpperWork_[i] = 0.0;
     2797                    columnLowerWork_[i] = 0.0;
     2798                  }
     2799                }
     2800              }
     2801            } else if (upperValue<1.0e20) {
     2802              columnLowerWork_[i]=-COIN_DBL_MAX;
     2803              columnUpperWork_[i]=upperValue*multiplier;
     2804            } else {
     2805              // free
     2806              columnLowerWork_[i]=-COIN_DBL_MAX;
     2807              columnUpperWork_[i]=COIN_DBL_MAX;
     2808            }
     2809          }
     2810        }
     2811        for (i=0;i<numberRows_;i++) {
     2812          double multiplier = rhsScale_*rowScale[i];
     2813          double lowerValue = rowLower_[i];
     2814          double upperValue = rowUpper_[i];
     2815          if (lowerValue>-1.0e20) {
     2816            rowLowerWork_[i]=lowerValue*multiplier;
     2817            if (upperValue>=1.0e20) {
     2818              rowUpperWork_[i]=COIN_DBL_MAX;
     2819            } else {
     2820              rowUpperWork_[i]=upperValue*multiplier;
     2821              if (rowUpperWork_[i]-rowLowerWork_[i]<=primalTolerance) {
     2822                if (rowLowerWork_[i]>=0.0) {
     2823                  rowUpperWork_[i] = rowLowerWork_[i];
     2824                } else if (rowUpperWork_[i]<=0.0) {
     2825                  rowLowerWork_[i] = rowUpperWork_[i];
     2826                } else {
     2827                  rowUpperWork_[i] = 0.0;
     2828                  rowLowerWork_[i] = 0.0;
     2829                }
     2830              }
     2831            }
     2832          } else if (upperValue<1.0e20) {
     2833            rowLowerWork_[i]=-COIN_DBL_MAX;
     2834            rowUpperWork_[i]=upperValue*multiplier;
     2835          } else {
     2836            // free
     2837            rowLowerWork_[i]=-COIN_DBL_MAX;
     2838            rowUpperWork_[i]=COIN_DBL_MAX;
     2839          }
     2840        }
     2841      } else if (rhsScale_!=1.0) {
    27452842        for (i=0;i<numberColumns_;i++) {
    2746           double multiplier = rhsScale_/columnScale_[i];
     2843          double lowerValue = columnLower_[i];
     2844          double upperValue = columnUpper_[i];
     2845          if (lowerValue>-1.0e20) {
     2846            columnLowerWork_[i]=lowerValue*rhsScale_;
     2847            if (upperValue>=1.0e20) {
     2848              columnUpperWork_[i]=COIN_DBL_MAX;
     2849            } else {
     2850              columnUpperWork_[i]=upperValue*rhsScale_;
     2851              if (columnUpperWork_[i]-columnLowerWork_[i]<=primalTolerance) {
     2852                if (columnLowerWork_[i]>=0.0) {
     2853                  columnUpperWork_[i] = columnLowerWork_[i];
     2854                } else if (columnUpperWork_[i]<=0.0) {
     2855                  columnLowerWork_[i] = columnUpperWork_[i];
     2856                } else {
     2857                  columnUpperWork_[i] = 0.0;
     2858                  columnLowerWork_[i] = 0.0;
     2859                }
     2860              }
     2861            }
     2862          } else if (upperValue<1.0e20) {
     2863            columnLowerWork_[i]=-COIN_DBL_MAX;
     2864            columnUpperWork_[i]=upperValue*rhsScale_;
     2865          } else {
     2866            // free
     2867            columnLowerWork_[i]=-COIN_DBL_MAX;
     2868            columnUpperWork_[i]=COIN_DBL_MAX;
     2869          }
     2870        }
     2871        for (i=0;i<numberRows_;i++) {
     2872          double lowerValue = rowLower_[i];
     2873          double upperValue = rowUpper_[i];
     2874          if (lowerValue>-1.0e20) {
     2875            rowLowerWork_[i]=lowerValue*rhsScale_;
     2876            if (upperValue>=1.0e20) {
     2877              rowUpperWork_[i]=COIN_DBL_MAX;
     2878            } else {
     2879              rowUpperWork_[i]=upperValue*rhsScale_;
     2880              if (rowUpperWork_[i]-rowLowerWork_[i]<=primalTolerance) {
     2881                if (rowLowerWork_[i]>=0.0) {
     2882                  rowUpperWork_[i] = rowLowerWork_[i];
     2883                } else if (rowUpperWork_[i]<=0.0) {
     2884                  rowLowerWork_[i] = rowUpperWork_[i];
     2885                } else {
     2886                  rowUpperWork_[i] = 0.0;
     2887                  rowLowerWork_[i] = 0.0;
     2888                }
     2889              }
     2890            }
     2891          } else if (upperValue<1.0e20) {
     2892            rowLowerWork_[i]=-COIN_DBL_MAX;
     2893            rowUpperWork_[i]=upperValue*rhsScale_;
     2894          } else {
     2895            // free
     2896            rowLowerWork_[i]=-COIN_DBL_MAX;
     2897            rowUpperWork_[i]=COIN_DBL_MAX;
     2898          }
     2899        }
     2900      } else {
     2901        for (i=0;i<numberColumns_;i++) {
     2902          double lowerValue = columnLower_[i];
     2903          double upperValue = columnUpper_[i];
     2904          if (lowerValue>-1.0e20) {
     2905            columnLowerWork_[i]=lowerValue;
     2906            if (upperValue>=1.0e20) {
     2907              columnUpperWork_[i]=COIN_DBL_MAX;
     2908            } else {
     2909              columnUpperWork_[i]=upperValue;
     2910              if (columnUpperWork_[i]-columnLowerWork_[i]<=primalTolerance) {
     2911                if (columnLowerWork_[i]>=0.0) {
     2912                  columnUpperWork_[i] = columnLowerWork_[i];
     2913                } else if (columnUpperWork_[i]<=0.0) {
     2914                  columnLowerWork_[i] = columnUpperWork_[i];
     2915                } else {
     2916                  columnUpperWork_[i] = 0.0;
     2917                  columnLowerWork_[i] = 0.0;
     2918                }
     2919              }
     2920            }
     2921          } else if (upperValue<1.0e20) {
     2922            columnLowerWork_[i]=-COIN_DBL_MAX;
     2923            columnUpperWork_[i]=upperValue;
     2924          } else {
     2925            // free
     2926            columnLowerWork_[i]=-COIN_DBL_MAX;
     2927            columnUpperWork_[i]=COIN_DBL_MAX;
     2928          }
     2929        }
     2930        for (i=0;i<numberRows_;i++) {
     2931          double lowerValue = rowLower_[i];
     2932          double upperValue = rowUpper_[i];
     2933          if (lowerValue>-1.0e20) {
     2934            rowLowerWork_[i]=lowerValue;
     2935            if (upperValue>=1.0e20) {
     2936              rowUpperWork_[i]=COIN_DBL_MAX;
     2937            } else {
     2938              rowUpperWork_[i]=upperValue;
     2939              if (rowUpperWork_[i]-rowLowerWork_[i]<=primalTolerance) {
     2940                if (rowLowerWork_[i]>=0.0) {
     2941                  rowUpperWork_[i] = rowLowerWork_[i];
     2942                } else if (rowUpperWork_[i]<=0.0) {
     2943                  rowLowerWork_[i] = rowUpperWork_[i];
     2944                } else {
     2945                  rowUpperWork_[i] = 0.0;
     2946                  rowLowerWork_[i] = 0.0;
     2947                }
     2948              }
     2949            }
     2950          } else if (upperValue<1.0e20) {
     2951            rowLowerWork_[i]=-COIN_DBL_MAX;
     2952            rowUpperWork_[i]=upperValue;
     2953          } else {
     2954            // free
     2955            rowLowerWork_[i]=-COIN_DBL_MAX;
     2956            rowUpperWork_[i]=COIN_DBL_MAX;
     2957          }
     2958        }
     2959      }
     2960    } else {
     2961      // auxiliary model
     2962      if (what!=63) {
     2963        // just copy
     2964        memcpy(lower_,auxiliaryModel_->lower_+numberTotal,numberTotal*sizeof(double));
     2965        memcpy(upper_,auxiliaryModel_->upper_+numberTotal,numberTotal*sizeof(double));
     2966      } else {
     2967        assert (rhsScale_==1.0);
     2968        assert (objectiveScale_==1.0);
     2969        if ((auxiliaryModel_->specialOptions_&2)==0) {
     2970          // need to do column bounds
     2971          if (columnScale) {
     2972            const double * inverseScale = columnScale+numberColumns_;
     2973            // scaled
     2974            for (i=0;i<numberColumns_;i++) {
     2975              double value;
     2976              value = columnLower_[i];
     2977              if (value>-1.0e20)
     2978                value *= inverseScale[i];
     2979              lower_[i] = value;
     2980              value = columnUpper_[i];
     2981              if (value<1.0e20)
     2982                value *= inverseScale[i];
     2983              upper_[i] = value;
     2984            }
     2985          } else {
     2986            for (i=0;i<numberColumns_;i++) {
     2987              lower_[i] = columnLower_[i];
     2988              upper_[i] = columnUpper_[i];
     2989            }
     2990          }
     2991          // save
     2992          memcpy(auxiliaryModel_->lower_+numberTotal,lower_,numberColumns_*sizeof(double));
     2993          memcpy(auxiliaryModel_->upper_+numberTotal,upper_,numberColumns_*sizeof(double));
     2994        } else {
     2995          // just copy
     2996          memcpy(lower_,auxiliaryModel_->lower_+numberTotal,numberColumns_*sizeof(double));
     2997          memcpy(upper_,auxiliaryModel_->upper_+numberTotal,numberColumns_*sizeof(double));
     2998        }
     2999        if ((auxiliaryModel_->specialOptions_&1)==0) {
     3000          // need to do row bounds
     3001          if (rowScale) {
     3002            // scaled
     3003            for (i=0;i<numberRows_;i++) {
     3004              double value;
     3005              value = rowLower_[i];
     3006              if (value>-1.0e20)
     3007                value *= rowScale[i];
     3008              lower_[i+numberColumns_] = value;
     3009              value = rowUpper_[i];
     3010              if (value<1.0e20)
     3011                value *= rowScale[i];
     3012              upper_[i+numberColumns_] = value;
     3013            }
     3014          } else {
     3015            for (i=0;i<numberRows_;i++) {
     3016              lower_[i+numberColumns_] = rowLower_[i];
     3017              upper_[i+numberColumns_] = rowUpper_[i];
     3018            }
     3019          }
     3020          // save
     3021          memcpy(auxiliaryModel_->lower_+numberTotal+numberColumns_,
     3022                 lower_+numberColumns_,numberRows_*sizeof(double));
     3023          memcpy(auxiliaryModel_->upper_+numberTotal+numberColumns_,
     3024                 upper_+numberColumns_,numberRows_*sizeof(double));
     3025        } else {
     3026          // just copy
     3027          memcpy(lower_+numberColumns_,
     3028                 auxiliaryModel_->lower_+numberTotal+numberColumns_,numberRows_*sizeof(double));
     3029          memcpy(upper_+numberColumns_,
     3030                 auxiliaryModel_->upper_+numberTotal+numberColumns_,numberRows_*sizeof(double));
     3031        }
     3032      }
     3033      if (what==63) {
     3034        // do basis
     3035        assert ((auxiliaryModel_->specialOptions_&8)!=0);
     3036        // clean solution trusting basis
     3037        for ( i=0;i<numberTotal;i++) {
     3038          Status status =getStatus(i);
     3039          double value=0.0;
     3040          if (status!=basic) {
     3041            if (upper_[i]==lower_[i]) {
     3042              setStatus(i,isFixed);
     3043              value=lower_[i];
     3044            } else if (status==atLowerBound) {
     3045              if (lower_[i]>-1.0e20) {
     3046                value=lower_[i];
     3047              } else {
     3048                if (upper_[i]<1.0e20) {
     3049                  value=upper_[i];
     3050                  setStatus(i,atUpperBound);
     3051                } else {
     3052                  setStatus(i,isFree);
     3053                }
     3054              }
     3055            } else if (status==atUpperBound) {
     3056              if (upper_[i]<1.0e20) {
     3057              value=upper_[i];
     3058              } else {
     3059                if (lower_[i]>-1.0e20) {
     3060                  value=lower_[i];
     3061                  setStatus(i,atLowerBound);
     3062                } else {
     3063                  setStatus(i,isFree);
     3064                }
     3065              }
     3066            }
     3067          }
     3068          solution_[i]=value;
     3069        }
     3070      }
     3071    }
     3072  }
     3073  if (what==63) {
     3074    // move information to work arrays
     3075    double direction = optimizationDirection_;
     3076    // direction is actually scale out not scale in
     3077    if (direction)
     3078      direction = 1.0/direction;
     3079    if (direction!=1.0&&(!auxiliaryModel_||(auxiliaryModel_->specialOptions_&8)==0)) {
     3080      // reverse all dual signs
     3081      for (i=0;i<numberColumns_;i++)
     3082        reducedCost_[i] *= direction;
     3083      for (i=0;i<numberRows_;i++)
     3084        dual_[i] *= direction;
     3085    }
     3086    for (i=0;i<numberRows_+numberColumns_;i++) {
     3087      setFakeBound(i,noFake);
     3088    }
     3089    if (!auxiliaryModel_) {
     3090      if (rowScale_) {
     3091        const double * obj = objective();
     3092        double direction = optimizationDirection_*objectiveScale_;
     3093        // clean up any mismatches on infinity
     3094        // and fix any variables with tiny gaps
     3095        double primalTolerance=dblParam_[ClpPrimalTolerance];
     3096        // on entry
     3097        for (i=0;i<numberColumns_;i++) {
     3098          CoinAssert(fabs(obj[i])<1.0e25);
     3099          double scaleFactor = columnScale_[i];
     3100          double multiplier = rhsScale_/scaleFactor;
     3101          scaleFactor *= direction;
     3102          objectiveWork_[i] = obj[i]*scaleFactor;
     3103          reducedCostWork_[i] = reducedCost_[i]*scaleFactor;
    27473104          double lowerValue = columnLower_[i];
    27483105          double upperValue = columnUpper_[i];
     
    27723129            columnUpperWork_[i]=COIN_DBL_MAX;
    27733130          }
    2774         }
    2775       }
    2776       for (i=0;i<numberRows_;i++) {
    2777         double multiplier = rhsScale_*rowScale_[i];
    2778         double lowerValue = rowLower_[i];
    2779         double upperValue = rowUpper_[i];
    2780         if (lowerValue>-1.0e20) {
    2781           rowLowerWork_[i]=lowerValue*multiplier;
    2782           if (upperValue>=1.0e20) {
    2783             rowUpperWork_[i]=COIN_DBL_MAX;
    2784           } else {
    2785             rowUpperWork_[i]=upperValue*multiplier;
    2786             if (rowUpperWork_[i]-rowLowerWork_[i]<=primalTolerance) {
    2787               if (rowLowerWork_[i]>=0.0) {
    2788                 rowUpperWork_[i] = rowLowerWork_[i];
    2789               } else if (rowUpperWork_[i]<=0.0) {
    2790                 rowLowerWork_[i] = rowUpperWork_[i];
    2791               } else {
    2792                 rowUpperWork_[i] = 0.0;
    2793                 rowLowerWork_[i] = 0.0;
    2794               }
     3131          double value = columnActivity_[i] * multiplier;
     3132          if (fabs(value)>1.0e20) {
     3133            //printf("bad value of %g for column %d\n",value,i);
     3134            setColumnStatus(i,superBasic);
     3135            if (columnUpperWork_[i]<0.0) {
     3136              value=columnUpperWork_[i];
     3137            } else if (columnLowerWork_[i]>0.0) {
     3138              value=columnLowerWork_[i];
     3139            } else {
     3140              value=0.0;
    27953141            }
    27963142          }
    2797         } else if (upperValue<1.0e20) {
    2798           rowLowerWork_[i]=-COIN_DBL_MAX;
    2799           rowUpperWork_[i]=upperValue*multiplier;
    2800         } else {
    2801           // free
    2802           rowLowerWork_[i]=-COIN_DBL_MAX;
    2803           rowUpperWork_[i]=COIN_DBL_MAX;
     3143          columnActivityWork_[i] = value;
    28043144        }
    2805       }
    2806     } else if (rhsScale_!=1.0) {
    2807       for (i=0;i<numberColumns_;i++) {
    2808         double lowerValue = columnLower_[i];
    2809         double upperValue = columnUpper_[i];
    2810         if (lowerValue>-1.0e20) {
    2811           columnLowerWork_[i]=lowerValue*rhsScale_;
    2812           if (upperValue>=1.0e20) {
    2813             columnUpperWork_[i]=COIN_DBL_MAX;
    2814           } else {
    2815             columnUpperWork_[i]=upperValue*rhsScale_;
    2816             if (columnUpperWork_[i]-columnLowerWork_[i]<=primalTolerance) {
    2817               if (columnLowerWork_[i]>=0.0) {
    2818                 columnUpperWork_[i] = columnLowerWork_[i];
    2819               } else if (columnUpperWork_[i]<=0.0) {
    2820                 columnLowerWork_[i] = columnUpperWork_[i];
    2821               } else {
    2822                 columnUpperWork_[i] = 0.0;
    2823                 columnLowerWork_[i] = 0.0;
    2824               }
     3145        for (i=0;i<numberRows_;i++) {
     3146          dual_[i] /= rowScale_[i];
     3147          dual_[i] *= objectiveScale_;
     3148          rowReducedCost_[i] = dual_[i];
     3149          double multiplier = rhsScale_*rowScale_[i];
     3150          double value = rowActivity_[i] * multiplier;
     3151          if (fabs(value)>1.0e20) {
     3152            //printf("bad value of %g for row %d\n",value,i);
     3153            setRowStatus(i,superBasic);
     3154            if (rowUpperWork_[i]<0.0) {
     3155              value=rowUpperWork_[i];
     3156            } else if (rowLowerWork_[i]>0.0) {
     3157              value=rowLowerWork_[i];
     3158            } else {
     3159              value=0.0;
    28253160            }
    28263161          }
    2827         } else if (upperValue<1.0e20) {
    2828           columnLowerWork_[i]=-COIN_DBL_MAX;
    2829           columnUpperWork_[i]=upperValue*rhsScale_;
    2830         } else {
    2831           // free
    2832           columnLowerWork_[i]=-COIN_DBL_MAX;
    2833           columnUpperWork_[i]=COIN_DBL_MAX;
     3162          rowActivityWork_[i] = value;
    28343163        }
    2835       }
    2836       for (i=0;i<numberRows_;i++) {
    2837         double lowerValue = rowLower_[i];
    2838         double upperValue = rowUpper_[i];
    2839         if (lowerValue>-1.0e20) {
    2840           rowLowerWork_[i]=lowerValue*rhsScale_;
    2841           if (upperValue>=1.0e20) {
    2842             rowUpperWork_[i]=COIN_DBL_MAX;
    2843           } else {
    2844             rowUpperWork_[i]=upperValue*rhsScale_;
    2845             if (rowUpperWork_[i]-rowLowerWork_[i]<=primalTolerance) {
    2846               if (rowLowerWork_[i]>=0.0) {
    2847                 rowUpperWork_[i] = rowLowerWork_[i];
    2848               } else if (rowUpperWork_[i]<=0.0) {
    2849                 rowLowerWork_[i] = rowUpperWork_[i];
    2850               } else {
    2851                 rowUpperWork_[i] = 0.0;
    2852                 rowLowerWork_[i] = 0.0;
    2853               }
     3164      } else if (objectiveScale_!=1.0||rhsScale_!=1.0) {
     3165        // on entry
     3166        for (i=0;i<numberColumns_;i++) {
     3167          double value = columnActivity_[i];
     3168          value *= rhsScale_;
     3169          if (fabs(value)>1.0e20) {
     3170            //printf("bad value of %g for column %d\n",value,i);
     3171            setColumnStatus(i,superBasic);
     3172            if (columnUpperWork_[i]<0.0) {
     3173              value=columnUpperWork_[i];
     3174            } else if (columnLowerWork_[i]>0.0) {
     3175              value=columnLowerWork_[i];
     3176            } else {
     3177              value=0.0;
    28543178            }
    28553179          }
    2856         } else if (upperValue<1.0e20) {
    2857           rowLowerWork_[i]=-COIN_DBL_MAX;
    2858           rowUpperWork_[i]=upperValue*rhsScale_;
    2859         } else {
    2860           // free
    2861           rowLowerWork_[i]=-COIN_DBL_MAX;
    2862           rowUpperWork_[i]=COIN_DBL_MAX;
     3180          columnActivityWork_[i] = value;
     3181          reducedCostWork_[i] = reducedCost_[i]*objectiveScale_;
    28633182        }
    2864       }
    2865     } else {
    2866       for (i=0;i<numberColumns_;i++) {
    2867         double lowerValue = columnLower_[i];
    2868         double upperValue = columnUpper_[i];
    2869         if (lowerValue>-1.0e20) {
    2870           columnLowerWork_[i]=lowerValue;
    2871           if (upperValue>=1.0e20) {
    2872             columnUpperWork_[i]=COIN_DBL_MAX;
    2873           } else {
    2874             columnUpperWork_[i]=upperValue;
    2875             if (columnUpperWork_[i]-columnLowerWork_[i]<=primalTolerance) {
    2876               if (columnLowerWork_[i]>=0.0) {
    2877                 columnUpperWork_[i] = columnLowerWork_[i];
    2878               } else if (columnUpperWork_[i]<=0.0) {
    2879                 columnLowerWork_[i] = columnUpperWork_[i];
    2880               } else {
    2881                 columnUpperWork_[i] = 0.0;
    2882                 columnLowerWork_[i] = 0.0;
    2883               }
     3183        for (i=0;i<numberRows_;i++) {
     3184          double value = rowActivity_[i];
     3185          value *= rhsScale_;
     3186          if (fabs(value)>1.0e20) {
     3187            //printf("bad value of %g for row %d\n",value,i);
     3188            setRowStatus(i,superBasic);
     3189            if (rowUpperWork_[i]<0.0) {
     3190              value=rowUpperWork_[i];
     3191            } else if (rowLowerWork_[i]>0.0) {
     3192              value=rowLowerWork_[i];
     3193            } else {
     3194              value=0.0;
    28843195            }
    28853196          }
    2886         } else if (upperValue<1.0e20) {
    2887           columnLowerWork_[i]=-COIN_DBL_MAX;
    2888           columnUpperWork_[i]=upperValue;
    2889         } else {
    2890           // free
    2891           columnLowerWork_[i]=-COIN_DBL_MAX;
    2892           columnUpperWork_[i]=COIN_DBL_MAX;
     3197          rowActivityWork_[i] = value;
     3198          dual_[i] *= objectiveScale_;
     3199          rowReducedCost_[i] = dual_[i];
    28933200        }
    2894       }
    2895       for (i=0;i<numberRows_;i++) {
    2896         double lowerValue = rowLower_[i];
    2897         double upperValue = rowUpper_[i];
    2898         if (lowerValue>-1.0e20) {
    2899           rowLowerWork_[i]=lowerValue;
    2900           if (upperValue>=1.0e20) {
    2901             rowUpperWork_[i]=COIN_DBL_MAX;
    2902           } else {
    2903             rowUpperWork_[i]=upperValue;
    2904             if (rowUpperWork_[i]-rowLowerWork_[i]<=primalTolerance) {
    2905               if (rowLowerWork_[i]>=0.0) {
    2906                 rowUpperWork_[i] = rowLowerWork_[i];
    2907               } else if (rowUpperWork_[i]<=0.0) {
    2908                 rowLowerWork_[i] = rowUpperWork_[i];
    2909               } else {
    2910                 rowUpperWork_[i] = 0.0;
    2911                 rowLowerWork_[i] = 0.0;
    2912               }
     3201      } else {
     3202        // on entry
     3203        for (i=0;i<numberColumns_;i++) {
     3204          double value = columnActivity_[i];
     3205          if (fabs(value)>1.0e20) {
     3206            //printf("bad value of %g for column %d\n",value,i);
     3207            setColumnStatus(i,superBasic);
     3208            if (columnUpperWork_[i]<0.0) {
     3209              value=columnUpperWork_[i];
     3210            } else if (columnLowerWork_[i]>0.0) {
     3211              value=columnLowerWork_[i];
     3212            } else {
     3213              value=0.0;
    29133214            }
    29143215          }
    2915         } else if (upperValue<1.0e20) {
    2916           rowLowerWork_[i]=-COIN_DBL_MAX;
    2917           rowUpperWork_[i]=upperValue;
    2918         } else {
    2919           // free
    2920           rowLowerWork_[i]=-COIN_DBL_MAX;
    2921           rowUpperWork_[i]=COIN_DBL_MAX;
     3216          columnActivityWork_[i] = value;
     3217          reducedCostWork_[i] = reducedCost_[i];
    29223218        }
    2923       }
    2924     }
    2925   }
    2926   if (what==63) {
    2927     // move information to work arrays
    2928     double direction = optimizationDirection_;
    2929     // direction is actually scale out not scale in
    2930     if (direction)
    2931       direction = 1.0/direction;
    2932     if (direction!=1.0) {
    2933       // reverse all dual signs
    2934       for (i=0;i<numberColumns_;i++)
    2935         reducedCost_[i] *= direction;
    2936       for (i=0;i<numberRows_;i++)
    2937         dual_[i] *= direction;
    2938     }
    2939     for (i=0;i<numberRows_+numberColumns_;i++) {
    2940       setFakeBound(i,noFake);
    2941     }
    2942     if (rowScale_) {
    2943       const double * obj = objective();
    2944       double direction = optimizationDirection_*objectiveScale_;
    2945       // clean up any mismatches on infinity
    2946       // and fix any variables with tiny gaps
    2947       double primalTolerance=dblParam_[ClpPrimalTolerance];
    2948       // on entry
    2949       for (i=0;i<numberColumns_;i++) {
    2950         CoinAssert(fabs(obj[i])<1.0e25);
    2951         double scaleFactor = columnScale_[i];
    2952         double multiplier = rhsScale_/scaleFactor;
    2953         scaleFactor *= direction;
    2954         objectiveWork_[i] = obj[i]*scaleFactor;
    2955         reducedCostWork_[i] = reducedCost_[i]*scaleFactor;
    2956         double lowerValue = columnLower_[i];
    2957         double upperValue = columnUpper_[i];
    2958         if (lowerValue>-1.0e20) {
    2959           columnLowerWork_[i]=lowerValue*multiplier;
    2960           if (upperValue>=1.0e20) {
    2961             columnUpperWork_[i]=COIN_DBL_MAX;
    2962           } else {
    2963             columnUpperWork_[i]=upperValue*multiplier;
    2964             if (columnUpperWork_[i]-columnLowerWork_[i]<=primalTolerance) {
    2965               if (columnLowerWork_[i]>=0.0) {
    2966                 columnUpperWork_[i] = columnLowerWork_[i];
    2967               } else if (columnUpperWork_[i]<=0.0) {
    2968                 columnLowerWork_[i] = columnUpperWork_[i];
    2969               } else {
    2970                 columnUpperWork_[i] = 0.0;
    2971                 columnLowerWork_[i] = 0.0;
    2972               }
     3219        for (i=0;i<numberRows_;i++) {
     3220          double value = rowActivity_[i];
     3221          if (fabs(value)>1.0e20) {
     3222            //printf("bad value of %g for row %d\n",value,i);
     3223            setRowStatus(i,superBasic);
     3224            if (rowUpperWork_[i]<0.0) {
     3225              value=rowUpperWork_[i];
     3226            } else if (rowLowerWork_[i]>0.0) {
     3227              value=rowLowerWork_[i];
     3228            } else {
     3229              value=0.0;
    29733230            }
    29743231          }
    2975         } else if (upperValue<1.0e20) {
    2976           columnLowerWork_[i]=-COIN_DBL_MAX;
    2977           columnUpperWork_[i]=upperValue*multiplier;
    2978         } else {
    2979           // free
    2980           columnLowerWork_[i]=-COIN_DBL_MAX;
    2981           columnUpperWork_[i]=COIN_DBL_MAX;
     3232          rowActivityWork_[i] = value;
     3233          rowReducedCost_[i] = dual_[i];
    29823234        }
    2983         double value = columnActivity_[i] * multiplier;
    2984         if (fabs(value)>1.0e20) {
    2985           //printf("bad value of %g for column %d\n",value,i);
    2986           setColumnStatus(i,superBasic);
    2987           if (columnUpperWork_[i]<0.0) {
    2988             value=columnUpperWork_[i];
    2989           } else if (columnLowerWork_[i]>0.0) {
    2990             value=columnLowerWork_[i];
    2991           } else {
    2992             value=0.0;
    2993           }
    2994         }
    2995         columnActivityWork_[i] = value;
    2996       }
    2997       for (i=0;i<numberRows_;i++) {
    2998         dual_[i] /= rowScale_[i];
    2999         dual_[i] *= objectiveScale_;
    3000         rowReducedCost_[i] = dual_[i];
    3001         double multiplier = rhsScale_*rowScale_[i];
    3002         double value = rowActivity_[i] * multiplier;
    3003         if (fabs(value)>1.0e20) {
    3004           //printf("bad value of %g for row %d\n",value,i);
    3005           setRowStatus(i,superBasic);
    3006           if (rowUpperWork_[i]<0.0) {
    3007             value=rowUpperWork_[i];
    3008           } else if (rowLowerWork_[i]>0.0) {
    3009             value=rowLowerWork_[i];
    3010           } else {
    3011             value=0.0;
    3012           }
    3013         }
    3014         rowActivityWork_[i] = value;
    3015       }
    3016     } else if (objectiveScale_!=1.0||rhsScale_!=1.0) {
    3017       // on entry
    3018       for (i=0;i<numberColumns_;i++) {
    3019         double value = columnActivity_[i];
    3020         value *= rhsScale_;
    3021         if (fabs(value)>1.0e20) {
    3022           //printf("bad value of %g for column %d\n",value,i);
    3023           setColumnStatus(i,superBasic);
    3024           if (columnUpperWork_[i]<0.0) {
    3025             value=columnUpperWork_[i];
    3026           } else if (columnLowerWork_[i]>0.0) {
    3027             value=columnLowerWork_[i];
    3028           } else {
    3029             value=0.0;
    3030           }
    3031         }
    3032         columnActivityWork_[i] = value;
    3033         reducedCostWork_[i] = reducedCost_[i]*objectiveScale_;
    3034       }
    3035       for (i=0;i<numberRows_;i++) {
    3036         double value = rowActivity_[i];
    3037         value *= rhsScale_;
    3038         if (fabs(value)>1.0e20) {
    3039           //printf("bad value of %g for row %d\n",value,i);
    3040           setRowStatus(i,superBasic);
    3041           if (rowUpperWork_[i]<0.0) {
    3042             value=rowUpperWork_[i];
    3043           } else if (rowLowerWork_[i]>0.0) {
    3044             value=rowLowerWork_[i];
    3045           } else {
    3046             value=0.0;
    3047           }
    3048         }
    3049         rowActivityWork_[i] = value;
    3050         dual_[i] *= objectiveScale_;
    3051         rowReducedCost_[i] = dual_[i];
    3052       }
    3053     } else {
    3054       // on entry
    3055       for (i=0;i<numberColumns_;i++) {
    3056         double value = columnActivity_[i];
    3057         if (fabs(value)>1.0e20) {
    3058           //printf("bad value of %g for column %d\n",value,i);
    3059           setColumnStatus(i,superBasic);
    3060           if (columnUpperWork_[i]<0.0) {
    3061             value=columnUpperWork_[i];
    3062           } else if (columnLowerWork_[i]>0.0) {
    3063             value=columnLowerWork_[i];
    3064           } else {
    3065             value=0.0;
    3066           }
    3067         }
    3068         columnActivityWork_[i] = value;
    3069         reducedCostWork_[i] = reducedCost_[i];
    3070       }
    3071       for (i=0;i<numberRows_;i++) {
    3072         double value = rowActivity_[i];
    3073         if (fabs(value)>1.0e20) {
    3074           //printf("bad value of %g for row %d\n",value,i);
    3075           setRowStatus(i,superBasic);
    3076           if (rowUpperWork_[i]<0.0) {
    3077             value=rowUpperWork_[i];
    3078           } else if (rowLowerWork_[i]>0.0) {
    3079             value=rowLowerWork_[i];
    3080           } else {
    3081             value=0.0;
    3082           }
    3083         }
    3084         rowActivityWork_[i] = value;
    3085         rowReducedCost_[i] = dual_[i];
    3086       }
    3087     }
    3088   }
    3089  
    3090   if (what==63&&doSanityCheck) {
     3235      }
     3236    }
     3237  }
     3238 
     3239  if (what==63&&doSanityCheck&&!auxiliaryModel_) {
    30913240    // check rim of problem okay
    30923241    if (!sanityCheck())
     
    31553304    }
    31563305  }
    3157 
     3306 
    31583307  if (what==63) {
    31593308    if (newArrays) {
     
    31813330          columnArray_[iColumn]->reserve(CoinMax(numberRows2,numberColumns_));
    31823331      }
     3332    } else {
     3333      int iRow,iColumn;
     3334      if (auxiliaryModel_) {
     3335        for (iRow=0;iRow<4;iRow++) {
     3336          assert (!rowArray_[iRow]);
     3337          rowArray_[iRow]=auxiliaryModel_->rowArray_[iRow];
     3338        }
     3339        for (iColumn=0;iColumn<2;iColumn++) {
     3340          assert (!columnArray_[iColumn]);
     3341          columnArray_[iColumn]=auxiliaryModel_->columnArray_[iColumn];
     3342        }
     3343        // do matrices
     3344        rowCopy_=auxiliaryModel_->rowCopy_;
     3345        ClpMatrixBase * temp = matrix_;
     3346        matrix_=auxiliaryModel_->matrix_;
     3347        auxiliaryModel_->matrix_=temp;
     3348      }
     3349#ifndef NDEBUG
     3350      for (iRow=0;iRow<4;iRow++) {
     3351        int length =numberRows2+factorization_->maximumPivots();
     3352        if (iRow==3||objective_->type()>1)
     3353          length += numberColumns_;
     3354        assert(rowArray_[iRow]->capacity()==length);
     3355        rowArray_[iRow]->checkClear();
     3356      }
     3357     
     3358      for (iColumn=0;iColumn<2;iColumn++) {
     3359        int length =numberColumns_;
     3360        if (iColumn)
     3361          length=CoinMax(numberRows2,numberColumns_);
     3362        assert(columnArray_[iColumn]->capacity()==length);
     3363        columnArray_[iColumn]->checkClear();
     3364      }
     3365#endif
    31833366    }   
    31843367  }
     
    31943377ClpSimplex::deleteRim(int getRidOfFactorizationData)
    31953378{
    3196   int i;
    3197   if (problemStatus_!=1&&problemStatus_!=2) {
    3198     delete [] ray_;
    3199     ray_=NULL;
    3200   }
    3201   // set upperOut_ to furthest away from bound so can use in dual for dualBound_
    3202   upperOut_=1.0;
     3379  if (!auxiliaryModel_) {
     3380    int i;
     3381    if (problemStatus_!=1&&problemStatus_!=2) {
     3382      delete [] ray_;
     3383      ray_=NULL;
     3384    }
     3385    // set upperOut_ to furthest away from bound so can use in dual for dualBound_
     3386    upperOut_=1.0;
    32033387#if 0
    3204   {
    3205     int nBad=0;
    3206     for (i=0;i<numberColumns_;i++) {
    3207       if (lower_[i]==upper_[i]&&getColumnStatus(i)==basic)
    3208         nBad++;
    3209     }
    3210     if (nBad)
    3211       printf("yy %d basic fixed\n",nBad);
    3212   }
     3388    {
     3389      int nBad=0;
     3390      for (i=0;i<numberColumns_;i++) {
     3391        if (lower_[i]==upper_[i]&&getColumnStatus(i)==basic)
     3392          nBad++;
     3393      }
     3394      if (nBad)
     3395        printf("yy %d basic fixed\n",nBad);
     3396    }
    32133397#endif
    3214   // ray may be null if in branch and bound
    3215   if (rowScale_) {
    3216     // Collect infeasibilities
    3217     int numberPrimalScaled=0;
    3218     int numberPrimalUnscaled=0;
    3219     int numberDualScaled=0;
    3220     int numberDualUnscaled=0;
    3221     double scaleC = 1.0/objectiveScale_;
    3222     double scaleR = 1.0/rhsScale_;
    3223     for (i=0;i<numberColumns_;i++) {
    3224       double scaleFactor = columnScale_[i];
    3225       double valueScaled = columnActivityWork_[i];
    3226       double lowerScaled = columnLowerWork_[i];
    3227       double upperScaled = columnUpperWork_[i];
    3228       if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
    3229         if (valueScaled<lowerScaled-primalTolerance_||
    3230             valueScaled>upperScaled+primalTolerance_)
    3231           numberPrimalScaled++;
    3232         else
    3233           upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
    3234       }
    3235       columnActivity_[i] = valueScaled*scaleFactor*scaleR;
    3236       double value = columnActivity_[i];
    3237       if (value<columnLower_[i]-primalTolerance_)
    3238         numberPrimalUnscaled++;
    3239       else if (value>columnUpper_[i]+primalTolerance_)
    3240         numberPrimalUnscaled++;
    3241       double valueScaledDual = reducedCostWork_[i];
    3242       if (valueScaled>columnLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
    3243         numberDualScaled++;
    3244       if (valueScaled<columnUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
    3245         numberDualScaled++;
    3246       reducedCost_[i] = (valueScaledDual*scaleC)/scaleFactor;
    3247       double valueDual = reducedCost_[i];
    3248       if (value>columnLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
    3249         numberDualUnscaled++;
    3250       if (value<columnUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
    3251         numberDualUnscaled++;
    3252     }
    3253     for (i=0;i<numberRows_;i++) {
    3254       double scaleFactor = rowScale_[i];
    3255       double valueScaled = rowActivityWork_[i];
    3256       double lowerScaled = rowLowerWork_[i];
    3257       double upperScaled = rowUpperWork_[i];
    3258       if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
    3259         if (valueScaled<lowerScaled-primalTolerance_||
    3260             valueScaled>upperScaled+primalTolerance_)
    3261           numberPrimalScaled++;
    3262         else
    3263           upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
    3264       }
    3265       rowActivity_[i] = (valueScaled*scaleR)/scaleFactor;
    3266       double value = rowActivity_[i];
    3267       if (value<rowLower_[i]-primalTolerance_)
    3268         numberPrimalUnscaled++;
    3269       else if (value>rowUpper_[i]+primalTolerance_)
    3270         numberPrimalUnscaled++;
    3271       double valueScaledDual = dual_[i]+rowObjectiveWork_[i];;
    3272       if (valueScaled>rowLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
    3273         numberDualScaled++;
    3274       if (valueScaled<rowUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
    3275         numberDualScaled++;
    3276       dual_[i] *= scaleFactor*scaleC;
    3277       double valueDual = dual_[i];
    3278       if (rowObjective_)
    3279         valueDual += rowObjective_[i];
    3280       if (value>rowLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
    3281         numberDualUnscaled++;
    3282       if (value<rowUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
    3283         numberDualUnscaled++;
    3284     }
    3285     if (!problemStatus_&&!secondaryStatus_) {
    3286       // See if we need to set secondary status
    3287       if (numberPrimalUnscaled) {
    3288         if (numberDualUnscaled)
    3289           secondaryStatus_=4;
    3290         else
    3291           secondaryStatus_=2;
     3398    // ray may be null if in branch and bound
     3399    if (rowScale_) {
     3400      // Collect infeasibilities
     3401      int numberPrimalScaled=0;
     3402      int numberPrimalUnscaled=0;
     3403      int numberDualScaled=0;
     3404      int numberDualUnscaled=0;
     3405      double scaleC = 1.0/objectiveScale_;
     3406      double scaleR = 1.0/rhsScale_;
     3407      for (i=0;i<numberColumns_;i++) {
     3408        double scaleFactor = columnScale_[i];
     3409        double valueScaled = columnActivityWork_[i];
     3410        double lowerScaled = columnLowerWork_[i];
     3411        double upperScaled = columnUpperWork_[i];
     3412        if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
     3413          if (valueScaled<lowerScaled-primalTolerance_||
     3414              valueScaled>upperScaled+primalTolerance_)
     3415            numberPrimalScaled++;
     3416          else
     3417            upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
     3418        }
     3419        columnActivity_[i] = valueScaled*scaleFactor*scaleR;
     3420        double value = columnActivity_[i];
     3421        if (value<columnLower_[i]-primalTolerance_)
     3422          numberPrimalUnscaled++;
     3423        else if (value>columnUpper_[i]+primalTolerance_)
     3424          numberPrimalUnscaled++;
     3425        double valueScaledDual = reducedCostWork_[i];
     3426        if (valueScaled>columnLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
     3427          numberDualScaled++;
     3428        if (valueScaled<columnUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
     3429          numberDualScaled++;
     3430        reducedCost_[i] = (valueScaledDual*scaleC)/scaleFactor;
     3431        double valueDual = reducedCost_[i];
     3432        if (value>columnLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
     3433          numberDualUnscaled++;
     3434        if (value<columnUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
     3435          numberDualUnscaled++;
     3436      }
     3437      for (i=0;i<numberRows_;i++) {
     3438        double scaleFactor = rowScale_[i];
     3439        double valueScaled = rowActivityWork_[i];
     3440        double lowerScaled = rowLowerWork_[i];
     3441        double upperScaled = rowUpperWork_[i];
     3442        if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
     3443          if (valueScaled<lowerScaled-primalTolerance_||
     3444              valueScaled>upperScaled+primalTolerance_)
     3445            numberPrimalScaled++;
     3446          else
     3447            upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
     3448        }
     3449        rowActivity_[i] = (valueScaled*scaleR)/scaleFactor;
     3450        double value = rowActivity_[i];
     3451        if (value<rowLower_[i]-primalTolerance_)
     3452          numberPrimalUnscaled++;
     3453        else if (value>rowUpper_[i]+primalTolerance_)
     3454          numberPrimalUnscaled++;
     3455        double valueScaledDual = dual_[i]+rowObjectiveWork_[i];;
     3456        if (valueScaled>rowLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
     3457          numberDualScaled++;
     3458        if (valueScaled<rowUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
     3459          numberDualScaled++;
     3460        dual_[i] *= scaleFactor*scaleC;
     3461        double valueDual = dual_[i];
     3462        if (rowObjective_)
     3463          valueDual += rowObjective_[i];
     3464        if (value>rowLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
     3465          numberDualUnscaled++;
     3466        if (value<rowUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
     3467          numberDualUnscaled++;
     3468      }
     3469      if (!problemStatus_&&!secondaryStatus_) {
     3470        // See if we need to set secondary status
     3471        if (numberPrimalUnscaled) {
     3472          if (numberDualUnscaled)
     3473            secondaryStatus_=4;
     3474          else
     3475            secondaryStatus_=2;
     3476        } else {
     3477          if (numberDualUnscaled)
     3478            secondaryStatus_=3;
     3479        }
     3480      }
     3481      if (problemStatus_==2) {
     3482        for (i=0;i<numberColumns_;i++) {
     3483          ray_[i] *= columnScale_[i];
     3484        }
     3485      } else if (problemStatus_==1&&ray_) {
     3486        for (i=0;i<numberRows_;i++) {
     3487          ray_[i] *= rowScale_[i];
     3488        }
     3489      }
     3490    } else if (rhsScale_!=1.0||objectiveScale_!=1.0) {
     3491      // Collect infeasibilities
     3492      int numberPrimalScaled=0;
     3493      int numberPrimalUnscaled=0;
     3494      int numberDualScaled=0;
     3495      int numberDualUnscaled=0;
     3496      double scaleC = 1.0/objectiveScale_;
     3497      double scaleR = 1.0/rhsScale_;
     3498      for (i=0;i<numberColumns_;i++) {
     3499        double valueScaled = columnActivityWork_[i];
     3500        double lowerScaled = columnLowerWork_[i];
     3501        double upperScaled = columnUpperWork_[i];
     3502        if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
     3503          if (valueScaled<lowerScaled-primalTolerance_||
     3504              valueScaled>upperScaled+primalTolerance_)
     3505            numberPrimalScaled++;
     3506          else
     3507            upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
     3508        }
     3509        columnActivity_[i] = valueScaled*scaleR;
     3510        double value = columnActivity_[i];
     3511        if (value<columnLower_[i]-primalTolerance_)
     3512          numberPrimalUnscaled++;
     3513        else if (value>columnUpper_[i]+primalTolerance_)
     3514          numberPrimalUnscaled++;
     3515        double valueScaledDual = reducedCostWork_[i];
     3516        if (valueScaled>columnLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
     3517          numberDualScaled++;
     3518        if (valueScaled<columnUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
     3519          numberDualScaled++;
     3520        reducedCost_[i] = valueScaledDual*scaleC;
     3521        double valueDual = reducedCost_[i];
     3522        if (value>columnLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
     3523          numberDualUnscaled++;
     3524        if (value<columnUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
     3525          numberDualUnscaled++;
     3526      }
     3527      for (i=0;i<numberRows_;i++) {
     3528        double valueScaled = rowActivityWork_[i];
     3529        double lowerScaled = rowLowerWork_[i];
     3530        double upperScaled = rowUpperWork_[i];
     3531        if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
     3532          if (valueScaled<lowerScaled-primalTolerance_||
     3533              valueScaled>upperScaled+primalTolerance_)
     3534            numberPrimalScaled++;
     3535          else
     3536            upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
     3537        }
     3538        rowActivity_[i] = valueScaled*scaleR;
     3539        double value = rowActivity_[i];
     3540        if (value<rowLower_[i]-primalTolerance_)
     3541          numberPrimalUnscaled++;
     3542        else if (value>rowUpper_[i]+primalTolerance_)
     3543          numberPrimalUnscaled++;
     3544        double valueScaledDual = dual_[i]+rowObjectiveWork_[i];;
     3545        if (valueScaled>rowLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
     3546          numberDualScaled++;
     3547        if (valueScaled<rowUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
     3548          numberDualScaled++;
     3549        dual_[i] *= scaleC;
     3550        double valueDual = dual_[i];
     3551        if (rowObjective_)
     3552          valueDual += rowObjective_[i];
     3553        if (value>rowLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
     3554          numberDualUnscaled++;
     3555        if (value<rowUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
     3556          numberDualUnscaled++;
     3557      }
     3558      if (!problemStatus_&&!secondaryStatus_) {
     3559        // See if we need to set secondary status
     3560        if (numberPrimalUnscaled) {
     3561          if (numberDualUnscaled)
     3562            secondaryStatus_=4;
     3563          else
     3564            secondaryStatus_=2;
     3565        } else {
     3566          if (numberDualUnscaled)
     3567            secondaryStatus_=3;
     3568        }
     3569      }
     3570    } else {
     3571      if (columnActivityWork_) {
     3572        for (i=0;i<numberColumns_;i++) {
     3573          double value = columnActivityWork_[i];
     3574          double lower = columnLowerWork_[i];
     3575          double upper = columnUpperWork_[i];
     3576          if (lower>-1.0e20||upper<1.0e20) {
     3577            if (value>lower&&value<upper)
     3578              upperOut_ = CoinMax(upperOut_,CoinMin(value-lower,upper-value));
     3579          }
     3580          columnActivity_[i] = columnActivityWork_[i];
     3581          reducedCost_[i] = reducedCostWork_[i];
     3582        }
     3583        for (i=0;i<numberRows_;i++) {
     3584          double value = rowActivityWork_[i];
     3585          double lower = rowLowerWork_[i];
     3586          double upper = rowUpperWork_[i];
     3587          if (lower>-1.0e20||upper<1.0e20) {
     3588            if (value>lower&&value<upper)
     3589              upperOut_ = CoinMax(upperOut_,CoinMin(value-lower,upper-value));
     3590          }
     3591          rowActivity_[i] = rowActivityWork_[i];
     3592        }
     3593      }
     3594    }
     3595    // switch off scalefactor if auto
     3596    if (automaticScale_) {
     3597      rhsScale_=1.0;
     3598      objectiveScale_=1.0;
     3599    }
     3600    if (optimizationDirection_!=1.0) {
     3601      // and modify all dual signs
     3602      for (i=0;i<numberColumns_;i++)
     3603        reducedCost_[i] *= optimizationDirection_;
     3604      for (i=0;i<numberRows_;i++)
     3605        dual_[i] *= optimizationDirection_;
     3606    }
     3607  } else {
     3608    // auxiliary model
     3609    cost_=NULL;
     3610    lower_=NULL;
     3611    upper_=NULL;
     3612    dj_=NULL;
     3613    solution_=NULL;
     3614    int iRow,iColumn;
     3615    for (iRow=0;iRow<4;iRow++) {
     3616      if (rowArray_[iRow]) {
     3617        rowArray_[iRow]->clear();
     3618        rowArray_[iRow]->checkClear();
     3619      }
     3620      rowArray_[iRow]=NULL;
     3621    }
     3622    for (iColumn=0;iColumn<2;iColumn++) {
     3623      if (columnArray_[iColumn]) {
     3624        columnArray_[iColumn]->clear();
     3625        columnArray_[iColumn]->checkClear();
     3626      }
     3627      columnArray_[iColumn]=NULL;
     3628    }
     3629    rowCopy_=NULL;
     3630    ClpMatrixBase * temp = matrix_;
     3631    matrix_=auxiliaryModel_->matrix_;
     3632    auxiliaryModel_->matrix_=temp;
     3633    assert ((auxiliaryModel_->specialOptions_&16+32)==16+32);
     3634    if (problemStatus_!=1) {
     3635      int i;
     3636      if (auxiliaryModel_->rowScale_) {
     3637        const double * scale = auxiliaryModel_->columnScale_;
     3638        const double * inverseScale = scale + numberColumns_;
     3639        for (i=0;i<numberColumns_;i++) {
     3640          columnActivity_[i] = auxiliaryModel_->columnActivityWork_[i]*scale[i];
     3641          reducedCost_[i] = auxiliaryModel_->reducedCostWork_[i]*inverseScale[i];
     3642        }
     3643        scale = auxiliaryModel_->rowScale_;
     3644        inverseScale = scale + numberRows_;
     3645        for (i=0;i<numberRows_;i++) {
     3646          rowActivity_[i] = auxiliaryModel_->rowActivityWork_[i]*inverseScale[i];
     3647        }
    32923648      } else {
    3293         if (numberDualUnscaled)
    3294           secondaryStatus_=3;
    3295       }
    3296     }
    3297     if (problemStatus_==2) {
    3298       for (i=0;i<numberColumns_;i++) {
    3299         ray_[i] *= columnScale_[i];
    3300       }
    3301     } else if (problemStatus_==1&&ray_) {
    3302       for (i=0;i<numberRows_;i++) {
    3303         ray_[i] *= rowScale_[i];
    3304       }
    3305     }
    3306   } else if (rhsScale_!=1.0||objectiveScale_!=1.0) {
    3307     // Collect infeasibilities
    3308     int numberPrimalScaled=0;
    3309     int numberPrimalUnscaled=0;
    3310     int numberDualScaled=0;
    3311     int numberDualUnscaled=0;
    3312     double scaleC = 1.0/objectiveScale_;
    3313     double scaleR = 1.0/rhsScale_;
    3314     for (i=0;i<numberColumns_;i++) {
    3315       double valueScaled = columnActivityWork_[i];
    3316       double lowerScaled = columnLowerWork_[i];
    3317       double upperScaled = columnUpperWork_[i];
    3318       if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
    3319         if (valueScaled<lowerScaled-primalTolerance_||
    3320             valueScaled>upperScaled+primalTolerance_)
    3321           numberPrimalScaled++;
    3322         else
    3323           upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
    3324       }
    3325       columnActivity_[i] = valueScaled*scaleR;
    3326       double value = columnActivity_[i];
    3327       if (value<columnLower_[i]-primalTolerance_)
    3328         numberPrimalUnscaled++;
    3329       else if (value>columnUpper_[i]+primalTolerance_)
    3330         numberPrimalUnscaled++;
    3331       double valueScaledDual = reducedCostWork_[i];
    3332       if (valueScaled>columnLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
    3333         numberDualScaled++;
    3334       if (valueScaled<columnUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
    3335         numberDualScaled++;
    3336       reducedCost_[i] = valueScaledDual*scaleC;
    3337       double valueDual = reducedCost_[i];
    3338       if (value>columnLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
    3339         numberDualUnscaled++;
    3340       if (value<columnUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
    3341         numberDualUnscaled++;
    3342     }
    3343     for (i=0;i<numberRows_;i++) {
    3344       double valueScaled = rowActivityWork_[i];
    3345       double lowerScaled = rowLowerWork_[i];
    3346       double upperScaled = rowUpperWork_[i];
    3347       if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
    3348         if (valueScaled<lowerScaled-primalTolerance_||
    3349             valueScaled>upperScaled+primalTolerance_)
    3350           numberPrimalScaled++;
    3351         else
    3352           upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
    3353       }
    3354       rowActivity_[i] = valueScaled*scaleR;
    3355       double value = rowActivity_[i];
    3356       if (value<rowLower_[i]-primalTolerance_)
    3357         numberPrimalUnscaled++;
    3358       else if (value>rowUpper_[i]+primalTolerance_)
    3359         numberPrimalUnscaled++;
    3360       double valueScaledDual = dual_[i]+rowObjectiveWork_[i];;
    3361       if (valueScaled>rowLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
    3362         numberDualScaled++;
    3363       if (valueScaled<rowUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
    3364         numberDualScaled++;
    3365       dual_[i] *= scaleC;
    3366       double valueDual = dual_[i];
    3367       if (rowObjective_)
    3368         valueDual += rowObjective_[i];
    3369       if (value>rowLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
    3370         numberDualUnscaled++;
    3371       if (value<rowUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
    3372         numberDualUnscaled++;
    3373     }
    3374     if (!problemStatus_&&!secondaryStatus_) {
    3375       // See if we need to set secondary status
    3376       if (numberPrimalUnscaled) {
    3377         if (numberDualUnscaled)
    3378           secondaryStatus_=4;
    3379         else
    3380           secondaryStatus_=2;
    3381       } else {
    3382         if (numberDualUnscaled)
    3383           secondaryStatus_=3;
    3384       }
    3385     }
    3386   } else {
    3387     if (columnActivityWork_) {
    3388       for (i=0;i<numberColumns_;i++) {
    3389         double value = columnActivityWork_[i];
    3390         double lower = columnLowerWork_[i];
    3391         double upper = columnUpperWork_[i];
    3392         if (lower>-1.0e20||upper<1.0e20) {
    3393           if (value>lower&&value<upper)
    3394             upperOut_ = CoinMax(upperOut_,CoinMin(value-lower,upper-value));
    3395         }
    3396         columnActivity_[i] = columnActivityWork_[i];
    3397         reducedCost_[i] = reducedCostWork_[i];
    3398       }
    3399       for (i=0;i<numberRows_;i++) {
    3400         double value = rowActivityWork_[i];
    3401         double lower = rowLowerWork_[i];
    3402         double upper = rowUpperWork_[i];
    3403         if (lower>-1.0e20||upper<1.0e20) {
    3404           if (value>lower&&value<upper)
    3405             upperOut_ = CoinMax(upperOut_,CoinMin(value-lower,upper-value));
    3406         }
    3407         rowActivity_[i] = rowActivityWork_[i];
    3408       }
    3409     }
    3410   }
    3411   // switch off scalefactor if auto
    3412   if (automaticScale_) {
    3413     rhsScale_=1.0;
    3414     objectiveScale_=1.0;
    3415   }
    3416   if (optimizationDirection_!=1.0) {
    3417     // and modify all dual signs
    3418     for (i=0;i<numberColumns_;i++)
    3419       reducedCost_[i] *= optimizationDirection_;
    3420     for (i=0;i<numberRows_;i++)
    3421       dual_[i] *= optimizationDirection_;
     3649        for (i=0;i<numberColumns_;i++) {
     3650          columnActivity_[i] = auxiliaryModel_->columnActivityWork_[i];
     3651          reducedCost_[i] = auxiliaryModel_->reducedCostWork_[i];
     3652        }
     3653        for (i=0;i<numberRows_;i++) {
     3654          rowActivity_[i] = auxiliaryModel_->rowActivityWork_[i];
     3655        }
     3656      }
     3657      if (optimizationDirection_!=1.0) {
     3658        // and modify reduced costs
     3659        for (i=0;i<numberColumns_;i++)
     3660          reducedCost_[i] *= optimizationDirection_;
     3661      }
     3662    }
    34223663  }
    34233664  // scaling may have been turned off
     
    74497690        objectiveWork_[elementIndex] = direction*elementValue;
    74507691      } else {
     7692        assert (!auxiliaryModel_);
    74517693        objectiveWork_[elementIndex] = direction*elementValue
    74527694          *columnScale_[elementIndex];
     
    74777719        rowLowerWork_[elementIndex] = elementValue * rhsScale_;
    74787720      } else {
     7721        assert (!auxiliaryModel_);
    74797722        rowLowerWork_[elementIndex] = elementValue * rhsScale_
    74807723          * rowScale_[elementIndex];
     
    75067749        rowUpperWork_[elementIndex] = elementValue * rhsScale_;
    75077750      } else {
     7751        assert (!auxiliaryModel_);
    75087752        rowUpperWork_[elementIndex] = elementValue * rhsScale_
    75097753          * rowScale_[elementIndex];
     
    75387782        rowLowerWork_[elementIndex] = lowerValue * rhsScale_;
    75397783      } else {
     7784        assert (!auxiliaryModel_);
    75407785        rowLowerWork_[elementIndex] = lowerValue * rhsScale_
    75417786          * rowScale_[elementIndex];
     
    75537798        rowUpperWork_[elementIndex] = upperValue * rhsScale_;
    75547799      } else {
     7800        assert (!auxiliaryModel_);
    75557801        rowUpperWork_[elementIndex] = upperValue * rhsScale_
    75567802          * rowScale_[elementIndex];
     
    76027848        rowLowerWork_[iRow] = rowLower_[iRow] * rhsScale_;
    76037849      } else {
     7850        assert (!auxiliaryModel_);
    76047851        rowLowerWork_[iRow] = rowLower_[iRow] * rhsScale_
    76057852          * rowScale_[iRow];
     
    76107857        rowUpperWork_[iRow] = rowUpper_[iRow] * rhsScale_;
    76117858      } else {
     7859        assert (!auxiliaryModel_);
    76127860        rowUpperWork_[iRow] = rowUpper_[iRow] * rhsScale_
    76137861          * rowScale_[iRow];
     
    76407888        columnLowerWork_[elementIndex] = elementValue * rhsScale_;
    76417889      } else {
     7890        assert (!auxiliaryModel_);
    76427891        columnLowerWork_[elementIndex] = elementValue * rhsScale_
    76437892          / columnScale_[elementIndex];
     
    76707919        columnUpperWork_[elementIndex] = elementValue * rhsScale_;
    76717920      } else {
     7921        assert (!auxiliaryModel_);
    76727922        columnUpperWork_[elementIndex] = elementValue * rhsScale_
    76737923          / columnScale_[elementIndex];
     
    77007950        columnLowerWork_[elementIndex] = lowerValue * rhsScale_;
    77017951      } else {
     7952        assert (!auxiliaryModel_);
    77027953        columnLowerWork_[elementIndex] = lowerValue * rhsScale_
    77037954          / columnScale_[elementIndex];
     
    77187969        columnUpperWork_[elementIndex] = upperValue * rhsScale_;
    77197970      } else {
     7971        assert (!auxiliaryModel_);
    77207972        columnUpperWork_[elementIndex] = upperValue * rhsScale_
    77217973          / columnScale_[elementIndex];
     
    77678019        columnLowerWork_[iColumn] = columnLower_[iColumn] * rhsScale_;
    77688020      } else {
     8021        assert (!auxiliaryModel_);
    77698022        columnLowerWork_[iColumn] = columnLower_[iColumn] * rhsScale_
    77708023          / columnScale_[iColumn];
     
    77758028        columnUpperWork_[iColumn] = columnUpper_[iColumn] * rhsScale_;
    77768029      } else {
     8030        assert (!auxiliaryModel_);
    77778031        columnUpperWork_[iColumn] = columnUpper_[iColumn] * rhsScale_
    77788032          / columnScale_[iColumn];
     
    78718125}
    78728126/*
     8127  If you are re-using the same matrix again and again then the setup time
     8128  to do scaling may be significant.  Also you may not want to initialize all values
     8129  or return all values (especially if infeasible).  While an auxiliary model exists
     8130  it will be faster.  If options -1 then model is switched off.  Otherwise switched on
     8131  with following options.
     8132  1 - rhs is constant
     8133  2 - bounds are constant
     8134  4 - objective is constant
     8135  8 - solution in by basis and no djs etc in
     8136  16 - no duals out (but reduced costs)
     8137  32 - no output if infeasible
     8138*/
     8139void
     8140ClpSimplex::auxiliaryModel(int options)
     8141{
     8142  delete auxiliaryModel_;
     8143  auxiliaryModel_=NULL;
     8144  if (options>=0) {
     8145    createRim(63,true,0);
     8146    auxiliaryModel_ = new ClpSimplex();
     8147    auxiliaryModel_->specialOptions_=options;
     8148    int i;
     8149    int numberRows2 = numberRows_+numberExtraRows_;
     8150    int numberTotal = numberRows2+numberColumns_;
     8151    if (rowScale_) {
     8152      auxiliaryModel_->rowScale_= new double [2*numberRows_];
     8153      for (i=0;i<numberRows_;i++) {
     8154        auxiliaryModel_->rowScale_[i]=rowScale_[i];
     8155        auxiliaryModel_->rowScale_[numberRows_+i]=1.0/rowScale_[i];
     8156      }
     8157      auxiliaryModel_->columnScale_= new double [2*numberColumns_];
     8158      for (i=0;i<numberColumns_;i++) {
     8159        auxiliaryModel_->columnScale_[i]=columnScale_[i];
     8160        auxiliaryModel_->columnScale_[numberColumns_+i]=1.0/columnScale_[i];
     8161      }
     8162    }
     8163    // copy anyway
     8164    auxiliaryModel_->lower_ = new double [2*numberTotal];
     8165    memcpy(auxiliaryModel_->lower_,lower_,numberTotal*sizeof(double));
     8166    memcpy(auxiliaryModel_->lower_+numberTotal,lower_,numberTotal*sizeof(double));
     8167    auxiliaryModel_->upper_ = new double [2*numberTotal];
     8168    memcpy(auxiliaryModel_->upper_,upper_,numberTotal*sizeof(double));
     8169    memcpy(auxiliaryModel_->upper_+numberTotal,upper_,numberTotal*sizeof(double));
     8170    auxiliaryModel_->cost_ = new double [2*numberTotal];
     8171    memcpy(auxiliaryModel_->cost_,cost_,numberTotal*sizeof(double));
     8172    memcpy(auxiliaryModel_->cost_+numberTotal,cost_,numberTotal*sizeof(double));
     8173    auxiliaryModel_->dj_ = new double [2*numberTotal];
     8174    memset(auxiliaryModel_->dj_,0,2*numberTotal*sizeof(double));
     8175    auxiliaryModel_->solution_ = new double [2*numberTotal];
     8176    memset(auxiliaryModel_->solution_,0,2*numberTotal*sizeof(double));
     8177    auxiliaryModel_->reducedCostWork_ = auxiliaryModel_->dj_;
     8178    auxiliaryModel_->rowReducedCost_ = auxiliaryModel_->dj_+numberColumns_;
     8179    auxiliaryModel_->columnActivityWork_ = auxiliaryModel_->solution_;
     8180    auxiliaryModel_->rowActivityWork_ = auxiliaryModel_->solution_+numberColumns_;
     8181    auxiliaryModel_->objectiveWork_ = auxiliaryModel_->cost_;
     8182    auxiliaryModel_->rowObjectiveWork_ = auxiliaryModel_->cost_+numberColumns_;
     8183    auxiliaryModel_->rowLowerWork_ = auxiliaryModel_->lower_+numberColumns_;
     8184    auxiliaryModel_->columnLowerWork_ = auxiliaryModel_->lower_;
     8185    auxiliaryModel_->rowUpperWork_ = auxiliaryModel_->upper_+numberColumns_;
     8186    auxiliaryModel_->columnUpperWork_ = auxiliaryModel_->upper_;
     8187    // delete in model
     8188    delete [] lower_;
     8189    lower_=NULL;
     8190    delete [] upper_;
     8191    upper_=NULL;
     8192    delete [] cost_;
     8193    cost_=NULL;
     8194    delete [] dj_;
     8195    dj_=NULL;
     8196    delete [] solution_;
     8197    solution_=NULL;
     8198    auxiliaryModel_->rowCopy_=rowCopy_;
     8199    ClpMatrixBase * columnCopy = rowCopy_->reverseOrderedCopy();
     8200    rowCopy_=NULL;
     8201    // Point to correct as will be restored
     8202    auxiliaryModel_->matrix_ = matrix_;
     8203    int iRow,iColumn;
     8204    for (iRow=0;iRow<4;iRow++) {
     8205      auxiliaryModel_->rowArray_[iRow]=rowArray_[iRow];
     8206      rowArray_[iRow]=NULL;
     8207    }
     8208    for (iColumn=0;iColumn<2;iColumn++) {
     8209      auxiliaryModel_->columnArray_[iColumn]=columnArray_[iColumn];
     8210      columnArray_[iColumn]=NULL;
     8211    }
     8212    deleteRim(1);
     8213    auxiliaryModel_->matrix_ = columnCopy;
     8214    // Get rid of scaling
     8215    delete [] rowScale_;
     8216    delete [] columnScale_;
     8217    rowScale_=NULL;
     8218    columnScale_=NULL;
     8219  }
     8220}
     8221/*
    78738222  When scaling is on it is possible that the scaled problem
    78748223  is feasible but the unscaled is not.  Clp returns a secondary
  • trunk/ClpSimplexDual.cpp

    r634 r636  
    124124#endif
    125125// dual
    126 int ClpSimplexDual::dual (int ifValuesPass , int startFinishOptions)
    127 {
    128126
    129127  /* *** Method
     
    212210
    213211  */
     212#define CLEAN_FIXED 0
     213// Startup part of dual (may be extended to other algorithms)
     214int
     215ClpSimplexDual::startupSolve(int ifValuesPass,double * saveDuals,int startFinishOptions)
     216{
     217  // If values pass then save given duals round check solution
     218  // sanity check
     219  // initialize - no values pass and algorithm_ is -1
     220  // put in standard form (and make row copy)
     221  // create modifiable copies of model rim and do optional scaling
     222  // If problem looks okay
     223  // Do initial factorization
     224  // If user asked for perturbation - do it
     225  if (!startup(0,startFinishOptions)) {
     226    // looks okay
     227    // Superbasic variables not allowed
     228    // If values pass then scale pi
     229    if (ifValuesPass) {
     230      if (problemStatus_&&perturbation_<100)
     231        perturb();
     232      int i;
     233      if (scalingFlag_>0) {
     234        for (i=0;i<numberRows_;i++) {
     235          dual_[i] = saveDuals[i]/rowScale_[i];
     236        }
     237      } else {
     238        memcpy(dual_,saveDuals,numberRows_*sizeof(double));
     239      }
     240      // now create my duals
     241      for (i=0;i<numberRows_;i++) {
     242        // slack
     243        double value = dual_[i];
     244        value += rowObjectiveWork_[i];
     245        saveDuals[i+numberColumns_]=value;
     246      }
     247      ClpDisjointCopyN(objectiveWork_,numberColumns_,saveDuals);
     248      transposeTimes(-1.0,dual_,saveDuals);
     249      // make reduced costs okay
     250      for (i=0;i<numberColumns_;i++) {
     251        if (getStatus(i)==atLowerBound) {
     252          if (saveDuals[i]<0.0) {
     253            //if (saveDuals[i]<-1.0e-3)
     254            //printf("bad dj at lb %d %g\n",i,saveDuals[i]);
     255            saveDuals[i]=0.0;
     256          }
     257        } else if (getStatus(i)==atUpperBound) {
     258          if (saveDuals[i]>0.0) {
     259            //if (saveDuals[i]>1.0e-3)
     260            //printf("bad dj at ub %d %g\n",i,saveDuals[i]);
     261            saveDuals[i]=0.0;
     262          }
     263        }
     264      }
     265      memcpy(dj_,saveDuals,(numberColumns_+numberRows_)*sizeof(double));
     266      // set up possible ones
     267      for (i=0;i<numberRows_+numberColumns_;i++)
     268        clearPivoted(i);
     269      int iRow;
     270      for (iRow=0;iRow<numberRows_;iRow++) {
     271        int iPivot=pivotVariable_[iRow];
     272        if (fabs(saveDuals[iPivot])>dualTolerance_) {
     273          if (getStatus(iPivot)!=isFree)
     274            setPivoted(iPivot);
     275        }
     276      }
     277    } else if ((specialOptions_&1024)!=0&&CLEAN_FIXED) {
     278      // set up possible ones
     279      for (int i=0;i<numberRows_+numberColumns_;i++)
     280        clearPivoted(i);
     281      int iRow;
     282      for (iRow=0;iRow<numberRows_;iRow++) {
     283        int iPivot=pivotVariable_[iRow];
     284        if (iPivot<numberColumns_&&lower_[iPivot]==upper_[iPivot]) {
     285          setPivoted(iPivot);
     286        }
     287      }
     288    }
     289
     290    double objectiveChange;
     291    numberFake_ =0; // Number of variables at fake bounds
     292    numberChanged_ =0; // Number of variables with changed costs
     293    changeBounds(true,NULL,objectiveChange);
     294   
     295    if (!ifValuesPass) {
     296      // Check optimal
     297      if (!numberDualInfeasibilities_&&!numberPrimalInfeasibilities_)
     298        problemStatus_=0;
     299    }
     300    if (problemStatus_<0&&perturbation_<100) {
     301      perturb();
     302      // Can't get here if values pass
     303      gutsOfSolution(NULL,NULL);
     304      if (handler_->logLevel()>2) {
     305        handler_->message(CLP_SIMPLEX_STATUS,messages_)
     306          <<numberIterations_<<objectiveValue();
     307        handler_->printing(sumPrimalInfeasibilities_>0.0)
     308          <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
     309        handler_->printing(sumDualInfeasibilities_>0.0)
     310          <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
     311        handler_->printing(numberDualInfeasibilitiesWithoutFree_
     312                           <numberDualInfeasibilities_)
     313                             <<numberDualInfeasibilitiesWithoutFree_;
     314        handler_->message()<<CoinMessageEol;
     315      }
     316    }
     317    return 0;
     318  } else {
     319    return 1;
     320  }
     321}
     322void
     323ClpSimplexDual::finishSolve(int startFinishOptions)
     324{
     325  assert (problemStatus_||!sumPrimalInfeasibilities_);
     326
     327  // clean up
     328  finish(startFinishOptions);
     329}
     330void
     331ClpSimplexDual::gutsOfDual(int ifValuesPass,double * & saveDuals,int initialStatus,
     332                           ClpDataSave & data)
     333{
     334  int lastCleaned=0; // last time objective or bounds cleaned up
     335 
     336  // This says whether to restore things etc
     337  // startup will have factorized so can skip
     338  int factorType=0;
     339  // Start check for cycles
     340  progress_->startCheck();
     341  // Say change made on first iteration
     342  changeMade_=1;
     343  /*
     344    Status of problem:
     345    0 - optimal
     346    1 - infeasible
     347    2 - unbounded
     348    -1 - iterating
     349    -2 - factorization wanted
     350    -3 - redo checking without factorization
     351    -4 - looks infeasible
     352  */
     353  while (problemStatus_<0) {
     354    int iRow, iColumn;
     355    // clear
     356    for (iRow=0;iRow<4;iRow++) {
     357      rowArray_[iRow]->clear();
     358    }   
     359   
     360    for (iColumn=0;iColumn<2;iColumn++) {
     361      columnArray_[iColumn]->clear();
     362    }   
     363   
     364    // give matrix (and model costs and bounds a chance to be
     365    // refreshed (normally null)
     366    matrix_->refresh(this);
     367    // If getting nowhere - why not give it a kick
     368    // does not seem to work too well - do some more work
     369    if (perturbation_<101&&numberIterations_>2*(numberRows_+numberColumns_)
     370        &&initialStatus!=10) {
     371      perturb();
     372      // Can't get here if values pass
     373      gutsOfSolution(NULL,NULL);
     374      if (handler_->logLevel()>2) {
     375        handler_->message(CLP_SIMPLEX_STATUS,messages_)
     376          <<numberIterations_<<objectiveValue();
     377        handler_->printing(sumPrimalInfeasibilities_>0.0)
     378          <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
     379        handler_->printing(sumDualInfeasibilities_>0.0)
     380          <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
     381        handler_->printing(numberDualInfeasibilitiesWithoutFree_
     382                           <numberDualInfeasibilities_)
     383                             <<numberDualInfeasibilitiesWithoutFree_;
     384        handler_->message()<<CoinMessageEol;
     385      }
     386    }
     387    // may factorize, checks if problem finished
     388    statusOfProblemInDual(lastCleaned,factorType,saveDuals,data,
     389                          ifValuesPass);
     390    // If values pass then do easy ones on first time
     391    if (ifValuesPass&&
     392        progress_->lastIterationNumber(0)<0&&saveDuals) {
     393      doEasyOnesInValuesPass(saveDuals);
     394    }
     395   
     396    // Say good factorization
     397    factorType=1;
     398    if (data.sparseThreshold_) {
     399      // use default at present
     400      factorization_->sparseThreshold(0);
     401      factorization_->goSparse();
     402    }
     403   
     404    // exit if victory declared
     405    if (problemStatus_>=0)
     406      break;
     407   
     408    // test for maximum iterations
     409    if (hitMaximumIterations()||(ifValuesPass==2&&!saveDuals)) {
     410      problemStatus_=3;
     411      break;
     412    }
     413    if (ifValuesPass&&!saveDuals) {
     414      // end of values pass
     415      ifValuesPass=0;
     416      int status = eventHandler_->event(ClpEventHandler::endOfValuesPass);
     417      if (status>=0) {
     418        problemStatus_=5;
     419        secondaryStatus_=ClpEventHandler::endOfValuesPass;
     420        break;
     421      }
     422    }
     423    // Check event
     424    {
     425      int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
     426      if (status>=0) {
     427        problemStatus_=5;
     428        secondaryStatus_=ClpEventHandler::endOfFactorization;
     429        break;
     430      }
     431    }
     432      // Do iterations
     433    whileIterating(saveDuals,ifValuesPass);
     434  }
     435}
     436int
     437ClpSimplexDual::dual(int ifValuesPass,int startFinishOptions)
     438{
    214439  algorithm_ = -1;
    215 
     440 
    216441  // save data
    217442  ClpDataSave data = saveData();
    218 #define CLEAN_FIXED 0
     443  double * saveDuals = NULL;
     444  if (ifValuesPass) {
     445    saveDuals = new double [numberRows_+numberColumns_];
     446    memcpy(saveDuals,dual_,numberRows_*sizeof(double));
     447  }
     448  int returnCode = startupSolve(ifValuesPass,saveDuals,startFinishOptions);
     449  // Save so can see if doing after primal
     450  int initialStatus=problemStatus_;
     451 
     452  if (!returnCode)
     453    gutsOfDual(ifValuesPass,saveDuals,initialStatus,data);
     454  finishSolve(startFinishOptions);
     455  delete [] saveDuals;
     456 
     457  // Restore any saved stuff
     458  restoreData(data);
     459  return problemStatus_;
     460}
     461// old way
     462#if 0
     463int ClpSimplexDual::dual (int ifValuesPass , int startFinishOptions)
     464{
     465  algorithm_ = -1;
     466
     467  // save data
     468  ClpDataSave data = saveData();
    219469  // Save so can see if doing after primal
    220470  int initialStatus=problemStatus_;
     
    439689  return problemStatus_;
    440690}
     691#endif
    441692//#define CHECK_ACCURACY
    442693#ifdef CHECK_ACCURACY
     
    36973948  if (getFakeBound(iSequence)!=noFake)
    36983949    numberFake_--;;
     3950  if (auxiliaryModel_) {
     3951    // just copy back
     3952    lower_[iSequence]=auxiliaryModel_->lowerRegion()[iSequence+numberRows_+numberColumns_];
     3953    upper_[iSequence]=auxiliaryModel_->upperRegion()[iSequence+numberRows_+numberColumns_];
     3954    return;
     3955  }
    36993956  if (iSequence>=numberColumns_) {
    37003957    // rows
     
    44974754  // create modifiable copies of model rim and do optional scaling
    44984755  int startFinishOptions;
    4499   if((specialOptions_&4096)==0) {
     4756  if((specialOptions_&4096)==0||auxiliaryModel_) {
    45004757    startFinishOptions=0;
    45014758  } else {
     
    45674824{
    45684825  int startFinishOptions;
    4569   if((specialOptions_&4096)==0) {
     4826  if((specialOptions_&4096)==0||auxiliaryModel_) {
    45704827    startFinishOptions=0;
    45714828  } else {
     
    50945351      const int * thisIndices = column+rowStart[iRow];
    50955352      if (rowScale_) {
     5353        assert (!auxiliaryModel_);
    50965354        // scale row
    50975355        double scale = rowScale_[iRow];
  • trunk/ClpSimplexOther.cpp

    r631 r636  
    88#include "CoinHelperFunctions.hpp"
    99#include "ClpSimplexOther.hpp"
     10#include "ClpSimplexDual.hpp"
     11#include "ClpSimplexPrimal.hpp"
     12#include "ClpEventHandler.hpp"
     13#include "ClpHelperFunctions.hpp"
    1014#include "ClpFactorization.hpp"
     15#include "ClpDualRowDantzig.hpp"
    1116#include "CoinPackedMatrix.hpp"
    1217#include "CoinIndexedVector.hpp"
     
    3237                            double * costDecreased, int * sequenceDecreased)
    3338{
    34   rowArray_[0]->clear();
    3539  rowArray_[1]->clear();
    3640  columnArray_[1]->clear();
    37   columnArray_[0]->clear();
    3841  // long enough for rows+columns
    3942  assert(rowArray_[3]->capacity()>=numberRows_+numberColumns_);
     
    5457  dualTolerance_ = dblParam_[ClpDualTolerance];
    5558  for ( i=0;i<numberCheck;i++) {
     59    rowArray_[0]->clear();
     60    //rowArray_[0]->checkClear();
     61    //rowArray_[1]->checkClear();
     62    //columnArray_[1]->checkClear();
     63    columnArray_[0]->clear();
     64    //columnArray_[0]->checkClear();
    5665    int iSequence = which[i];
    5766    double costIncrease=COIN_DBL_MAX;
     
    8594              // we are going to use for cutoff so be exact
    8695              costIncrease = fabs(djValue/alphaIncrease);
    87               if(sequenceIncrease<numberColumns_&&integerType_[sequenceIncrease]) {
     96              /* Not sure this is good idea as I don't think correct e.g.
     97                 suppose a continuous variable has dj slightly greater. */
     98              if(false&&sequenceIncrease<numberColumns_&&integerType_[sequenceIncrease]) {
    8899                // can improve
    89100                double movement = (columnScale_==NULL) ? 1.0 :
     
    132143    }
    133144    double scaleFactor;
    134     if (rowScale_) {
    135       if (iSequence<numberColumns_)
    136         scaleFactor = 1.0/(objectiveScale_*columnScale_[iSequence]);
    137       else
    138         scaleFactor = rowScale_[iSequence-numberColumns_]/objectiveScale_;
     145    if (!auxiliaryModel_) {
     146      if (rowScale_) {
     147        if (iSequence<numberColumns_)
     148          scaleFactor = 1.0/(objectiveScale_*columnScale_[iSequence]);
     149        else
     150          scaleFactor = rowScale_[iSequence-numberColumns_]/objectiveScale_;
     151      } else {
     152        scaleFactor = 1.0/objectiveScale_;
     153      }
    139154    } else {
    140       scaleFactor = 1.0/objectiveScale_;
     155      if (auxiliaryModel_->rowScale()) {
     156        if (iSequence<numberColumns_)
     157          scaleFactor = 1.0/(objectiveScale_*auxiliaryModel_->columnScale()[iSequence]);
     158        else
     159          scaleFactor = auxiliaryModel_->rowScale()[iSequence-numberColumns_]/objectiveScale_;
     160      } else {
     161        scaleFactor = 1.0/objectiveScale_;
     162      }
    141163    }
    142164    if (costIncrease<1.0e30)
     
    164186    }
    165187  }
     188  //rowArray_[0]->clear();
     189  //rowArray_[1]->clear();
     190  //columnArray_[1]->clear();
     191  //columnArray_[0]->clear();
     192  //rowArray_[3]->clear();
    166193  if (!optimizationDirection_)
    167194    printf("*** ????? Ranging with zero optimization costs\n");
     
    944971ClpSimplex *
    945972ClpSimplexOther::crunch(double * rhs, int * whichRow, int * whichColumn,
    946                         int & nBound, bool moreBounds)
     973                        int & nBound, bool moreBounds, bool tightenBounds)
    947974{
    948975#if 0
     
    11381165      delete small;
    11391166      small=NULL;
     1167    } else if (tightenBounds) {
     1168      // See if we can tighten any bounds
     1169      // use rhs for upper and small duals for lower
     1170      double * up = rhs;
     1171      double * lo = small->dualRowSolution();
     1172      const double * element = small->clpMatrix()->getElements();
     1173      const int * row = small->clpMatrix()->getIndices();
     1174      const CoinBigIndex * columnStart = small->clpMatrix()->getVectorStarts();
     1175      //const int * columnLength = small->clpMatrix()->getVectorLengths();
     1176      CoinZeroN(lo,numberRows2);
     1177      CoinZeroN(up,numberRows2);
     1178      for (int iColumn=0;iColumn<numberColumns2;iColumn++) {
     1179        double upper=columnUpper2[iColumn];
     1180        double lower=columnLower2[iColumn];
     1181        //assert (columnLength[iColumn]==columnStart[iColumn+1]-columnStart[iColumn]);
     1182        for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn+1];j++) {
     1183          int iRow=row[j];
     1184          double value = element[j];
     1185          if (value>0.0) {
     1186            if (upper<1.0e20)
     1187              up[iRow] += upper*value;
     1188            else
     1189              up[iRow] = COIN_DBL_MAX;
     1190            if (lower>-1.0e20)
     1191              lo[iRow] += lower*value;
     1192            else
     1193              lo[iRow] = -COIN_DBL_MAX;
     1194          } else {
     1195            if (upper<1.0e20)
     1196              lo[iRow] += upper*value;
     1197            else
     1198              lo[iRow] = -COIN_DBL_MAX;
     1199            if (lower>-1.0e20)
     1200              up[iRow] += lower*value;
     1201            else
     1202              up[iRow] = COIN_DBL_MAX;
     1203          }
     1204        }
     1205      }
     1206      double * rowLower2 = small->rowLower();
     1207      double * rowUpper2 = small->rowUpper();
     1208      bool feasible=true;
     1209      // make safer
     1210      for (int iRow=0;iRow<numberRows2;iRow++) {
     1211        double lower = lo[iRow];
     1212        if (lower>rowUpper2[iRow]+tolerance) {
     1213          feasible=false;
     1214          break;
     1215        } else {
     1216          lo[iRow] = CoinMin(lower-rowUpper2[iRow],0.0)-tolerance;
     1217        }
     1218        double upper = up[iRow];
     1219        if (upper<rowLower2[iRow]-tolerance) {
     1220          feasible=false;
     1221          break;
     1222        } else {
     1223          up[iRow] = CoinMax(upper-rowLower2[iRow],0.0)+tolerance;
     1224        }
     1225      }
     1226      if (!feasible) {
     1227        delete small;
     1228        small=NULL;
     1229      } else {
     1230        // and tighten
     1231        for (int iColumn=0;iColumn<numberColumns2;iColumn++) {
     1232          if (integerInformation[whichColumn[iColumn]]) {
     1233            double upper=columnUpper2[iColumn];
     1234            double lower=columnLower2[iColumn];
     1235            double newUpper = upper;
     1236            double newLower = lower;
     1237            double difference = upper-lower;
     1238            if (lower>-1000.0&&upper<1000.0) {
     1239              for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn+1];j++) {
     1240                int iRow=row[j];
     1241                double value = element[j];
     1242                if (value>0.0) {
     1243                  double upWithOut = up[iRow] - value*difference;
     1244                  if (upWithOut<0.0) {
     1245                    newLower = CoinMax(newLower,lower-(upWithOut+tolerance)/value);
     1246                  }
     1247                  double lowWithOut = lo[iRow] + value*difference;
     1248                  if (lowWithOut>0.0) {
     1249                    newUpper = CoinMin(newUpper,upper-(lowWithOut-tolerance)/value);
     1250                  }
     1251                } else {
     1252                  double upWithOut = up[iRow] + value*difference;
     1253                  if (upWithOut<0.0) {
     1254                    newUpper = CoinMin(newUpper,upper-(upWithOut+tolerance)/value);
     1255                  }
     1256                  double lowWithOut = lo[iRow] - value*difference;
     1257                  if (lowWithOut>0.0) {
     1258                    newLower = CoinMax(newLower,lower-(lowWithOut-tolerance)/value);
     1259                  }
     1260                }
     1261              }
     1262              if (newLower>lower||newUpper<upper) {
     1263                if (fabs(newUpper-floor(newUpper+0.5))>1.0e-6)
     1264                  newUpper = floor(newUpper);
     1265                else
     1266                  newUpper = floor(newUpper+0.5);
     1267                if (fabs(newLower-ceil(newLower-0.5))>1.0e-6)
     1268                  newLower = ceil(newLower);
     1269                else
     1270                  newLower = ceil(newLower-0.5);
     1271                // change may be too small - check
     1272                if (newLower>lower||newUpper<upper) {
     1273                  if (newUpper>=newLower) {
     1274                    // Could also tighten in this
     1275                    //printf("%d bounds %g %g tightened to %g %g\n",
     1276                    //     iColumn,columnLower2[iColumn],columnUpper2[iColumn],
     1277                    //     newLower,newUpper);
     1278#if 1
     1279                    columnUpper2[iColumn]=newUpper;
     1280                    columnLower2[iColumn]=newLower;
     1281                    columnUpper_[whichColumn[iColumn]]=newUpper;
     1282                    columnLower_[whichColumn[iColumn]]=newLower;
     1283#endif
     1284                    // and adjust bounds on rows
     1285                    newUpper -= upper;
     1286                    newLower -= lower;
     1287                    for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn+1];j++) {
     1288                      int iRow=row[j];
     1289                      double value = element[j];
     1290                      if (value>0.0) {
     1291                        up[iRow] += newUpper*value;
     1292                        lo[iRow] += newLower*value;
     1293                      } else {
     1294                        lo[iRow] += newUpper*value;
     1295                        up[iRow] += newLower*value;
     1296                      }
     1297                    }
     1298                  } else {
     1299                    // infeasible
     1300                    //printf("%d bounds infeasible %g %g tightened to %g %g\n",
     1301                    //     iColumn,columnLower2[iColumn],columnUpper2[iColumn],
     1302                    //     newLower,newUpper);
     1303#if 1
     1304                    delete small;
     1305                    small=NULL;
     1306                    break;
     1307#endif
     1308                  }
     1309                }
     1310              }
     1311            }
     1312          }
     1313        }
     1314      }
    11401315    }
    11411316  }
     
    12201395#endif
    12211396}
     1397/* Tightens integer bounds - returns number tightened or -1 if infeasible
     1398 */
     1399int
     1400ClpSimplexOther::tightenIntegerBounds(double * rhsSpace)
     1401{
     1402  // See if we can tighten any bounds
     1403  // use rhs for upper and small duals for lower
     1404  double * up = rhsSpace;
     1405  double * lo = dual_;
     1406  const double * element = matrix_->getElements();
     1407  const int * row = matrix_->getIndices();
     1408  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
     1409  const int * columnLength = matrix_->getVectorLengths();
     1410  CoinZeroN(lo,numberRows_);
     1411  CoinZeroN(up,numberRows_);
     1412  for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
     1413    double upper=columnUpper_[iColumn];
     1414    double lower=columnLower_[iColumn];
     1415    //assert (columnLength[iColumn]==columnStart[iColumn+1]-columnStart[iColumn]);
     1416    for (CoinBigIndex j=columnStart[iColumn];
     1417         j<columnStart[iColumn]+columnLength[iColumn];j++) {
     1418      int iRow=row[j];
     1419      double value = element[j];
     1420      if (value>0.0) {
     1421        if (upper<1.0e20)
     1422          up[iRow] += upper*value;
     1423        else
     1424          up[iRow] = COIN_DBL_MAX;
     1425        if (lower>-1.0e20)
     1426          lo[iRow] += lower*value;
     1427        else
     1428          lo[iRow] = -COIN_DBL_MAX;
     1429      } else {
     1430        if (upper<1.0e20)
     1431          lo[iRow] += upper*value;
     1432        else
     1433          lo[iRow] = -COIN_DBL_MAX;
     1434        if (lower>-1.0e20)
     1435          up[iRow] += lower*value;
     1436        else
     1437          up[iRow] = COIN_DBL_MAX;
     1438      }
     1439    }
     1440  }
     1441  bool feasible=true;
     1442  // make safer
     1443  double tolerance = primalTolerance();
     1444  for (int iRow=0;iRow<numberRows_;iRow++) {
     1445    double lower = lo[iRow];
     1446    if (lower>rowUpper_[iRow]+tolerance) {
     1447      feasible=false;
     1448      break;
     1449    } else {
     1450      lo[iRow] = CoinMin(lower-rowUpper_[iRow],0.0)-tolerance;
     1451    }
     1452    double upper = up[iRow];
     1453    if (upper<rowLower_[iRow]-tolerance) {
     1454      feasible=false;
     1455      break;
     1456    } else {
     1457      up[iRow] = CoinMax(upper-rowLower_[iRow],0.0)+tolerance;
     1458    }
     1459  }
     1460  int numberTightened=0;
     1461  if (!feasible) {
     1462    return -1;
     1463  } else {
     1464    // and tighten
     1465    for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
     1466      if (integerType_[iColumn]) {
     1467        double upper=columnUpper_[iColumn];
     1468        double lower=columnLower_[iColumn];
     1469        double newUpper = upper;
     1470        double newLower = lower;
     1471        double difference = upper-lower;
     1472        if (lower>-1000.0&&upper<1000.0) {
     1473          for (CoinBigIndex j=columnStart[iColumn];
     1474               j<columnStart[iColumn]+columnLength[iColumn];j++) {
     1475            int iRow=row[j];
     1476            double value = element[j];
     1477            if (value>0.0) {
     1478              double upWithOut = up[iRow] - value*difference;
     1479              if (upWithOut<0.0) {
     1480                newLower = CoinMax(newLower,lower-(upWithOut+tolerance)/value);
     1481              }
     1482              double lowWithOut = lo[iRow] + value*difference;
     1483              if (lowWithOut>0.0) {
     1484                newUpper = CoinMin(newUpper,upper-(lowWithOut-tolerance)/value);
     1485              }
     1486            } else {
     1487              double upWithOut = up[iRow] + value*difference;
     1488              if (upWithOut<0.0) {
     1489                newUpper = CoinMin(newUpper,upper-(upWithOut+tolerance)/value);
     1490              }
     1491              double lowWithOut = lo[iRow] - value*difference;
     1492              if (lowWithOut>0.0) {
     1493                newLower = CoinMax(newLower,lower-(lowWithOut-tolerance)/value);
     1494              }
     1495            }
     1496          }
     1497          if (newLower>lower||newUpper<upper) {
     1498            if (fabs(newUpper-floor(newUpper+0.5))>1.0e-6)
     1499              newUpper = floor(newUpper);
     1500            else
     1501              newUpper = floor(newUpper+0.5);
     1502            if (fabs(newLower-ceil(newLower-0.5))>1.0e-6)
     1503              newLower = ceil(newLower);
     1504            else
     1505              newLower = ceil(newLower-0.5);
     1506            // change may be too small - check
     1507            if (newLower>lower||newUpper<upper) {
     1508              if (newUpper>=newLower) {
     1509                numberTightened++;
     1510                //printf("%d bounds %g %g tightened to %g %g\n",
     1511                //     iColumn,columnLower_[iColumn],columnUpper_[iColumn],
     1512                //     newLower,newUpper);
     1513                columnUpper_[iColumn]=newUpper;
     1514                columnLower_[iColumn]=newLower;
     1515                // and adjust bounds on rows
     1516                newUpper -= upper;
     1517                newLower -= lower;
     1518                for (CoinBigIndex j=columnStart[iColumn];
     1519                     j<columnStart[iColumn]+columnLength[iColumn];j++) {
     1520                  int iRow=row[j];
     1521                  double value = element[j];
     1522                  if (value>0.0) {
     1523                    up[iRow] += newUpper*value;
     1524                    lo[iRow] += newLower*value;
     1525                  } else {
     1526                    lo[iRow] += newUpper*value;
     1527                    up[iRow] += newLower*value;
     1528                  }
     1529                }
     1530              } else {
     1531                // infeasible
     1532                //printf("%d bounds infeasible %g %g tightened to %g %g\n",
     1533                //     iColumn,columnLower_[iColumn],columnUpper_[iColumn],
     1534                //     newLower,newUpper);
     1535                return -1;
     1536              }
     1537            }
     1538          }
     1539        }
     1540      }
     1541    }
     1542  }
     1543  return numberTightened;
     1544}
     1545/* Parametrics
     1546   This is an initial slow version.
     1547   The code uses current bounds + theta * change (if change array not NULL)
     1548   and similarly for objective.
     1549   It starts at startingTheta and returns ending theta in endingTheta.
     1550   If reportIncrement 0.0 it will report on any movement
     1551   If reportIncrement >0.0 it will report at startingTheta+k*reportIncrement.
     1552   If it can not reach input endingTheta return code will be 1 for infeasible,
     1553   2 for unbounded, if error on ranges -1,  otherwise 0.
     1554   Normal report is just theta and objective but
     1555   if event handler exists it may do more
     1556   On exit endingTheta is maximum reached (can be used for next startingTheta)
     1557*/
     1558int
     1559ClpSimplexOther::parametrics(double startingTheta, double & endingTheta,double reportIncrement,
     1560                             const double * changeLowerBound, const double * changeUpperBound,
     1561                             const double * changeLowerRhs, const double * changeUpperRhs,
     1562                             const double * changeObjective)
     1563{
     1564  bool needToDoSomething=true;
     1565  bool canTryQuick = (reportIncrement) ? true : false;
     1566  // Save copy of model
     1567  ClpSimplex copyModel = *this;
     1568  int savePerturbation = perturbation_;
     1569  perturbation_=102; // switch off
     1570  while (needToDoSomething) {
     1571    needToDoSomething=false;
     1572    algorithm_ = -1;
     1573   
     1574    // save data
     1575    ClpDataSave data = saveData();
     1576    int returnCode = ((ClpSimplexDual *) this)->startupSolve(0,NULL,0);
     1577    int iRow,iColumn;
     1578    double * chgUpper=NULL;
     1579    double * chgLower=NULL;
     1580    double * chgObjective=NULL;
     1581   
     1582    // Dantzig (as will not be used) (out later)
     1583    ClpDualRowPivot * savePivot = dualRowPivot_;
     1584    dualRowPivot_ = new ClpDualRowDantzig();
     1585   
     1586    if (!returnCode) {
     1587      // Find theta when bounds will cross over and create arrays
     1588      int numberTotal = numberRows_+numberColumns_;
     1589      chgLower = new double[numberTotal];
     1590      memset(chgLower,0,numberTotal*sizeof(double));
     1591      chgUpper = new double[numberTotal];
     1592      memset(chgUpper,0,numberTotal*sizeof(double));
     1593      chgObjective = new double[numberTotal];
     1594      memset(chgObjective,0,numberTotal*sizeof(double));
     1595      assert (!rowScale_);
     1596      double maxTheta=1.0e50;
     1597      if (changeLowerRhs||changeUpperRhs) {
     1598        for (iRow=0;iRow<numberRows_;iRow++) {
     1599          double lower = rowLower_[iRow];
     1600          double upper = rowUpper_[iRow];
     1601          if (lower>upper) {
     1602            maxTheta=-1.0;
     1603            break;
     1604          }
     1605          double changeLower = (changeLowerRhs) ? changeLowerRhs[iRow] : 0.0;
     1606          double changeUpper = (changeUpperRhs) ? changeUpperRhs[iRow] : 0.0;
     1607          if (lower>-1.0e20&&upper<1.0e20) {
     1608            if (lower+maxTheta*changeLower>upper+maxTheta*changeUpper) {
     1609              maxTheta = (upper-lower)/(changeLower-changeUpper);
     1610            }
     1611          }
     1612          if (lower>-1.0e20) {
     1613            lower_[numberColumns_+iRow] += startingTheta*changeLower;
     1614            chgLower[numberColumns_+iRow]=changeLower;
     1615          }
     1616          if (upper<1.0e20) {
     1617            upper_[numberColumns_+iRow] += startingTheta*changeUpper;
     1618            chgUpper[numberColumns_+iRow]=changeUpper;
     1619          }
     1620        }
     1621      }
     1622      if (maxTheta>0.0) {
     1623        if (changeLowerBound||changeUpperBound) {
     1624          for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     1625            double lower = columnLower_[iColumn];
     1626            double upper = columnUpper_[iColumn];
     1627            if (lower>upper) {
     1628              maxTheta=-1.0;
     1629              break;
     1630            }
     1631            double changeLower = (changeLowerBound) ? changeLowerBound[iColumn] : 0.0;
     1632            double changeUpper = (changeUpperBound) ? changeUpperBound[iColumn] : 0.0;
     1633            if (lower>-1.0e20&&upper<1.0e20) {
     1634              if (lower+maxTheta*changeLower>upper+maxTheta*changeUpper) {
     1635                maxTheta = (upper-lower)/(changeLower-changeUpper);
     1636              }
     1637            }
     1638            if (lower>-1.0e20) {
     1639              lower_[iColumn] += startingTheta*changeLower;
     1640              chgLower[iColumn]=changeLower;
     1641            }
     1642            if (upper<1.0e20) {
     1643              upper_[iColumn] += startingTheta*changeUpper;
     1644              chgUpper[iColumn]=changeUpper;
     1645            }
     1646          }
     1647        }
     1648        if (maxTheta==1.0e50)
     1649          maxTheta = COIN_DBL_MAX;
     1650      }
     1651      if (maxTheta<0.0) {
     1652        // bad ranges or initial
     1653        returnCode = -1;
     1654      }
     1655      endingTheta = CoinMin(endingTheta,maxTheta);
     1656      if (endingTheta<startingTheta) {
     1657        // bad initial
     1658        returnCode = -2;
     1659      }
     1660    }
     1661    double saveEndingTheta=endingTheta;
     1662    if (!returnCode) {
     1663      if (changeObjective) {
     1664        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     1665          chgObjective[iColumn] = changeObjective[iColumn];
     1666          cost_[iColumn] += startingTheta*changeObjective[iColumn];
     1667        }
     1668      }
     1669      double * saveDuals=NULL;
     1670      ((ClpSimplexDual *) this)->gutsOfDual(0,saveDuals,-1,data);
     1671      assert (!problemStatus_);
     1672      // Now do parametrics
     1673      printf("at starting theta of %g objective value is %g\n",startingTheta,
     1674             objectiveValue());
     1675      while (!returnCode) {
     1676        assert (reportIncrement);
     1677        returnCode = parametricsLoop(startingTheta,endingTheta,reportIncrement,
     1678                                     chgLower,chgUpper,chgObjective,data,
     1679                                     canTryQuick);
     1680        if (!returnCode) {
     1681          //double change = endingTheta-startingTheta;
     1682          startingTheta=endingTheta;
     1683          endingTheta = saveEndingTheta;
     1684          //for (int i=0;i<numberTotal;i++) {
     1685          //lower_[i] += change*chgLower[i];
     1686          //upper_[i] += change*chgUpper[i];
     1687          //cost_[i] += change*chgObjective[i];
     1688          //}
     1689          printf("at theta of %g objective value is %g\n",startingTheta,
     1690                 objectiveValue());
     1691          if (startingTheta>=endingTheta)
     1692            break;
     1693        } else if (returnCode==-1) {
     1694          // trouble - do external solve
     1695          needToDoSomething=true;
     1696        } else {
     1697          abort();
     1698        }
     1699      }
     1700    }
     1701    ((ClpSimplexDual *) this)->finishSolve(0);
     1702   
     1703    delete dualRowPivot_;
     1704    dualRowPivot_ = savePivot;
     1705    // Restore any saved stuff
     1706    restoreData(data);
     1707    if (needToDoSomething) {
     1708      double saveStartingTheta=startingTheta; // known to be feasible
     1709      int cleanedUp=1;
     1710      while (cleanedUp) {
     1711        // tweak
     1712        if (cleanedUp==1) {
     1713          if (!reportIncrement)
     1714            startingTheta = CoinMin(startingTheta+1.0e-5,saveEndingTheta);
     1715          else
     1716            startingTheta = CoinMin(startingTheta+reportIncrement,saveEndingTheta);
     1717        } else {
     1718          // restoring to go slowly
     1719          startingTheta=saveStartingTheta;
     1720        }
     1721        // only works if not scaled
     1722        int i;
     1723        const double * obj1 = objective();
     1724        double * obj2 = copyModel.objective();
     1725        const double * lower1 = columnLower_;
     1726        double * lower2 = copyModel.columnLower();
     1727        const double * upper1 = columnUpper_;
     1728        double * upper2 = copyModel.columnUpper();
     1729        for (i=0;i<numberColumns_;i++) {
     1730          obj2[i] = obj1[i] + startingTheta*chgObjective[i];
     1731          lower2[i] = lower1[i] + startingTheta*chgLower[i];
     1732          upper2[i] = upper1[i] + startingTheta*chgUpper[i];
     1733        }
     1734        lower1 = rowLower_;
     1735        lower2 = copyModel.rowLower();
     1736        upper1 = rowUpper_;
     1737        upper2 = copyModel.rowUpper();
     1738        for (i=0;i<numberRows_;i++) {
     1739          lower2[i] = lower1[i] + startingTheta*chgLower[i+numberColumns_];
     1740          upper2[i] = upper1[i] + startingTheta*chgUpper[i+numberColumns_];
     1741        }
     1742        copyModel.dual();
     1743        if (copyModel.problemStatus()) {
     1744          printf("Can not get to theta of %g\n",startingTheta);
     1745          canTryQuick=false; // do slowly to get exact amount
     1746          // back to last known good
     1747          if (cleanedUp==1)
     1748            cleanedUp=2;
     1749          else
     1750            abort();
     1751        } else {
     1752          // and move stuff back
     1753          int numberTotal = numberRows_+numberColumns_;
     1754          memcpy(status_,copyModel.statusArray(),numberTotal);
     1755          memcpy(columnActivity_,copyModel.primalColumnSolution(),numberColumns_*sizeof(double));
     1756          memcpy(rowActivity_,copyModel.primalRowSolution(),numberRows_*sizeof(double));
     1757          cleanedUp=0;
     1758        }
     1759      }
     1760    }
     1761    delete [] chgLower;
     1762    delete [] chgUpper;
     1763    delete [] chgObjective;
     1764  }
     1765  perturbation_ = savePerturbation;
     1766  return problemStatus_;
     1767}
     1768int
     1769ClpSimplexOther::parametricsLoop(double startingTheta, double & endingTheta,double reportIncrement,
     1770                                 const double * changeLower, const double * changeUpper,
     1771                                 const double * changeObjective, ClpDataSave & data,
     1772                                 bool canTryQuick)
     1773{
     1774  // stuff is already at starting
     1775  // For this crude version just try and go to end
     1776  double change=0.0;
     1777  if (reportIncrement&&canTryQuick) {
     1778    endingTheta = CoinMin(endingTheta,startingTheta+reportIncrement);
     1779    change = endingTheta-startingTheta;
     1780  }
     1781  int numberTotal = numberRows_+numberColumns_;
     1782  int i;
     1783  for ( i=0;i<numberTotal;i++) {
     1784    lower_[i] += change*changeLower[i];
     1785    upper_[i] += change*changeUpper[i];
     1786    switch(getStatus(i)) {
     1787     
     1788    case basic:
     1789    case isFree:
     1790    case superBasic:
     1791      break;
     1792    case isFixed:
     1793    case atUpperBound:
     1794      solution_[i]=upper_[i];
     1795      break;
     1796    case atLowerBound:
     1797      solution_[i]=lower_[i];
     1798      break;
     1799    }
     1800    cost_[i] += change*changeObjective[i];
     1801  }
     1802  problemStatus_=-1;
     1803 
     1804  // This says whether to restore things etc
     1805  // startup will have factorized so can skip
     1806  int factorType=0;
     1807  // Start check for cycles
     1808  progress_->startCheck();
     1809  // Say change made on first iteration
     1810  changeMade_=1;
     1811  /*
     1812    Status of problem:
     1813    0 - optimal
     1814    1 - infeasible
     1815    2 - unbounded
     1816    -1 - iterating
     1817    -2 - factorization wanted
     1818    -3 - redo checking without factorization
     1819    -4 - looks infeasible
     1820  */
     1821  while (problemStatus_<0) {
     1822    int iRow, iColumn;
     1823    // clear
     1824    for (iRow=0;iRow<4;iRow++) {
     1825      rowArray_[iRow]->clear();
     1826    }   
     1827   
     1828    for (iColumn=0;iColumn<2;iColumn++) {
     1829      columnArray_[iColumn]->clear();
     1830    }   
     1831   
     1832    // give matrix (and model costs and bounds a chance to be
     1833    // refreshed (normally null)
     1834    matrix_->refresh(this);
     1835    // may factorize, checks if problem finished
     1836    statusOfProblemInParametrics(factorType,data);
     1837    // Say good factorization
     1838    factorType=1;
     1839    if (data.sparseThreshold_) {
     1840      // use default at present
     1841      factorization_->sparseThreshold(0);
     1842      factorization_->goSparse();
     1843    }
     1844   
     1845    // exit if victory declared
     1846    if (problemStatus_>=0)
     1847      break;
     1848   
     1849    // test for maximum iterations
     1850    if (hitMaximumIterations()) {
     1851      problemStatus_=3;
     1852      break;
     1853    }
     1854    // Check event
     1855    {
     1856      int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
     1857      if (status>=0) {
     1858        problemStatus_=5;
     1859        secondaryStatus_=ClpEventHandler::endOfFactorization;
     1860        break;
     1861      }
     1862    }
     1863    // Do iterations
     1864    if (canTryQuick) {
     1865      double * saveDuals=NULL;
     1866      ((ClpSimplexDual *)this)->whileIterating(saveDuals,0);
     1867    } else {
     1868      whileIterating(startingTheta,  endingTheta, reportIncrement,
     1869                     changeLower, changeUpper,
     1870                     changeObjective);
     1871    }
     1872  }
     1873  if (!problemStatus_) {
     1874    theta_=change+startingTheta;
     1875    eventHandler_->event(ClpEventHandler::theta);
     1876    return 0;
     1877  } else if (problemStatus_==10) {
     1878    return -1;
     1879  } else {
     1880    return problemStatus_;
     1881  }
     1882}
     1883/* Checks if finished.  Updates status */
     1884void
     1885ClpSimplexOther::statusOfProblemInParametrics(int type, ClpDataSave & saveData)
     1886{
     1887  if (type==2) {
     1888    // trouble - go to recovery
     1889    problemStatus_=10;
     1890    return;
     1891  }
     1892  if (problemStatus_>-3||factorization_->pivots()) {
     1893    // factorize
     1894    // later on we will need to recover from singularities
     1895    // also we could skip if first time
     1896    if (type) {
     1897      // is factorization okay?
     1898      if (internalFactorize(1)) {
     1899        // trouble - go to recovery
     1900        problemStatus_=10;
     1901        return;
     1902      }
     1903    }
     1904    if (problemStatus_!=-4||factorization_->pivots()>10)
     1905      problemStatus_=-3;
     1906  }
     1907  // at this stage status is -3 or -4 if looks infeasible
     1908  // get primal and dual solutions
     1909  gutsOfSolution(NULL,NULL);
     1910  double realDualInfeasibilities=sumDualInfeasibilities_;
     1911  // If bad accuracy treat as singular
     1912  if ((largestPrimalError_>1.0e15||largestDualError_>1.0e15)&&numberIterations_) {
     1913    // trouble - go to recovery
     1914    problemStatus_=10;
     1915    return;
     1916  } else if (largestPrimalError_<1.0e-7&&largestDualError_<1.0e-7) {
     1917    // Can reduce tolerance
     1918    double newTolerance = CoinMax(0.99*factorization_->pivotTolerance(),saveData.pivotTolerance_);
     1919    factorization_->pivotTolerance(newTolerance);
     1920  }
     1921  // Check if looping
     1922  int loop;
     1923  if (type!=2)
     1924    loop = progress_->looping();
     1925  else
     1926    loop=-1;
     1927  if (loop>=0) {
     1928    problemStatus_ = loop; //exit if in loop
     1929    if (!problemStatus_) {
     1930      // declaring victory
     1931      numberPrimalInfeasibilities_ = 0;
     1932      sumPrimalInfeasibilities_ = 0.0;
     1933    } else {
     1934      problemStatus_ = 10; // instead - try other algorithm
     1935    }
     1936    return;
     1937  } else if (loop<-1) {
     1938    // something may have changed
     1939    gutsOfSolution(NULL,NULL);
     1940  }
     1941  progressFlag_ = 0; //reset progress flag
     1942  if (handler_->detail(CLP_SIMPLEX_STATUS,messages_)<100) {
     1943    handler_->message(CLP_SIMPLEX_STATUS,messages_)
     1944      <<numberIterations_<<objectiveValue();
     1945    handler_->printing(sumPrimalInfeasibilities_>0.0)
     1946      <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
     1947    handler_->printing(sumDualInfeasibilities_>0.0)
     1948      <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
     1949    handler_->printing(numberDualInfeasibilitiesWithoutFree_
     1950                       <numberDualInfeasibilities_)
     1951                         <<numberDualInfeasibilitiesWithoutFree_;
     1952    handler_->message()<<CoinMessageEol;
     1953  }
     1954  /* If we are primal feasible and any dual infeasibilities are on
     1955     free variables then it is better to go to primal */
     1956  if (!numberPrimalInfeasibilities_&&!numberDualInfeasibilitiesWithoutFree_&&
     1957      numberDualInfeasibilities_) {
     1958    problemStatus_=10;
     1959    return;
     1960  }
     1961 
     1962  // check optimal
     1963  // give code benefit of doubt
     1964  if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
     1965      sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
     1966    // say optimal (with these bounds etc)
     1967    numberDualInfeasibilities_ = 0;
     1968    sumDualInfeasibilities_ = 0.0;
     1969    numberPrimalInfeasibilities_ = 0;
     1970    sumPrimalInfeasibilities_ = 0.0;
     1971  }
     1972  if (dualFeasible()||problemStatus_==-4) {
     1973    progress_->modifyObjective(objectiveValue_
     1974                               -sumDualInfeasibilities_*dualBound_);
     1975  }
     1976  if (numberPrimalInfeasibilities_) {
     1977    if (problemStatus_==-4||problemStatus_==-5) {
     1978      problemStatus_=1; // infeasible
     1979    }
     1980  } else if (numberDualInfeasibilities_) {
     1981    // clean up
     1982    problemStatus_=10;
     1983  } else {
     1984    problemStatus_=0;
     1985  }
     1986  lastGoodIteration_ = numberIterations_;
     1987  if (problemStatus_<0) {
     1988    sumDualInfeasibilities_=realDualInfeasibilities; // back to say be careful
     1989    if (sumDualInfeasibilities_)
     1990      numberDualInfeasibilities_=1;
     1991  }
     1992}
     1993/* This has the flow between re-factorizations
     1994   Reasons to come out:
     1995   -1 iterations etc
     1996   -2 inaccuracy
     1997   -3 slight inaccuracy (and done iterations)
     1998   +0 looks optimal (might be unbounded - but we will investigate)
     1999   +1 looks infeasible
     2000   +3 max iterations
     2001*/
     2002int
     2003ClpSimplexOther::whileIterating(double startingTheta, double & endingTheta,double reportIncrement,
     2004                                const double * changeLower, const double * changeUpper,
     2005                                const double * changeObjective)
     2006{
     2007  {
     2008    int i;
     2009    for (i=0;i<4;i++) {
     2010      rowArray_[i]->clear();
     2011    }   
     2012    for (i=0;i<2;i++) {
     2013      columnArray_[i]->clear();
     2014    }   
     2015  }     
     2016  // if can't trust much and long way from optimal then relax
     2017  if (largestPrimalError_>10.0)
     2018    factorization_->relaxAccuracyCheck(CoinMin(1.0e2,largestPrimalError_/10.0));
     2019  else
     2020    factorization_->relaxAccuracyCheck(1.0);
     2021  // status stays at -1 while iterating, >=0 finished, -2 to invert
     2022  // status -3 to go to top without an invert
     2023  int returnCode = -1;
     2024  double saveSumDual = sumDualInfeasibilities_; // so we know to be careful
     2025  //double useTheta = startingTheta;
     2026  double * primalChange = new double[numberRows_];
     2027  double * dualChange = new double[numberColumns_];
     2028  int numberTotal = numberColumns_+numberRows_;
     2029
     2030  while (problemStatus_==-1) {
     2031
     2032    // Get theta for bounds - we know can't crossover
     2033    // get change
     2034    int iSequence;
     2035    for (iSequence=0;iSequence<numberTotal;iSequence++) {
     2036      primalChange[iSequence]=0.0;
     2037      switch(getStatus(iSequence)) {
     2038       
     2039      case basic:
     2040      case isFree:
     2041      case superBasic:
     2042        break;
     2043      case isFixed:
     2044      case atUpperBound:
     2045        primalChange[iSequence]=changeUpper[iSequence];
     2046        break;
     2047      case atLowerBound:
     2048        primalChange[iSequence]=changeLower[iSequence];
     2049        break;
     2050      }
     2051    }
     2052    // use rowActivity
     2053    memset(rowActivity_,0,numberRows_*sizeof(double));
     2054    times(1.0,primalChange,rowActivity_);
     2055   
     2056    // choose row to go out
     2057    // dualRow will go to virtual row pivot choice algorithm
     2058    ((ClpSimplexDual *) this)->dualRow(-1);
     2059    if (pivotRow_>=0) {
     2060      // we found a pivot row
     2061      if (handler_->detail(CLP_SIMPLEX_PIVOTROW,messages_)<100) {
     2062        handler_->message(CLP_SIMPLEX_PIVOTROW,messages_)
     2063          <<pivotRow_
     2064          <<CoinMessageEol;
     2065      }
     2066      // check accuracy of weights
     2067      dualRowPivot_->checkAccuracy();
     2068      // Get good size for pivot
     2069      // Allow first few iterations to take tiny
     2070      double acceptablePivot=1.0e-9;
     2071      if (numberIterations_>100)
     2072        acceptablePivot=1.0e-8;
     2073      if (factorization_->pivots()>10||
     2074          (factorization_->pivots()&&saveSumDual))
     2075        acceptablePivot=1.0e-5; // if we have iterated be more strict
     2076      else if (factorization_->pivots()>5)
     2077        acceptablePivot=1.0e-6; // if we have iterated be slightly more strict
     2078      else if (factorization_->pivots())
     2079        acceptablePivot=1.0e-8; // relax
     2080      double bestPossiblePivot=1.0;
     2081      // get sign for finding row of tableau
     2082      // normal iteration
     2083      // create as packed
     2084      double direction=directionOut_;
     2085      rowArray_[0]->createPacked(1,&pivotRow_,&direction);
     2086      factorization_->updateColumnTranspose(rowArray_[1],rowArray_[0]);
     2087      // put row of tableau in rowArray[0] and columnArray[0]
     2088      matrix_->transposeTimes(this,-1.0,
     2089                              rowArray_[0],rowArray_[3],columnArray_[0]);
     2090      // do ratio test for normal iteration
     2091      bestPossiblePivot = ((ClpSimplexDual *) this)->dualColumn(rowArray_[0],
     2092                                                                columnArray_[0],columnArray_[1],
     2093                                                                rowArray_[3],acceptablePivot,NULL);
     2094      if (sequenceIn_>=0) {
     2095        // normal iteration
     2096        // update the incoming column
     2097        double btranAlpha = -alpha_*directionOut_; // for check
     2098        unpackPacked(rowArray_[1]);
     2099        factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
     2100        // and update dual weights (can do in parallel - with extra array)
     2101        alpha_ = dualRowPivot_->updateWeights(rowArray_[0],
     2102                                              rowArray_[2],
     2103                                              rowArray_[1]);
     2104        // see if update stable
     2105#ifdef CLP_DEBUG
     2106        if ((handler_->logLevel()&32))
     2107          printf("btran alpha %g, ftran alpha %g\n",btranAlpha,alpha_);
     2108#endif
     2109        double checkValue=1.0e-7;
     2110        // if can't trust much and long way from optimal then relax
     2111        if (largestPrimalError_>10.0)
     2112          checkValue = CoinMin(1.0e-4,1.0e-8*largestPrimalError_);
     2113        if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
     2114            fabs(btranAlpha-alpha_)>checkValue*(1.0+fabs(alpha_))) {
     2115          handler_->message(CLP_DUAL_CHECK,messages_)
     2116            <<btranAlpha
     2117            <<alpha_
     2118            <<CoinMessageEol;
     2119          if (factorization_->pivots()) {
     2120            dualRowPivot_->unrollWeights();
     2121            problemStatus_=-2; // factorize now
     2122            rowArray_[0]->clear();
     2123            rowArray_[1]->clear();
     2124            columnArray_[0]->clear();
     2125            returnCode=-2;
     2126            break;
     2127          } else {
     2128            // take on more relaxed criterion
     2129            double test;
     2130            if (fabs(btranAlpha)<1.0e-8||fabs(alpha_)<1.0e-8)
     2131              test = 1.0e-1*fabs(alpha_);
     2132            else
     2133              test = 1.0e-4*(1.0+fabs(alpha_));
     2134            if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
     2135                fabs(btranAlpha-alpha_)>test) {
     2136              dualRowPivot_->unrollWeights();
     2137              // need to reject something
     2138              char x = isColumn(sequenceOut_) ? 'C' :'R';
     2139              handler_->message(CLP_SIMPLEX_FLAG,messages_)
     2140                <<x<<sequenceWithin(sequenceOut_)
     2141                <<CoinMessageEol;
     2142              setFlagged(sequenceOut_);
     2143              progress_->clearBadTimes();
     2144              lastBadIteration_ = numberIterations_; // say be more cautious
     2145              rowArray_[0]->clear();
     2146              rowArray_[1]->clear();
     2147              columnArray_[0]->clear();
     2148              if (fabs(alpha_)<1.0e-10&&fabs(btranAlpha)<1.0e-8&&numberIterations_>100) {
     2149                //printf("I think should declare infeasible\n");
     2150                problemStatus_=1;
     2151                returnCode=1;
     2152                break;
     2153              }
     2154              continue;
     2155            }
     2156          }
     2157        }
     2158        // update duals BEFORE replaceColumn so can do updateColumn
     2159        double objectiveChange=0.0;
     2160        // do duals first as variables may flip bounds
     2161        // rowArray_[0] and columnArray_[0] may have flips
     2162        // so use rowArray_[3] for work array from here on
     2163        int nswapped = 0;
     2164        //rowArray_[0]->cleanAndPackSafe(1.0e-60);
     2165        //columnArray_[0]->cleanAndPackSafe(1.0e-60);
     2166        nswapped = ((ClpSimplexDual *) this)->updateDualsInDual(rowArray_[0],columnArray_[0],
     2167                                     rowArray_[2],theta_,
     2168                                     objectiveChange,false);
     2169
     2170        // which will change basic solution
     2171        if (nswapped) {
     2172          factorization_->updateColumn(rowArray_[3],rowArray_[2]);
     2173          dualRowPivot_->updatePrimalSolution(rowArray_[2],
     2174                                              1.0,objectiveChange);
     2175          // recompute dualOut_
     2176          valueOut_ = solution_[sequenceOut_];
     2177          if (directionOut_<0) {
     2178            dualOut_ = valueOut_ - upperOut_;
     2179          } else {
     2180            dualOut_ = lowerOut_ - valueOut_;
     2181          }
     2182        }
     2183        // amount primal will move
     2184        double movement = -dualOut_*directionOut_/alpha_;
     2185        // so objective should increase by fabs(dj)*movement
     2186        // but we already have objective change - so check will be good
     2187        if (objectiveChange+fabs(movement*dualIn_)<-1.0e-5) {
     2188#ifdef CLP_DEBUG
     2189          if (handler_->logLevel()&32)
     2190            printf("movement %g, swap change %g, rest %g  * %g\n",
     2191                   objectiveChange+fabs(movement*dualIn_),
     2192                   objectiveChange,movement,dualIn_);
     2193#endif
     2194          if(factorization_->pivots()) {
     2195            // going backwards - factorize
     2196            dualRowPivot_->unrollWeights();
     2197            problemStatus_=-2; // factorize now
     2198            returnCode=-2;
     2199            break;
     2200          }
     2201        }
     2202        CoinAssert(fabs(dualOut_)<1.0e50);
     2203        // if stable replace in basis
     2204        int updateStatus = factorization_->replaceColumn(this,
     2205                                                         rowArray_[2],
     2206                                                         rowArray_[1],
     2207                                                         pivotRow_,
     2208                                                         alpha_);
     2209        // if no pivots, bad update but reasonable alpha - take and invert
     2210        if (updateStatus==2&&
     2211                   !factorization_->pivots()&&fabs(alpha_)>1.0e-5)
     2212          updateStatus=4;
     2213        if (updateStatus==1||updateStatus==4) {
     2214          // slight error
     2215          if (factorization_->pivots()>5||updateStatus==4) {
     2216            problemStatus_=-2; // factorize now
     2217            returnCode=-3;
     2218          }
     2219        } else if (updateStatus==2) {
     2220          // major error
     2221          dualRowPivot_->unrollWeights();
     2222          // later we may need to unwind more e.g. fake bounds
     2223          if (factorization_->pivots()) {
     2224            problemStatus_=-2; // factorize now
     2225            returnCode=-2;
     2226            break;
     2227          } else {
     2228            // need to reject something
     2229            char x = isColumn(sequenceOut_) ? 'C' :'R';
     2230            handler_->message(CLP_SIMPLEX_FLAG,messages_)
     2231              <<x<<sequenceWithin(sequenceOut_)
     2232              <<CoinMessageEol;
     2233            setFlagged(sequenceOut_);
     2234            progress_->clearBadTimes();
     2235            lastBadIteration_ = numberIterations_; // say be more cautious
     2236            rowArray_[0]->clear();
     2237            rowArray_[1]->clear();
     2238            columnArray_[0]->clear();
     2239            // make sure dual feasible
     2240            // look at all rows and columns
     2241            double objectiveChange=0.0;
     2242            ((ClpSimplexDual *) this)->updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[1],
     2243                              0.0,objectiveChange,true);
     2244            continue;
     2245          }
     2246        } else if (updateStatus==3) {
     2247          // out of memory
     2248          // increase space if not many iterations
     2249          if (factorization_->pivots()<
     2250              0.5*factorization_->maximumPivots()&&
     2251              factorization_->pivots()<200)
     2252            factorization_->areaFactor(
     2253                                       factorization_->areaFactor() * 1.1);
     2254          problemStatus_=-2; // factorize now
     2255        } else if (updateStatus==5) {
     2256          problemStatus_=-2; // factorize now
     2257        }
     2258        // update primal solution
     2259        if (theta_<0.0) {
     2260#ifdef CLP_DEBUG
     2261          if (handler_->logLevel()&32)
     2262            printf("negative theta %g\n",theta_);
     2263#endif
     2264          theta_=0.0;
     2265        }
     2266        // do actual flips
     2267        ((ClpSimplexDual *) this)->flipBounds(rowArray_[0],columnArray_[0],theta_);
     2268        //rowArray_[1]->expand();
     2269        dualRowPivot_->updatePrimalSolution(rowArray_[1],
     2270                                            movement,
     2271                                            objectiveChange);
     2272        // modify dualout
     2273        dualOut_ /= alpha_;
     2274        dualOut_ *= -directionOut_;
     2275        //setStatus(sequenceIn_,basic);
     2276        dj_[sequenceIn_]=0.0;
     2277        double oldValue=valueIn_;
     2278        if (directionIn_==-1) {
     2279          // as if from upper bound
     2280          valueIn_ = upperIn_+dualOut_;
     2281        } else {
     2282          // as if from lower bound
     2283          valueIn_ = lowerIn_+dualOut_;
     2284        }
     2285        objectiveChange += cost_[sequenceIn_]*(valueIn_-oldValue);
     2286        // outgoing
     2287        // set dj to zero unless values pass
     2288        if (directionOut_>0) {
     2289          valueOut_ = lowerOut_;
     2290          dj_[sequenceOut_] = theta_;
     2291        } else {
     2292          valueOut_ = upperOut_;
     2293          dj_[sequenceOut_] = -theta_;
     2294        }
     2295        solution_[sequenceOut_]=valueOut_;
     2296        int whatNext=housekeeping(objectiveChange);
     2297        // and set bounds correctly
     2298        ((ClpSimplexDual *) this)->originalBound(sequenceIn_);
     2299        ((ClpSimplexDual *) this)->changeBound(sequenceOut_);
     2300        if (whatNext==1) {
     2301          problemStatus_ =-2; // refactorize
     2302        } else if (whatNext==2) {
     2303          // maximum iterations or equivalent
     2304          problemStatus_= 3;
     2305          returnCode=3;
     2306          break;
     2307        }
     2308        // Check event
     2309        {
     2310          int status = eventHandler_->event(ClpEventHandler::endOfIteration);
     2311          if (status>=0) {
     2312            problemStatus_=5;
     2313            secondaryStatus_=ClpEventHandler::endOfIteration;
     2314            returnCode=4;
     2315            break;
     2316          }
     2317        }
     2318      } else {
     2319        // no incoming column is valid
     2320        pivotRow_=-1;
     2321#ifdef CLP_DEBUG
     2322        if (handler_->logLevel()&32)
     2323          printf("** no column pivot\n");
     2324#endif
     2325        if (factorization_->pivots()<5) {
     2326          // If not in branch and bound etc save ray
     2327          if ((specialOptions_&(1024|4096))==0) {
     2328            // create ray anyway
     2329            delete [] ray_;
     2330            ray_ = new double [ numberRows_];
     2331            rowArray_[0]->expand(); // in case packed
     2332            ClpDisjointCopyN(rowArray_[0]->denseVector(),numberRows_,ray_);
     2333          }
     2334          // If we have just factorized and infeasibility reasonable say infeas
     2335          if (((specialOptions_&4096)!=0||bestPossiblePivot<1.0e-11)&&dualBound_>1.0e8) {
     2336            if (valueOut_>upperOut_+1.0e-3||valueOut_<lowerOut_-1.0e-3
     2337                || (specialOptions_&64)==0) {
     2338              // say infeasible
     2339              problemStatus_=1;
     2340              // unless primal feasible!!!!
     2341              //printf("%d %g %d %g\n",numberPrimalInfeasibilities_,sumPrimalInfeasibilities_,
     2342              //     numberDualInfeasibilities_,sumDualInfeasibilities_);
     2343              if (numberDualInfeasibilities_)
     2344                problemStatus_=10;
     2345              rowArray_[0]->clear();
     2346              columnArray_[0]->clear();
     2347              returnCode=1;
     2348              break;
     2349            }
     2350          }
     2351          // If special option set - put off as long as possible
     2352          if ((specialOptions_&64)==0) {
     2353            problemStatus_=-4; //say looks infeasible
     2354          } else {
     2355            // flag
     2356            char x = isColumn(sequenceOut_) ? 'C' :'R';
     2357            handler_->message(CLP_SIMPLEX_FLAG,messages_)
     2358              <<x<<sequenceWithin(sequenceOut_)
     2359              <<CoinMessageEol;
     2360            setFlagged(sequenceOut_);
     2361            if (!factorization_->pivots()) {
     2362              rowArray_[0]->clear();
     2363              columnArray_[0]->clear();
     2364              continue;
     2365            }
     2366          }
     2367        }
     2368        rowArray_[0]->clear();
     2369        columnArray_[0]->clear();
     2370        returnCode=1;
     2371        break;
     2372      }
     2373    } else {
     2374      // no pivot row
     2375#ifdef CLP_DEBUG
     2376      if (handler_->logLevel()&32)
     2377        printf("** no row pivot\n");
     2378#endif
     2379      int numberPivots = factorization_->pivots();
     2380      bool specialCase;
     2381      int useNumberFake;
     2382      returnCode=0;
     2383      if (numberPivots<20&&
     2384          (specialOptions_&2048)!=0&&!numberChanged_&&perturbation_>=100
     2385          &&dualBound_>1.0e8) {
     2386        specialCase=true;
     2387        // as dual bound high - should be okay
     2388        useNumberFake=0;
     2389      } else {
     2390        specialCase=false;
     2391        useNumberFake=numberFake_;
     2392      }
     2393      if (!numberPivots||specialCase) {
     2394        // may have crept through - so may be optimal
     2395        // check any flagged variables
     2396        int iRow;
     2397        for (iRow=0;iRow<numberRows_;iRow++) {
     2398          int iPivot=pivotVariable_[iRow];
     2399          if (flagged(iPivot))
     2400            break;
     2401        }
     2402        if (iRow<numberRows_&&numberPivots) {
     2403          // try factorization
     2404          returnCode=-2;
     2405        }
     2406       
     2407        if (useNumberFake||numberDualInfeasibilities_) {
     2408          // may be dual infeasible
     2409          problemStatus_=-5;
     2410        } else {
     2411          if (iRow<numberRows_) {
     2412            problemStatus_=-5;
     2413          } else {
     2414            if (numberPivots) {
     2415              // objective may be wrong
     2416              objectiveValue_ = innerProduct(cost_,
     2417                                                                        numberColumns_+numberRows_,
     2418                                                                        solution_);
     2419              objectiveValue_ += objective_->nonlinearOffset();
     2420              objectiveValue_ /= (objectiveScale_*rhsScale_);
     2421              if ((specialOptions_&16384)==0) {
     2422                // and dual_ may be wrong (i.e. for fixed or basic)
     2423                CoinIndexedVector * arrayVector = rowArray_[1];
     2424                arrayVector->clear();
     2425                int iRow;
     2426                double * array = arrayVector->denseVector();
     2427                /* Use dual_ instead of array
     2428                   Even though dual_ is only numberRows_ long this is
     2429                   okay as gets permuted to longer rowArray_[2]
     2430                */
     2431                arrayVector->setDenseVector(dual_);
     2432                int * index = arrayVector->getIndices();
     2433                int number=0;
     2434                for (iRow=0;iRow<numberRows_;iRow++) {
     2435                  int iPivot=pivotVariable_[iRow];
     2436                  double value = cost_[iPivot];
     2437                  dual_[iRow]=value;
     2438                  if (value) {
     2439                    index[number++]=iRow;
     2440                  }
     2441                }
     2442                arrayVector->setNumElements(number);
     2443                // Extended duals before "updateTranspose"
     2444                matrix_->dualExpanded(this,arrayVector,NULL,0);
     2445                // Btran basic costs
     2446                rowArray_[2]->clear();
     2447                factorization_->updateColumnTranspose(rowArray_[2],arrayVector);
     2448                // and return vector
     2449                arrayVector->setDenseVector(array);
     2450              }
     2451            }
     2452            problemStatus_=0;
     2453            sumPrimalInfeasibilities_=0.0;
     2454            if ((specialOptions_&(1024+16384))!=0) {
     2455              CoinIndexedVector * arrayVector = rowArray_[1];
     2456              arrayVector->clear();
     2457              double * rhs = arrayVector->denseVector();
     2458              times(1.0,solution_,rhs);
     2459              bool bad2=false;
     2460              int i;
     2461              for ( i=0;i<numberRows_;i++) {
     2462                if (rhs[i]<rowLowerWork_[i]-primalTolerance_||
     2463                    rhs[i]>rowUpperWork_[i]+primalTolerance_) {
     2464                  bad2=true;
     2465                } else if (fabs(rhs[i]-rowActivityWork_[i])>1.0e-3) {
     2466                }
     2467                rhs[i]=0.0;
     2468              }
     2469              for ( i=0;i<numberColumns_;i++) {
     2470                if (solution_[i]<columnLowerWork_[i]-primalTolerance_||
     2471                    solution_[i]>columnUpperWork_[i]+primalTolerance_) {
     2472                  bad2=true;
     2473                }
     2474              }
     2475              if (bad2) {
     2476                problemStatus_=-3;
     2477                returnCode=-2;
     2478                // Force to re-factorize early next time
     2479                int numberPivots = factorization_->pivots();
     2480                forceFactorization_=CoinMin(forceFactorization_,(numberPivots+1)>>1);
     2481              }
     2482            }
     2483          }
     2484        }
     2485      } else {
     2486        problemStatus_=-3;
     2487        returnCode=-2;
     2488        // Force to re-factorize early next time
     2489        int numberPivots = factorization_->pivots();
     2490        forceFactorization_=CoinMin(forceFactorization_,(numberPivots+1)>>1);
     2491      }
     2492      break;
     2493    }
     2494  }
     2495  delete [] primalChange;
     2496  delete [] dualChange;
     2497  return returnCode;
     2498}
  • trunk/include/ClpEventHandler.hpp

    r635 r636  
    3535    node, // for Cbc
    3636    treeStatus, // for Cbc
    37     solution // for Cbc
     37    solution, // for Cbc
     38    theta // hit in parametrics
    3839  };
    3940  /**@name Virtual method that the derived classe should provide.
  • trunk/include/ClpSimplex.hpp

    r627 r636  
    105105      Only to be used with mini constructor */
    106106  void originalModel(ClpSimplex * miniModel);
     107  /**
     108     If you are re-using the same matrix again and again then the setup time
     109     to do scaling may be significant.  Also you may not want to initialize all values
     110     or return all values (especially if infeasible).  While an auxiliary model exists
     111     it will be faster.  If options -1 then model is switched off.  Otherwise switched on
     112     with following options.
     113     1 - rhs is constant
     114     2 - bounds are constant
     115     4 - objective is constant
     116     8 - solution in by basis and no djs etc in
     117     16 - no duals out (but reduced costs)
     118     32 - no output if infeasible
     119  */
     120  void auxiliaryModel(int options);
    107121  /// Assignment operator. This copies the data
    108122    ClpSimplex & operator=(const ClpSimplex & rhs);
     
    11381152  /// Column activities - working copy
    11391153  double * columnActivityWork_;
     1154  /// Auxiliary model
     1155  ClpSimplex * auxiliaryModel_;
    11401156  /// Number of dual infeasibilities
    11411157  int numberDualInfeasibilities_;
  • trunk/include/ClpSimplexDual.hpp

    r627 r636  
    280280  /** Get next free , -1 if none */
    281281  int nextSuperBasic();
     282  /** Startup part of dual (may be extended to other algorithms)
     283      returns 0 if good, 1 if bad */
     284  int startupSolve(int ifValuesPass,double * saveDuals,int startFinishOptions);
     285  void finishSolve(int startFinishOptions);
     286  void gutsOfDual(int ifValuesPass,double * & saveDuals,int initialStatus,
     287                  ClpDataSave & saveData);
     288  //int dual2(int ifValuesPass,int startFinishOptions=0);
    282289 
    283290  //@}
  • trunk/include/ClpSimplexOther.hpp

    r631 r636  
    5252                  double * valueIncrease, int * sequenceIncrease,
    5353                  double * valueDecrease, int * sequenceDecrease);
     54  /** Parametrics
     55      This is an initial slow version.
     56      The code uses current bounds + theta * change (if change array not NULL)
     57      and similarly for objective.
     58      It starts at startingTheta and returns ending theta in endingTheta.
     59      If reportIncrement 0.0 it will report on any movement
     60      If reportIncrement >0.0 it will report at startingTheta+k*reportIncrement.
     61      If it can not reach input endingTheta return code will be 1 for infeasible,
     62      2 for unbounded, if error on ranges -1,  otherwise 0.
     63      Normal report is just theta and objective but
     64      if event handler exists it may do more
     65      On exit endingTheta is maximum reached (can be used for next startingTheta)
     66  */
     67  int parametrics(double startingTheta, double & endingTheta,double reportIncrement,
     68                  const double * changeLowerBound, const double * changeUpperBound,
     69                  const double * changeLowerRhs, const double * changeUpperRhs,
     70                  const double * changeObjective);
     71  /** Parametrics - inner loop
     72      This first attempt is when reportIncrement non zero and may
     73      not report endingTheta correctly
     74      If it can not reach input endingTheta return code will be 1 for infeasible,
     75      2 for unbounded,  otherwise 0.
     76      Normal report is just theta and objective but
     77      if event handler exists it may do more
     78  */
     79  int parametricsLoop(double startingTheta, double & endingTheta,double reportIncrement,
     80                      const double * changeLower, const double * changeUpper,
     81                      const double * changeObjective, ClpDataSave & data,
     82                      bool canTryQuick);
     83  /**  Refactorizes if necessary
     84       Checks if finished.  Updates status.
     85
     86       type - 0 initial so set up save arrays etc
     87            - 1 normal -if good update save
     88            - 2 restoring from saved
     89  */
     90  void statusOfProblemInParametrics(int type,ClpDataSave & saveData);
     91  /** This has the flow between re-factorizations
     92
     93      Reasons to come out:
     94      -1 iterations etc
     95      -2 inaccuracy
     96      -3 slight inaccuracy (and done iterations)
     97      +0 looks optimal (might be unbounded - but we will investigate)
     98      +1 looks infeasible
     99      +3 max iterations
     100   */
     101  int whileIterating(double startingTheta, double & endingTheta,double reportIncrement,
     102                      const double * changeLower, const double * changeUpper,
     103                      const double * changeObjective);
    54104  /**
    55105      Row array has row part of pivot row
     
    94144  */
    95145  ClpSimplex * crunch(double * rhs, int * whichRows, int * whichColumns,
    96                       int & nBound, bool moreBounds=false);
     146                      int & nBound, bool moreBounds=false, bool tightenBounds=false);
    97147  /** After very cursory presolve.
    98148      rhs is numberRows, whichRows is 3*numberRows and whichColumns is 2*numberColumns.
     
    101151                   const int * whichRows, const int * whichColumns,
    102152                   int nBound);
     153  /** Tightens integer bounds - returns number tightened or -1 if infeasible
     154  */
     155  int tightenIntegerBounds(double * rhsSpace);
    103156  //@}
    104157};
Note: See TracChangeset for help on using the changeset viewer.