source: trunk/Clp/src/ClpInterior.cpp @ 1321

Last change on this file since 1321 was 1321, checked in by forrest, 11 years ago

out compiler warnings and stability improvements

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 34.9 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4
5
6
7#include "CoinPragma.hpp"
8
9#include <math.h>
10
11#include "CoinHelperFunctions.hpp"
12#include "ClpInterior.hpp"
13#include "ClpMatrixBase.hpp"
14#ifdef COIN_DO_PDCO
15#include "ClpLsqr.hpp"
16#include "ClpPdcoBase.hpp"
17#endif
18#include "CoinDenseVector.hpp"
19#include "ClpMessage.hpp"
20#include "ClpHelperFunctions.hpp"
21#include "ClpCholeskyDense.hpp"
22#include "ClpLinearObjective.hpp"
23#include "ClpQuadraticObjective.hpp"
24#include <cfloat>
25
26#include <string>
27#include <stdio.h>
28#include <iostream>
29//#############################################################################
30
31ClpInterior::ClpInterior () :
32
33  ClpModel(),
34  largestPrimalError_(0.0),
35  largestDualError_(0.0),
36  sumDualInfeasibilities_(0.0),
37  sumPrimalInfeasibilities_(0.0),
38  worstComplementarity_(0.0),
39  xsize_(0.0),
40  zsize_(0.0),
41  lower_(NULL),
42  rowLowerWork_(NULL),
43  columnLowerWork_(NULL),
44  upper_(NULL),
45  rowUpperWork_(NULL),
46  columnUpperWork_(NULL),
47  cost_(NULL),
48  rhs_(NULL),
49   x_(NULL),
50   y_(NULL),
51  dj_(NULL),
52  lsqrObject_(NULL),
53  pdcoStuff_(NULL),
54  mu_(0.0),
55  objectiveNorm_(1.0e-12),
56  rhsNorm_(1.0e-12),
57  solutionNorm_(1.0e-12),
58  dualObjective_(0.0),
59  primalObjective_(0.0),
60  diagonalNorm_(1.0e-12),
61  stepLength_(0.995),
62  linearPerturbation_(1.0e-12),
63  diagonalPerturbation_(1.0e-15),
64  gamma_(0.0),
65  delta_(0),
66  targetGap_(1.0e-12),
67  projectionTolerance_(1.0e-7),
68  maximumRHSError_(0.0),
69  maximumBoundInfeasibility_(0.0),
70  maximumDualError_(0.0),
71  diagonalScaleFactor_(0.0),
72  scaleFactor_(1.0),
73  actualPrimalStep_(0.0),
74  actualDualStep_(0.0),
75  smallestInfeasibility_(0.0),
76  complementarityGap_(0.0),
77  baseObjectiveNorm_(0.0),
78  worstDirectionAccuracy_(0.0),
79  maximumRHSChange_(0.0),
80  errorRegion_(NULL),
81  rhsFixRegion_(NULL),
82  upperSlack_(NULL),
83  lowerSlack_(NULL),
84  diagonal_(NULL),
85  solution_(NULL),
86  workArray_(NULL),
87  deltaX_(NULL),
88  deltaY_(NULL),
89  deltaZ_(NULL),
90  deltaW_(NULL),
91  deltaSU_(NULL),
92  deltaSL_(NULL),
93  primalR_(NULL),
94  dualR_(NULL),
95  rhsB_(NULL),
96  rhsU_(NULL),
97  rhsL_(NULL),
98  rhsZ_(NULL),
99  rhsW_(NULL),
100  rhsC_(NULL),
101  zVec_(NULL),
102  wVec_(NULL),
103  cholesky_(NULL),
104  numberComplementarityPairs_(0),
105  numberComplementarityItems_(0),
106  maximumBarrierIterations_(200),
107  gonePrimalFeasible_(false),
108  goneDualFeasible_(false),
109  algorithm_(-1)
110{
111  memset(historyInfeasibility_,0,LENGTH_HISTORY*sizeof(double));
112  solveType_=3; // say interior based life form
113  cholesky_ = new ClpCholeskyDense(); // put in placeholder
114}
115
116// Subproblem constructor
117ClpInterior::ClpInterior ( const ClpModel * rhs,
118                           int numberRows, const int * whichRow,
119                           int numberColumns, const int * whichColumn,
120                           bool dropNames, bool dropIntegers)
121  : ClpModel(rhs, numberRows, whichRow,
122             numberColumns,whichColumn,dropNames,dropIntegers),
123    largestPrimalError_(0.0),
124    largestDualError_(0.0),
125    sumDualInfeasibilities_(0.0),
126    sumPrimalInfeasibilities_(0.0),
127    worstComplementarity_(0.0),
128    xsize_(0.0),
129    zsize_(0.0),
130    lower_(NULL),
131    rowLowerWork_(NULL),
132    columnLowerWork_(NULL),
133    upper_(NULL),
134    rowUpperWork_(NULL),
135    columnUpperWork_(NULL),
136    cost_(NULL),
137    rhs_(NULL),
138    x_(NULL),
139    y_(NULL),
140    dj_(NULL),
141    lsqrObject_(NULL),
142    pdcoStuff_(NULL),
143    mu_(0.0),
144    objectiveNorm_(1.0e-12),
145    rhsNorm_(1.0e-12),
146    solutionNorm_(1.0e-12),
147    dualObjective_(0.0),
148    primalObjective_(0.0),
149    diagonalNorm_(1.0e-12),
150    stepLength_(0.99995),
151    linearPerturbation_(1.0e-12),
152    diagonalPerturbation_(1.0e-15),
153    gamma_(0.0),
154    delta_(0),
155    targetGap_(1.0e-12),
156    projectionTolerance_(1.0e-7),
157    maximumRHSError_(0.0),
158    maximumBoundInfeasibility_(0.0),
159    maximumDualError_(0.0),
160    diagonalScaleFactor_(0.0),
161    scaleFactor_(0.0),
162    actualPrimalStep_(0.0),
163    actualDualStep_(0.0),
164    smallestInfeasibility_(0.0),
165    complementarityGap_(0.0),
166    baseObjectiveNorm_(0.0),
167    worstDirectionAccuracy_(0.0),
168    maximumRHSChange_(0.0),
169    errorRegion_(NULL),
170    rhsFixRegion_(NULL),
171    upperSlack_(NULL),
172    lowerSlack_(NULL),
173    diagonal_(NULL),
174    solution_(NULL),
175    workArray_(NULL),
176    deltaX_(NULL),
177    deltaY_(NULL),
178    deltaZ_(NULL),
179    deltaW_(NULL),
180    deltaSU_(NULL),
181    deltaSL_(NULL),
182    primalR_(NULL),
183    dualR_(NULL),
184    rhsB_(NULL),
185    rhsU_(NULL),
186    rhsL_(NULL),
187    rhsZ_(NULL),
188    rhsW_(NULL),
189    rhsC_(NULL),
190    zVec_(NULL),
191    wVec_(NULL),
192    cholesky_(NULL),
193    numberComplementarityPairs_(0),
194    numberComplementarityItems_(0),
195    maximumBarrierIterations_(200),
196    gonePrimalFeasible_(false),
197    goneDualFeasible_(false),
198    algorithm_(-1)
199{
200  memset(historyInfeasibility_,0,LENGTH_HISTORY*sizeof(double));
201  solveType_=3; // say interior based life form
202  cholesky_= new ClpCholeskyDense();
203}
204
205//-----------------------------------------------------------------------------
206
207ClpInterior::~ClpInterior ()
208{
209  gutsOfDelete();
210}
211//#############################################################################
212/*
213   This does housekeeping
214*/
215int 
216ClpInterior::housekeeping()
217{
218  numberIterations_++;
219  return 0;
220}
221// Copy constructor.
222ClpInterior::ClpInterior(const ClpInterior &rhs) :
223  ClpModel(rhs),
224  largestPrimalError_(0.0),
225  largestDualError_(0.0),
226  sumDualInfeasibilities_(0.0),
227  sumPrimalInfeasibilities_(0.0),
228  worstComplementarity_(0.0),
229  xsize_(0.0),
230  zsize_(0.0),
231  lower_(NULL),
232  rowLowerWork_(NULL),
233  columnLowerWork_(NULL),
234  upper_(NULL),
235  rowUpperWork_(NULL),
236  columnUpperWork_(NULL),
237  cost_(NULL),
238  rhs_(NULL),
239   x_(NULL),
240   y_(NULL),
241  dj_(NULL),
242  lsqrObject_(NULL),
243  pdcoStuff_(NULL),
244  errorRegion_(NULL),
245  rhsFixRegion_(NULL),
246  upperSlack_(NULL),
247  lowerSlack_(NULL),
248  diagonal_(NULL),
249  solution_(NULL),
250  workArray_(NULL),
251  deltaX_(NULL),
252  deltaY_(NULL),
253  deltaZ_(NULL),
254  deltaW_(NULL),
255  deltaSU_(NULL),
256  deltaSL_(NULL),
257  primalR_(NULL),
258  dualR_(NULL),
259  rhsB_(NULL),
260  rhsU_(NULL),
261  rhsL_(NULL),
262  rhsZ_(NULL),
263  rhsW_(NULL),
264  rhsC_(NULL),
265  zVec_(NULL),
266  wVec_(NULL),
267  cholesky_(NULL)
268{
269  gutsOfDelete();
270  gutsOfCopy(rhs);
271  solveType_=3; // say interior based life form
272}
273// Copy constructor from model
274ClpInterior::ClpInterior(const ClpModel &rhs) :
275  ClpModel(rhs),
276  largestPrimalError_(0.0),
277  largestDualError_(0.0),
278  sumDualInfeasibilities_(0.0),
279  sumPrimalInfeasibilities_(0.0),
280  worstComplementarity_(0.0),
281  xsize_(0.0),
282  zsize_(0.0),
283  lower_(NULL),
284  rowLowerWork_(NULL),
285  columnLowerWork_(NULL),
286  upper_(NULL),
287  rowUpperWork_(NULL),
288  columnUpperWork_(NULL),
289  cost_(NULL),
290  rhs_(NULL),
291   x_(NULL),
292   y_(NULL),
293  dj_(NULL),
294  lsqrObject_(NULL),
295  pdcoStuff_(NULL),
296  mu_(0.0),
297  objectiveNorm_(1.0e-12),
298  rhsNorm_(1.0e-12),
299  solutionNorm_(1.0e-12),
300  dualObjective_(0.0),
301  primalObjective_(0.0),
302  diagonalNorm_(1.0e-12),
303  stepLength_(0.99995),
304  linearPerturbation_(1.0e-12),
305  diagonalPerturbation_(1.0e-15),
306  gamma_(0.0),
307  delta_(0),
308  targetGap_(1.0e-12),
309  projectionTolerance_(1.0e-7),
310  maximumRHSError_(0.0),
311  maximumBoundInfeasibility_(0.0),
312  maximumDualError_(0.0),
313  diagonalScaleFactor_(0.0),
314  scaleFactor_(0.0),
315  actualPrimalStep_(0.0),
316  actualDualStep_(0.0),
317  smallestInfeasibility_(0.0),
318  complementarityGap_(0.0),
319  baseObjectiveNorm_(0.0),
320  worstDirectionAccuracy_(0.0),
321  maximumRHSChange_(0.0),
322  errorRegion_(NULL),
323  rhsFixRegion_(NULL),
324  upperSlack_(NULL),
325  lowerSlack_(NULL),
326  diagonal_(NULL),
327  solution_(NULL),
328  workArray_(NULL),
329  deltaX_(NULL),
330  deltaY_(NULL),
331  deltaZ_(NULL),
332  deltaW_(NULL),
333  deltaSU_(NULL),
334  deltaSL_(NULL),
335  primalR_(NULL),
336  dualR_(NULL),
337  rhsB_(NULL),
338  rhsU_(NULL),
339  rhsL_(NULL),
340  rhsZ_(NULL),
341  rhsW_(NULL),
342  rhsC_(NULL),
343  zVec_(NULL),
344  wVec_(NULL),
345  cholesky_(NULL),
346  numberComplementarityPairs_(0),
347  numberComplementarityItems_(0),
348  maximumBarrierIterations_(200),
349  gonePrimalFeasible_(false),
350  goneDualFeasible_(false),
351  algorithm_(-1)
352{
353  memset(historyInfeasibility_,0,LENGTH_HISTORY*sizeof(double));
354  solveType_=3; // say interior based life form
355  cholesky_= new ClpCholeskyDense();
356}
357// Assignment operator. This copies the data
358ClpInterior & 
359ClpInterior::operator=(const ClpInterior & rhs)
360{
361  if (this != &rhs) {
362    gutsOfDelete();
363    ClpModel::operator=(rhs);
364    gutsOfCopy(rhs);
365  }
366  return *this;
367}
368void 
369ClpInterior::gutsOfCopy(const ClpInterior & rhs)
370{
371  lower_ = ClpCopyOfArray(rhs.lower_,numberColumns_+numberRows_);
372  rowLowerWork_ = lower_+numberColumns_;
373  columnLowerWork_ = lower_;
374  upper_ = ClpCopyOfArray(rhs.upper_,numberColumns_+numberRows_);
375  rowUpperWork_ = upper_+numberColumns_;
376  columnUpperWork_ = upper_;
377  //cost_ = ClpCopyOfArray(rhs.cost_,2*(numberColumns_+numberRows_));
378  cost_ = ClpCopyOfArray(rhs.cost_,numberColumns_);
379  rhs_ = ClpCopyOfArray(rhs.rhs_,numberRows_);
380   x_ = ClpCopyOfArray(rhs.x_,numberColumns_);
381   y_ = ClpCopyOfArray(rhs.y_,numberRows_);
382  dj_ = ClpCopyOfArray(rhs.dj_,numberColumns_+numberRows_);
383#ifdef COIN_DO_PDCO
384  lsqrObject_= new ClpLsqr(*rhs.lsqrObject_);
385  pdcoStuff_ = rhs.pdcoStuff_->clone();
386#endif
387  largestPrimalError_ = rhs.largestPrimalError_;
388  largestDualError_ = rhs.largestDualError_;
389  sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_;
390  sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
391  worstComplementarity_ = rhs.worstComplementarity_;
392  xsize_ = rhs.xsize_;
393  zsize_ = rhs.zsize_;
394  solveType_=rhs.solveType_;
395  mu_ = rhs.mu_;
396  objectiveNorm_ = rhs.objectiveNorm_;
397  rhsNorm_ = rhs.rhsNorm_;
398  solutionNorm_ = rhs.solutionNorm_;
399  dualObjective_ = rhs.dualObjective_;
400  primalObjective_ = rhs.primalObjective_;
401  diagonalNorm_ = rhs.diagonalNorm_;
402  stepLength_ = rhs.stepLength_;
403  linearPerturbation_ = rhs.linearPerturbation_;
404  diagonalPerturbation_ = rhs.diagonalPerturbation_;
405  gamma_=rhs.gamma_;
406  delta_=rhs.delta_;
407  targetGap_ = rhs.targetGap_;
408  projectionTolerance_ = rhs.projectionTolerance_;
409  maximumRHSError_ = rhs.maximumRHSError_;
410  maximumBoundInfeasibility_ = rhs.maximumBoundInfeasibility_;
411  maximumDualError_ = rhs.maximumDualError_;
412  diagonalScaleFactor_ = rhs.diagonalScaleFactor_;
413  scaleFactor_ = rhs.scaleFactor_;
414  actualPrimalStep_ = rhs.actualPrimalStep_;
415  actualDualStep_ = rhs.actualDualStep_;
416  smallestInfeasibility_ = rhs.smallestInfeasibility_;
417  complementarityGap_ = rhs.complementarityGap_;
418  baseObjectiveNorm_ = rhs.baseObjectiveNorm_;
419  worstDirectionAccuracy_ = rhs.worstDirectionAccuracy_;
420  maximumRHSChange_ = rhs.maximumRHSChange_;
421  errorRegion_ = ClpCopyOfArray(rhs.errorRegion_,numberRows_);
422  rhsFixRegion_ = ClpCopyOfArray(rhs.rhsFixRegion_,numberRows_);
423  deltaY_ = ClpCopyOfArray(rhs.deltaY_,numberRows_);
424  upperSlack_ = ClpCopyOfArray(rhs.upperSlack_,numberRows_+numberColumns_);
425  lowerSlack_ = ClpCopyOfArray(rhs.lowerSlack_,numberRows_+numberColumns_);
426  diagonal_ = ClpCopyOfArray(rhs.diagonal_,numberRows_+numberColumns_);
427  deltaX_ = ClpCopyOfArray(rhs.deltaX_,numberRows_+numberColumns_);
428  deltaZ_ = ClpCopyOfArray(rhs.deltaZ_,numberRows_+numberColumns_);
429  deltaW_ = ClpCopyOfArray(rhs.deltaW_,numberRows_+numberColumns_);
430  deltaSU_ = ClpCopyOfArray(rhs.deltaSU_,numberRows_+numberColumns_);
431  deltaSL_ = ClpCopyOfArray(rhs.deltaSL_,numberRows_+numberColumns_);
432  primalR_ = ClpCopyOfArray(rhs.primalR_,numberRows_+numberColumns_);
433  dualR_ = ClpCopyOfArray(rhs.dualR_,numberRows_+numberColumns_);
434  rhsB_ = ClpCopyOfArray(rhs.rhsB_,numberRows_);
435  rhsU_ = ClpCopyOfArray(rhs.rhsU_,numberRows_+numberColumns_);
436  rhsL_ = ClpCopyOfArray(rhs.rhsL_,numberRows_+numberColumns_);
437  rhsZ_ = ClpCopyOfArray(rhs.rhsZ_,numberRows_+numberColumns_);
438  rhsW_ = ClpCopyOfArray(rhs.rhsW_,numberRows_+numberColumns_);
439  rhsC_ = ClpCopyOfArray(rhs.rhsC_,numberRows_+numberColumns_);
440  solution_ = ClpCopyOfArray(rhs.solution_,numberRows_+numberColumns_);
441  workArray_ = ClpCopyOfArray(rhs.workArray_,numberRows_+numberColumns_);
442  zVec_ = ClpCopyOfArray(rhs.zVec_,numberRows_+numberColumns_);
443  wVec_ = ClpCopyOfArray(rhs.wVec_,numberRows_+numberColumns_);
444  cholesky_ = rhs.cholesky_->clone();
445  numberComplementarityPairs_ = rhs.numberComplementarityPairs_;
446  numberComplementarityItems_ = rhs.numberComplementarityItems_;
447  maximumBarrierIterations_ = rhs.maximumBarrierIterations_;
448  gonePrimalFeasible_ = rhs.gonePrimalFeasible_;
449  goneDualFeasible_ = rhs.goneDualFeasible_;
450  algorithm_ = rhs.algorithm_;
451}
452
453void 
454ClpInterior::gutsOfDelete()
455{
456  delete [] lower_;
457  lower_=NULL;
458  rowLowerWork_=NULL;
459  columnLowerWork_=NULL;
460  delete [] upper_;
461  upper_=NULL;
462  rowUpperWork_=NULL;
463  columnUpperWork_=NULL;
464  delete [] cost_;
465  cost_=NULL;
466  delete [] rhs_;
467  rhs_ = NULL;
468  delete [] x_;
469  x_ = NULL;
470  delete [] y_;
471  y_ = NULL;
472  delete [] dj_;
473  dj_ = NULL;
474#ifdef COIN_DO_PDCO
475  delete lsqrObject_;
476  lsqrObject_ = NULL;
477  //delete pdcoStuff_;
478  pdcoStuff_=NULL;
479#endif
480  delete [] errorRegion_;
481  errorRegion_ = NULL;
482  delete [] rhsFixRegion_;
483  rhsFixRegion_ = NULL;
484  delete [] deltaY_;
485  deltaY_ = NULL;
486  delete [] upperSlack_;
487  upperSlack_ = NULL;
488  delete [] lowerSlack_;
489  lowerSlack_ = NULL;
490  delete [] diagonal_;
491  diagonal_ = NULL;
492  delete [] deltaX_;
493  deltaX_ = NULL;
494  delete [] deltaZ_;
495  deltaZ_ = NULL;
496  delete [] deltaW_;
497  deltaW_ = NULL;
498  delete [] deltaSU_;
499  deltaSU_ = NULL;
500  delete [] deltaSL_;
501  deltaSL_ = NULL;
502  delete [] primalR_;
503  primalR_=NULL;
504  delete [] dualR_;
505  dualR_=NULL;
506  delete [] rhsB_;
507  rhsB_ = NULL;
508  delete [] rhsU_;
509  rhsU_ = NULL;
510  delete [] rhsL_;
511  rhsL_ = NULL;
512  delete [] rhsZ_;
513  rhsZ_ = NULL;
514  delete [] rhsW_;
515  rhsW_ = NULL;
516  delete [] rhsC_;
517  rhsC_ = NULL;
518  delete [] solution_;
519  solution_ = NULL;
520  delete [] workArray_;
521  workArray_ = NULL;
522  delete [] zVec_;
523  zVec_ = NULL;
524  delete [] wVec_;
525  wVec_ = NULL;
526  delete cholesky_;
527}
528bool
529ClpInterior::createWorkingData()
530{
531  bool goodMatrix=true;
532  //check matrix
533  if (!matrix_->allElementsInRange(this,1.0e-12,1.0e20)) {
534    problemStatus_=4;
535    goodMatrix= false;
536  }
537  int nTotal = numberRows_+numberColumns_;
538  delete [] solution_;
539  solution_ = new double[nTotal];
540  CoinMemcpyN(columnActivity_,  numberColumns_,solution_);
541  CoinMemcpyN(rowActivity_,     numberRows_,solution_+numberColumns_);
542  delete [] cost_;
543  cost_ = new double[nTotal];
544  int i;
545  double direction = optimizationDirection_*objectiveScale_;
546  // direction is actually scale out not scale in
547  if (direction)
548    direction = 1.0/direction;
549  const double * obj = objective();
550  for (i=0;i<numberColumns_;i++)
551    cost_[i] = direction*obj[i];
552  memset(cost_+numberColumns_,0,numberRows_*sizeof(double));
553  // do scaling if needed
554  if (scalingFlag_>0&&!rowScale_) {
555    if (matrix_->scale(this))
556      scalingFlag_=-scalingFlag_; // not scaled after all
557  }
558  delete [] lower_;
559  delete [] upper_;
560  lower_ = new double[nTotal];
561  upper_ = new double[nTotal];
562  rowLowerWork_ = lower_+numberColumns_;
563  columnLowerWork_ = lower_;
564  rowUpperWork_ = upper_+numberColumns_;
565  columnUpperWork_ = upper_;
566  CoinMemcpyN(rowLower_,numberRows_,rowLowerWork_);
567  CoinMemcpyN(rowUpper_,numberRows_,rowUpperWork_);
568  CoinMemcpyN(columnLower_,numberColumns_,columnLowerWork_);
569  CoinMemcpyN(columnUpper_,numberColumns_,columnUpperWork_);
570  // clean up any mismatches on infinity
571  for (i=0;i<numberColumns_;i++) {
572    if (columnLowerWork_[i]<-1.0e30)
573      columnLowerWork_[i] = -COIN_DBL_MAX;
574    if (columnUpperWork_[i]>1.0e30)
575      columnUpperWork_[i] = COIN_DBL_MAX;
576  }
577  // clean up any mismatches on infinity
578  for (i=0;i<numberRows_;i++) {
579    if (rowLowerWork_[i]<-1.0e30)
580      rowLowerWork_[i] = -COIN_DBL_MAX;
581    if (rowUpperWork_[i]>1.0e30)
582      rowUpperWork_[i] = COIN_DBL_MAX;
583  }
584  // check rim of problem okay
585  if (!sanityCheck())
586    goodMatrix=false;
587  if(rowScale_) {
588    for (i=0;i<numberColumns_;i++) {
589      double multiplier = rhsScale_/columnScale_[i];
590      cost_[i] *= columnScale_[i];
591      if (columnLowerWork_[i]>-1.0e50)
592        columnLowerWork_[i] *= multiplier;
593      if (columnUpperWork_[i]<1.0e50)
594        columnUpperWork_[i] *= multiplier;
595     
596    }
597    for (i=0;i<numberRows_;i++) {
598      double multiplier = rhsScale_*rowScale_[i];
599      if (rowLowerWork_[i]>-1.0e50)
600        rowLowerWork_[i] *= multiplier;
601      if (rowUpperWork_[i]<1.0e50)
602        rowUpperWork_[i] *= multiplier;
603    }
604  } else if (rhsScale_!=1.0) {
605    for (i=0;i<numberColumns_+numberRows_;i++) {
606      if (lower_[i]>-1.0e50)
607        lower_[i] *= rhsScale_;
608      if (upper_[i]<1.0e50)
609        upper_[i] *= rhsScale_;
610     
611    }
612  }
613  assert (!errorRegion_);
614  errorRegion_ = new double [numberRows_];
615  assert (!rhsFixRegion_);
616  rhsFixRegion_ = new double [numberRows_];
617  assert (!deltaY_);
618  deltaY_ = new double [numberRows_];
619  CoinZeroN(deltaY_,numberRows_);
620  assert (!upperSlack_);
621  upperSlack_ = new double [nTotal];
622  assert (!lowerSlack_);
623  lowerSlack_ = new double [nTotal];
624  assert (!diagonal_);
625  diagonal_ = new double [nTotal];
626  assert (!deltaX_);
627  deltaX_ = new double [nTotal];
628  CoinZeroN(deltaX_,nTotal);
629  assert (!deltaZ_);
630  deltaZ_ = new double [nTotal];
631  CoinZeroN(deltaZ_,nTotal);
632  assert (!deltaW_);
633  deltaW_ = new double [nTotal];
634  CoinZeroN(deltaW_,nTotal);
635  assert (!deltaSU_);
636  deltaSU_ = new double [nTotal];
637  CoinZeroN(deltaSU_,nTotal);
638  assert (!deltaSL_);
639  deltaSL_ = new double [nTotal];
640  CoinZeroN(deltaSL_,nTotal);
641  assert (!primalR_);
642  assert (!dualR_);
643  // create arrays if we are doing KKT
644  if (cholesky_->type()>=20) {
645    primalR_ = new double [nTotal];
646    CoinZeroN(primalR_,nTotal);
647    dualR_ = new double [numberRows_];
648    CoinZeroN(dualR_,numberRows_);
649  }
650  assert (!rhsB_);
651  rhsB_ = new double [numberRows_];
652  CoinZeroN(rhsB_,numberRows_);
653  assert (!rhsU_);
654  rhsU_ = new double [nTotal];
655  CoinZeroN(rhsU_,nTotal);
656  assert (!rhsL_);
657  rhsL_ = new double [nTotal];
658  CoinZeroN(rhsL_,nTotal);
659  assert (!rhsZ_);
660  rhsZ_ = new double [nTotal];
661  CoinZeroN(rhsZ_,nTotal);
662  assert (!rhsW_);
663  rhsW_ = new double [nTotal];
664  CoinZeroN(rhsW_,nTotal);
665  assert (!rhsC_);
666  rhsC_ = new double [nTotal];
667  CoinZeroN(rhsC_,nTotal);
668  assert (!workArray_);
669  workArray_ = new double [nTotal];
670  CoinZeroN(workArray_,nTotal);
671  assert (!zVec_);
672  zVec_ = new double [nTotal];
673  CoinZeroN(zVec_,nTotal);
674  assert (!wVec_);
675  wVec_ = new double [nTotal];
676  CoinZeroN(wVec_,nTotal);
677  assert (!dj_);
678  dj_ = new double [nTotal];
679  if (!status_)
680    status_ = new unsigned char [numberRows_+numberColumns_];
681  memset(status_,0,numberRows_+numberColumns_);
682  return goodMatrix;
683}
684void
685ClpInterior::deleteWorkingData()
686{
687  int i;
688  if (optimizationDirection_!=1.0||objectiveScale_!=1.0) {
689    double scaleC = optimizationDirection_/objectiveScale_;
690    // and modify all dual signs
691    for (i=0;i<numberColumns_;i++) 
692      reducedCost_[i] = scaleC*dj_[i];
693    for (i=0;i<numberRows_;i++) 
694      dual_[i] *= scaleC;
695  }
696  if (rowScale_) {
697    double scaleR = 1.0/rhsScale_;
698    for (i=0;i<numberColumns_;i++) {
699      double scaleFactor = columnScale_[i];
700      double valueScaled = columnActivity_[i];
701      columnActivity_[i] = valueScaled*scaleFactor*scaleR;
702      double valueScaledDual = reducedCost_[i];
703      reducedCost_[i] = valueScaledDual/scaleFactor;
704    }
705    for (i=0;i<numberRows_;i++) {
706      double scaleFactor = rowScale_[i];
707      double valueScaled = rowActivity_[i];
708      rowActivity_[i] = (valueScaled*scaleR)/scaleFactor;
709      double valueScaledDual = dual_[i];
710      dual_[i] = valueScaledDual*scaleFactor;
711    }
712  } else if (rhsScale_!=1.0) {
713    double scaleR = 1.0/rhsScale_;
714    for (i=0;i<numberColumns_;i++) {
715      double valueScaled = columnActivity_[i];
716      columnActivity_[i] = valueScaled*scaleR;
717    }
718    for (i=0;i<numberRows_;i++) {
719      double valueScaled = rowActivity_[i];
720      rowActivity_[i] = valueScaled*scaleR;
721    }
722  }
723  delete [] cost_;
724  cost_ = NULL;
725  delete [] solution_;
726  solution_ = NULL;
727  delete [] lower_;
728  lower_ = NULL;
729  delete [] upper_;
730  upper_ = NULL;
731  delete [] errorRegion_;
732  errorRegion_ = NULL;
733  delete [] rhsFixRegion_;
734  rhsFixRegion_ = NULL;
735  delete [] deltaY_;
736  deltaY_ = NULL;
737  delete [] upperSlack_;
738  upperSlack_ = NULL;
739  delete [] lowerSlack_;
740  lowerSlack_ = NULL;
741  delete [] diagonal_;
742  diagonal_ = NULL;
743  delete [] deltaX_;
744  deltaX_ = NULL;
745  delete [] workArray_;
746  workArray_ = NULL;
747  delete [] zVec_;
748  zVec_ = NULL;
749  delete [] wVec_;
750  wVec_ = NULL;
751  delete [] dj_;
752  dj_ = NULL;
753}
754// Sanity check on input data - returns true if okay
755bool 
756ClpInterior::sanityCheck()
757{
758  // bad if empty
759  if (!numberColumns_||((!numberRows_||!matrix_->getNumElements())&&objective_->type()<2)) {
760    problemStatus_=emptyProblem();
761    return false;
762  }
763  int numberBad ;
764  double largestBound, smallestBound, minimumGap;
765  double smallestObj, largestObj;
766  int firstBad;
767  int modifiedBounds=0;
768  int i;
769  numberBad=0;
770  firstBad=-1;
771  minimumGap=1.0e100;
772  smallestBound=1.0e100;
773  largestBound=0.0;
774  smallestObj=1.0e100;
775  largestObj=0.0;
776  // If bounds are too close - fix
777  double fixTolerance = 1.1*primalTolerance();
778  for (i=numberColumns_;i<numberColumns_+numberRows_;i++) {
779    double value;
780    value = fabs(cost_[i]);
781    if (value>1.0e50) {
782      numberBad++;
783      if (firstBad<0)
784        firstBad=i;
785    } else if (value) {
786      if (value>largestObj)
787        largestObj=value;
788      if (value<smallestObj)
789        smallestObj=value;
790    }
791    value=upper_[i]-lower_[i];
792    if (value<-primalTolerance()) {
793      numberBad++;
794      if (firstBad<0)
795        firstBad=i;
796    } else if (value<=fixTolerance) {
797      if (value) {
798        // modify
799        upper_[i] = lower_[i];
800        modifiedBounds++;
801      }
802    } else {
803      if (value<minimumGap)
804        minimumGap=value;
805    }
806    if (lower_[i]>-1.0e100&&lower_[i]) {
807      value = fabs(lower_[i]);
808      if (value>largestBound)
809        largestBound=value;
810      if (value<smallestBound)
811        smallestBound=value;
812    }
813    if (upper_[i]<1.0e100&&upper_[i]) {
814      value = fabs(upper_[i]);
815      if (value>largestBound)
816        largestBound=value;
817      if (value<smallestBound)
818        smallestBound=value;
819    }
820  }
821  if (largestBound)
822    handler_->message(CLP_RIMSTATISTICS3,messages_)
823      <<smallestBound
824      <<largestBound
825      <<minimumGap
826      <<CoinMessageEol;
827  minimumGap=1.0e100;
828  smallestBound=1.0e100;
829  largestBound=0.0;
830  for (i=0;i<numberColumns_;i++) {
831    double value;
832    value = fabs(cost_[i]);
833    if (value>1.0e50) {
834      numberBad++;
835      if (firstBad<0)
836        firstBad=i;
837    } else if (value) {
838      if (value>largestObj)
839        largestObj=value;
840      if (value<smallestObj)
841        smallestObj=value;
842    }
843    value=upper_[i]-lower_[i];
844    if (value<-primalTolerance()) {
845      numberBad++;
846      if (firstBad<0)
847        firstBad=i;
848    } else if (value<=fixTolerance) {
849      if (value) {
850        // modify
851        upper_[i] = lower_[i];
852        modifiedBounds++;
853      }
854    } else {
855      if (value<minimumGap)
856        minimumGap=value;
857    }
858    if (lower_[i]>-1.0e100&&lower_[i]) {
859      value = fabs(lower_[i]);
860      if (value>largestBound)
861        largestBound=value;
862      if (value<smallestBound)
863        smallestBound=value;
864    }
865    if (upper_[i]<1.0e100&&upper_[i]) {
866      value = fabs(upper_[i]);
867      if (value>largestBound)
868        largestBound=value;
869      if (value<smallestBound)
870        smallestBound=value;
871    }
872  }
873  char rowcol[]={'R','C'};
874  if (numberBad) {
875    handler_->message(CLP_BAD_BOUNDS,messages_)
876      <<numberBad
877      <<rowcol[isColumn(firstBad)]<<sequenceWithin(firstBad)
878      <<CoinMessageEol;
879    problemStatus_=4;
880    return false;
881  }
882  if (modifiedBounds)
883    handler_->message(CLP_MODIFIEDBOUNDS,messages_)
884      <<modifiedBounds
885      <<CoinMessageEol;
886  handler_->message(CLP_RIMSTATISTICS1,messages_)
887    <<smallestObj
888    <<largestObj
889    <<CoinMessageEol;  if (largestBound)
890    handler_->message(CLP_RIMSTATISTICS2,messages_)
891      <<smallestBound
892      <<largestBound
893      <<minimumGap
894      <<CoinMessageEol;
895  return true;
896}
897/* Loads a problem (the constraints on the
898   rows are given by lower and upper bounds). If a pointer is 0 then the
899   following values are the default:
900   <ul>
901   <li> <code>colub</code>: all columns have upper bound infinity
902   <li> <code>collb</code>: all columns have lower bound 0
903   <li> <code>rowub</code>: all rows have upper bound infinity
904   <li> <code>rowlb</code>: all rows have lower bound -infinity
905   <li> <code>obj</code>: all variables have 0 objective coefficient
906   </ul>
907*/
908void 
909ClpInterior::loadProblem (  const ClpMatrixBase& matrix,
910                    const double* collb, const double* colub,   
911                    const double* obj,
912                    const double* rowlb, const double* rowub,
913                    const double * rowObjective)
914{
915  ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
916                        rowObjective);
917}
918void 
919ClpInterior::loadProblem (  const CoinPackedMatrix& matrix,
920                    const double* collb, const double* colub,   
921                    const double* obj,
922                    const double* rowlb, const double* rowub,
923                    const double * rowObjective)
924{
925  ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
926                        rowObjective);
927}
928
929/* Just like the other loadProblem() method except that the matrix is
930   given in a standard column major ordered format (without gaps). */
931void 
932ClpInterior::loadProblem (  const int numcols, const int numrows,
933                    const CoinBigIndex* start, const int* index,
934                    const double* value,
935                    const double* collb, const double* colub,   
936                    const double* obj,
937                    const double* rowlb, const double* rowub,
938                    const double * rowObjective)
939{
940  ClpModel::loadProblem(numcols, numrows, start, index, value,
941                          collb, colub, obj, rowlb, rowub,
942                          rowObjective);
943}
944void 
945ClpInterior::loadProblem (  const int numcols, const int numrows,
946                           const CoinBigIndex* start, const int* index,
947                           const double* value,const int * length,
948                           const double* collb, const double* colub,   
949                           const double* obj,
950                           const double* rowlb, const double* rowub,
951                           const double * rowObjective)
952{
953  ClpModel::loadProblem(numcols, numrows, start, index, value, length,
954                          collb, colub, obj, rowlb, rowub,
955                          rowObjective);
956}
957// Read an mps file from the given filename
958int 
959ClpInterior::readMps(const char *filename,
960            bool keepNames,
961            bool ignoreErrors)
962{
963  int status = ClpModel::readMps(filename,keepNames,ignoreErrors);
964  return status;
965}
966#ifdef COIN_DO_PDCO
967#include "ClpPdco.hpp"
968/* Pdco algorithm - see ClpPdco.hpp for method */
969int 
970ClpInterior::pdco()
971{
972  return ((ClpPdco *) this)->pdco();
973}
974// ** Temporary version
975int 
976ClpInterior::pdco( ClpPdcoBase * stuff, Options &options, Info &info, Outfo &outfo)
977{
978  return ((ClpPdco *) this)->pdco(stuff,options,info,outfo);
979}
980#endif
981#include "ClpPredictorCorrector.hpp"
982// Primal-Dual Predictor-Corrector barrier
983int 
984ClpInterior::primalDual()
985{ 
986  return (static_cast<ClpPredictorCorrector *> (this))->solve();
987}
988
989void 
990ClpInterior::checkSolution()
991{
992  int iRow,iColumn;
993  CoinMemcpyN(cost_,numberColumns_,reducedCost_);
994  matrix_->transposeTimes(-1.0,dual_,reducedCost_);
995  // Now modify reduced costs for quadratic
996  double quadraticOffset=quadraticDjs(reducedCost_,
997                                      solution_,scaleFactor_);
998
999  objectiveValue_ = 0.0;
1000  // now look at solution
1001  sumPrimalInfeasibilities_=0.0;
1002  sumDualInfeasibilities_=0.0;
1003  double dualTolerance =  10.0*dblParam_[ClpDualTolerance];
1004  double primalTolerance =  dblParam_[ClpPrimalTolerance];
1005  double primalTolerance2 =  10.0*dblParam_[ClpPrimalTolerance];
1006  worstComplementarity_=0.0;
1007  complementarityGap_=0.0;
1008
1009  // Done scaled - use permanent regions for output
1010  // but internal for bounds
1011  const double * lower = lower_+numberColumns_;
1012  const double * upper = upper_+numberColumns_;
1013  for (iRow=0;iRow<numberRows_;iRow++) {
1014    double infeasibility=0.0;
1015    double distanceUp = CoinMin(upper[iRow]-
1016      rowActivity_[iRow],1.0e10);
1017    double distanceDown = CoinMin(rowActivity_[iRow] -
1018      lower[iRow],1.0e10);
1019    if (distanceUp>primalTolerance2) {
1020      double value = dual_[iRow];
1021      // should not be negative
1022      if (value<-dualTolerance) {
1023        sumDualInfeasibilities_ += -dualTolerance-value;
1024        value = - value*distanceUp;
1025        if (value>worstComplementarity_) 
1026          worstComplementarity_=value;
1027        complementarityGap_ += value;
1028      }
1029    }
1030    if (distanceDown>primalTolerance2) {
1031      double value = dual_[iRow];
1032      // should not be positive
1033      if (value>dualTolerance) {
1034        sumDualInfeasibilities_ += value-dualTolerance;
1035        value =  value*distanceDown;
1036        if (value>worstComplementarity_) 
1037          worstComplementarity_=value;
1038        complementarityGap_ += value;
1039      }
1040    }
1041    if (rowActivity_[iRow]>upper[iRow]) {
1042      infeasibility=rowActivity_[iRow]-upper[iRow];
1043    } else if (rowActivity_[iRow]<lower[iRow]) {
1044      infeasibility=lower[iRow]-rowActivity_[iRow];
1045    }
1046    if (infeasibility>primalTolerance) {
1047      sumPrimalInfeasibilities_ += infeasibility-primalTolerance;
1048    }
1049  }
1050  lower = lower_;
1051  upper = upper_;
1052  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1053    double infeasibility=0.0;
1054    objectiveValue_ += cost_[iColumn]*columnActivity_[iColumn];
1055    double distanceUp = CoinMin(upper[iColumn]-
1056      columnActivity_[iColumn],1.0e10);
1057    double distanceDown = CoinMin(columnActivity_[iColumn] -
1058      lower[iColumn],1.0e10);
1059    if (distanceUp>primalTolerance2) {
1060      double value = reducedCost_[iColumn];
1061      // should not be negative
1062      if (value<-dualTolerance) {
1063        sumDualInfeasibilities_ += -dualTolerance-value;
1064        value = - value*distanceUp;
1065        if (value>worstComplementarity_) 
1066          worstComplementarity_=value;
1067        complementarityGap_ += value;
1068      }
1069    }
1070    if (distanceDown>primalTolerance2) {
1071      double value = reducedCost_[iColumn];
1072      // should not be positive
1073      if (value>dualTolerance) {
1074        sumDualInfeasibilities_ += value-dualTolerance;
1075        value =  value*distanceDown;
1076        if (value>worstComplementarity_) 
1077          worstComplementarity_=value;
1078        complementarityGap_ += value;
1079      }
1080    }
1081    if (columnActivity_[iColumn]>upper[iColumn]) {
1082      infeasibility=columnActivity_[iColumn]-upper[iColumn];
1083    } else if (columnActivity_[iColumn]<lower[iColumn]) {
1084      infeasibility=lower[iColumn]-columnActivity_[iColumn];
1085    }
1086    if (infeasibility>primalTolerance) {
1087      sumPrimalInfeasibilities_ += infeasibility-primalTolerance;
1088    }
1089  }
1090  // add in offset
1091  objectiveValue_ += 0.5*quadraticOffset;
1092}
1093// Set cholesky (and delete present one)
1094void 
1095ClpInterior::setCholesky(ClpCholeskyBase * cholesky)
1096{
1097  delete cholesky_;
1098  cholesky_= cholesky;
1099}
1100/* Borrow model.  This is so we dont have to copy large amounts
1101   of data around.  It assumes a derived class wants to overwrite
1102   an empty model with a real one - while it does an algorithm.
1103   This is same as ClpModel one. */
1104void 
1105ClpInterior::borrowModel(ClpModel & otherModel) 
1106{
1107  ClpModel::borrowModel(otherModel);
1108}
1109/* Return model - updates any scalars */
1110void 
1111ClpInterior::returnModel(ClpModel & otherModel)
1112{
1113  ClpModel::returnModel(otherModel);
1114}
1115// Return number fixed to see if worth presolving
1116int 
1117ClpInterior::numberFixed() const
1118{
1119  int i;
1120  int nFixed=0;
1121  for (i=0;i<numberColumns_;i++) {
1122    if (columnUpper_[i]<1.0e20||columnLower_[i]>-1.0e20) { 
1123      if (columnUpper_[i]>columnLower_[i]) { 
1124        if (fixedOrFree(i))
1125          nFixed++;
1126      }
1127    }
1128  }
1129  for (i=0;i<numberRows_;i++) {
1130    if (rowUpper_[i]<1.0e20||rowLower_[i]>-1.0e20) { 
1131      if (rowUpper_[i]>rowLower_[i]) { 
1132        if (fixedOrFree(i+numberColumns_))
1133          nFixed++;
1134      }
1135    }
1136  }
1137  return nFixed;
1138}
1139// fix variables interior says should be
1140void 
1141ClpInterior::fixFixed(bool reallyFix)
1142{
1143  // Arrays for change in columns and rhs
1144  double * columnChange = new double[numberColumns_];
1145  double * rowChange = new double[numberRows_];
1146  CoinZeroN(columnChange,numberColumns_);
1147  CoinZeroN(rowChange,numberRows_);
1148  matrix_->times(1.0,columnChange,rowChange);
1149  int i;
1150  double tolerance = primalTolerance();
1151  for (i=0;i<numberColumns_;i++) {
1152    if (columnUpper_[i]<1.0e20||columnLower_[i]>-1.0e20) { 
1153      if (columnUpper_[i]>columnLower_[i]) { 
1154        if (fixedOrFree(i)) {
1155          if (columnActivity_[i]-columnLower_[i]<columnUpper_[i]-columnActivity_[i]) {
1156            double change = columnLower_[i]-columnActivity_[i]; 
1157            if (fabs(change)<tolerance) {
1158              if (reallyFix)
1159                columnUpper_[i]=columnLower_[i];
1160              columnChange[i] = change;
1161              columnActivity_[i]=columnLower_[i];
1162            }
1163          } else {
1164            double change = columnUpper_[i]-columnActivity_[i]; 
1165            if (fabs(change)<tolerance) {
1166              if (reallyFix)
1167                columnLower_[i]=columnUpper_[i];
1168              columnChange[i] = change;
1169              columnActivity_[i]=columnUpper_[i];
1170            }
1171          }
1172        }
1173      }
1174    }
1175  }
1176  CoinZeroN(rowChange,numberRows_);
1177  matrix_->times(1.0,columnChange,rowChange);
1178  // If makes mess of things then don't do
1179  double newSum=0.0;
1180  for (i=0;i<numberRows_;i++) {
1181    double value=rowActivity_[i]+rowChange[i];
1182    if (value>rowUpper_[i]+tolerance)
1183      newSum += value-rowUpper_[i]-tolerance;
1184    else if (value<rowLower_[i]-tolerance)
1185      newSum -= value-rowLower_[i]+tolerance;
1186  }
1187  if (newSum>1.0e-5+1.5*sumPrimalInfeasibilities_) {
1188    // put back and skip changes
1189    for (i=0;i<numberColumns_;i++) 
1190      columnActivity_[i] -= columnChange[i];
1191  } else {
1192    CoinZeroN(rowActivity_,numberRows_);
1193    matrix_->times(1.0,columnActivity_,rowActivity_);
1194    if (reallyFix) {
1195      for (i=0;i<numberRows_;i++) {
1196        if (rowUpper_[i]<1.0e20||rowLower_[i]>-1.0e20) { 
1197          if (rowUpper_[i]>rowLower_[i]) { 
1198            if (fixedOrFree(i+numberColumns_)) {
1199              if (rowActivity_[i]-rowLower_[i]<rowUpper_[i]-rowActivity_[i]) {
1200                double change = rowLower_[i]-rowActivity_[i]; 
1201                if (fabs(change)<tolerance) {
1202                  if (reallyFix)
1203                    rowUpper_[i]=rowLower_[i];
1204                  rowActivity_[i]=rowLower_[i];
1205                }
1206              } else {
1207                double change = rowLower_[i]-rowActivity_[i]; 
1208                if (fabs(change)<tolerance) {
1209                  if (reallyFix)
1210                    rowLower_[i]=rowUpper_[i];
1211                  rowActivity_[i]=rowUpper_[i];
1212                }
1213              }
1214            }
1215          }
1216        }
1217      }
1218    }
1219  }
1220  delete [] rowChange;
1221  delete [] columnChange;
1222}
1223/* Modifies djs to allow for quadratic.
1224   returns quadratic offset */
1225double 
1226ClpInterior::quadraticDjs(double * djRegion, const double * solution,
1227                          double scaleFactor)
1228{
1229  double quadraticOffset=0.0;
1230#ifndef NO_RTTI
1231  ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
1232#else
1233  ClpQuadraticObjective * quadraticObj = NULL;
1234  if (objective_->type()==2)
1235    quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
1236#endif
1237  if (quadraticObj) {
1238    CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
1239    const int * columnQuadratic = quadratic->getIndices();
1240    const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
1241    const int * columnQuadraticLength = quadratic->getVectorLengths();
1242    double * quadraticElement = quadratic->getMutableElements();
1243    int numberColumns = quadratic->getNumCols();
1244    for (int iColumn=0;iColumn<numberColumns;iColumn++) {
1245      double value = 0.0;
1246      for (CoinBigIndex j=columnQuadraticStart[iColumn];
1247           j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
1248        int jColumn = columnQuadratic[j];
1249        double valueJ = solution[jColumn];
1250        double elementValue = quadraticElement[j];
1251        //value += valueI*valueJ*elementValue;
1252        value += valueJ*elementValue;
1253        quadraticOffset += solution[iColumn]*valueJ*elementValue;
1254      }
1255      djRegion[iColumn] += scaleFactor*value;
1256    }
1257  }
1258  return quadraticOffset;
1259}
Note: See TracBrowser for help on using the repository browser.