source: tags/move-to-subversion/ClpInterior.cpp @ 753

Last change on this file since 753 was 739, checked in by forrest, 14 years ago

trying to improve

  • 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 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 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 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  memcpy(solution_,columnActivity_,
541         numberColumns_*sizeof(double));
542  memcpy(solution_+numberColumns_,rowActivity_,
543         numberRows_*sizeof(double));
544  delete [] cost_;
545  cost_ = new double[nTotal];
546  int i;
547  double direction = optimizationDirection_*objectiveScale_;
548  // direction is actually scale out not scale in
549  if (direction)
550    direction = 1.0/direction;
551  const double * obj = objective();
552  for (i=0;i<numberColumns_;i++)
553    cost_[i] = direction*obj[i];
554  memset(cost_+numberColumns_,0,numberRows_*sizeof(double));
555  // do scaling if needed
556  if (scalingFlag_>0&&!rowScale_) {
557    if (matrix_->scale(this))
558      scalingFlag_=-scalingFlag_; // not scaled after all
559  }
560  delete [] lower_;
561  delete [] upper_;
562  lower_ = new double[nTotal];
563  upper_ = new double[nTotal];
564  rowLowerWork_ = lower_+numberColumns_;
565  columnLowerWork_ = lower_;
566  rowUpperWork_ = upper_+numberColumns_;
567  columnUpperWork_ = upper_;
568  memcpy(rowLowerWork_,rowLower_,numberRows_*sizeof(double));
569  memcpy(rowUpperWork_,rowUpper_,numberRows_*sizeof(double));
570  memcpy(columnLowerWork_,columnLower_,numberColumns_*sizeof(double));
571  memcpy(columnUpperWork_,columnUpper_,numberColumns_*sizeof(double));
572  // clean up any mismatches on infinity
573  for (i=0;i<numberColumns_;i++) {
574    if (columnLowerWork_[i]<-1.0e30)
575      columnLowerWork_[i] = -COIN_DBL_MAX;
576    if (columnUpperWork_[i]>1.0e30)
577      columnUpperWork_[i] = COIN_DBL_MAX;
578  }
579  // clean up any mismatches on infinity
580  for (i=0;i<numberRows_;i++) {
581    if (rowLowerWork_[i]<-1.0e30)
582      rowLowerWork_[i] = -COIN_DBL_MAX;
583    if (rowUpperWork_[i]>1.0e30)
584      rowUpperWork_[i] = COIN_DBL_MAX;
585  }
586  // check rim of problem okay
587  if (!sanityCheck())
588    goodMatrix=false;
589  if(rowScale_) {
590    for (i=0;i<numberColumns_;i++) {
591      double multiplier = rhsScale_/columnScale_[i];
592      cost_[i] *= columnScale_[i];
593      if (columnLowerWork_[i]>-1.0e50)
594        columnLowerWork_[i] *= multiplier;
595      if (columnUpperWork_[i]<1.0e50)
596        columnUpperWork_[i] *= multiplier;
597     
598    }
599    for (i=0;i<numberRows_;i++) {
600      double multiplier = rhsScale_*rowScale_[i];
601      if (rowLowerWork_[i]>-1.0e50)
602        rowLowerWork_[i] *= multiplier;
603      if (rowUpperWork_[i]<1.0e50)
604        rowUpperWork_[i] *= multiplier;
605    }
606  } else if (rhsScale_!=1.0) {
607    for (i=0;i<numberColumns_+numberRows_;i++) {
608      if (lower_[i]>-1.0e50)
609        lower_[i] *= rhsScale_;
610      if (upper_[i]<1.0e50)
611        upper_[i] *= rhsScale_;
612     
613    }
614  }
615  assert (!errorRegion_);
616  errorRegion_ = new double [numberRows_];
617  assert (!rhsFixRegion_);
618  rhsFixRegion_ = new double [numberRows_];
619  assert (!deltaY_);
620  deltaY_ = new double [numberRows_];
621  CoinZeroN(deltaY_,numberRows_);
622  assert (!upperSlack_);
623  upperSlack_ = new double [nTotal];
624  assert (!lowerSlack_);
625  lowerSlack_ = new double [nTotal];
626  assert (!diagonal_);
627  diagonal_ = new double [nTotal];
628  assert (!deltaX_);
629  deltaX_ = new double [nTotal];
630  CoinZeroN(deltaX_,nTotal);
631  assert (!deltaZ_);
632  deltaZ_ = new double [nTotal];
633  CoinZeroN(deltaZ_,nTotal);
634  assert (!deltaW_);
635  deltaW_ = new double [nTotal];
636  CoinZeroN(deltaW_,nTotal);
637  assert (!deltaSU_);
638  deltaSU_ = new double [nTotal];
639  CoinZeroN(deltaSU_,nTotal);
640  assert (!deltaSL_);
641  deltaSL_ = new double [nTotal];
642  CoinZeroN(deltaSL_,nTotal);
643  assert (!primalR_);
644  assert (!dualR_);
645  // create arrays if we are doing KKT
646  if (cholesky_->type()>=20) {
647    primalR_ = new double [nTotal];
648    CoinZeroN(primalR_,nTotal);
649    dualR_ = new double [numberRows_];
650    CoinZeroN(dualR_,numberRows_);
651  }
652  assert (!rhsB_);
653  rhsB_ = new double [numberRows_];
654  CoinZeroN(rhsB_,numberRows_);
655  assert (!rhsU_);
656  rhsU_ = new double [nTotal];
657  CoinZeroN(rhsU_,nTotal);
658  assert (!rhsL_);
659  rhsL_ = new double [nTotal];
660  CoinZeroN(rhsL_,nTotal);
661  assert (!rhsZ_);
662  rhsZ_ = new double [nTotal];
663  CoinZeroN(rhsZ_,nTotal);
664  assert (!rhsW_);
665  rhsW_ = new double [nTotal];
666  CoinZeroN(rhsW_,nTotal);
667  assert (!rhsC_);
668  rhsC_ = new double [nTotal];
669  CoinZeroN(rhsC_,nTotal);
670  assert (!workArray_);
671  workArray_ = new double [nTotal];
672  CoinZeroN(workArray_,nTotal);
673  assert (!zVec_);
674  zVec_ = new double [nTotal];
675  CoinZeroN(zVec_,nTotal);
676  assert (!wVec_);
677  wVec_ = new double [nTotal];
678  CoinZeroN(wVec_,nTotal);
679  assert (!dj_);
680  dj_ = new double [nTotal];
681  if (!status_)
682    status_ = new unsigned char [numberRows_+numberColumns_];
683  memset(status_,0,numberRows_+numberColumns_);
684  return goodMatrix;
685}
686void
687ClpInterior::deleteWorkingData()
688{
689  int i;
690  if (optimizationDirection_!=1.0||objectiveScale_!=1.0) {
691    double scaleC = optimizationDirection_/objectiveScale_;
692    // and modify all dual signs
693    for (i=0;i<numberColumns_;i++) 
694      reducedCost_[i] = scaleC*dj_[i];
695    for (i=0;i<numberRows_;i++) 
696      dual_[i] *= scaleC;
697  }
698  if (rowScale_) {
699    double scaleR = 1.0/rhsScale_;
700    for (i=0;i<numberColumns_;i++) {
701      double scaleFactor = columnScale_[i];
702      double valueScaled = columnActivity_[i];
703      columnActivity_[i] = valueScaled*scaleFactor*scaleR;
704      double valueScaledDual = reducedCost_[i];
705      reducedCost_[i] = valueScaledDual/scaleFactor;
706    }
707    for (i=0;i<numberRows_;i++) {
708      double scaleFactor = rowScale_[i];
709      double valueScaled = rowActivity_[i];
710      rowActivity_[i] = (valueScaled*scaleR)/scaleFactor;
711      double valueScaledDual = dual_[i];
712      dual_[i] = valueScaledDual*scaleFactor;
713    }
714  } else if (rhsScale_!=1.0) {
715    double scaleR = 1.0/rhsScale_;
716    for (i=0;i<numberColumns_;i++) {
717      double valueScaled = columnActivity_[i];
718      columnActivity_[i] = valueScaled*scaleR;
719    }
720    for (i=0;i<numberRows_;i++) {
721      double valueScaled = rowActivity_[i];
722      rowActivity_[i] = valueScaled*scaleR;
723    }
724  }
725  delete [] cost_;
726  cost_ = NULL;
727  delete [] solution_;
728  solution_ = NULL;
729  delete [] lower_;
730  lower_ = NULL;
731  delete [] upper_;
732  upper_ = NULL;
733  delete [] errorRegion_;
734  errorRegion_ = NULL;
735  delete [] rhsFixRegion_;
736  rhsFixRegion_ = NULL;
737  delete [] deltaY_;
738  deltaY_ = NULL;
739  delete [] upperSlack_;
740  upperSlack_ = NULL;
741  delete [] lowerSlack_;
742  lowerSlack_ = NULL;
743  delete [] diagonal_;
744  diagonal_ = NULL;
745  delete [] deltaX_;
746  deltaX_ = NULL;
747  delete [] workArray_;
748  workArray_ = NULL;
749  delete [] zVec_;
750  zVec_ = NULL;
751  delete [] wVec_;
752  wVec_ = NULL;
753  delete [] dj_;
754  dj_ = NULL;
755}
756// Sanity check on input data - returns true if okay
757bool 
758ClpInterior::sanityCheck()
759{
760  // bad if empty
761  if (!numberColumns_||((!numberRows_||!matrix_->getNumElements())&&objective_->type()<2)) {
762    problemStatus_=emptyProblem();
763    return false;
764  }
765  int numberBad ;
766  double largestBound, smallestBound, minimumGap;
767  double smallestObj, largestObj;
768  int firstBad;
769  int modifiedBounds=0;
770  int i;
771  numberBad=0;
772  firstBad=-1;
773  minimumGap=1.0e100;
774  smallestBound=1.0e100;
775  largestBound=0.0;
776  smallestObj=1.0e100;
777  largestObj=0.0;
778  // If bounds are too close - fix
779  double fixTolerance = 1.1*primalTolerance();
780  for (i=numberColumns_;i<numberColumns_+numberRows_;i++) {
781    double value;
782    value = fabs(cost_[i]);
783    if (value>1.0e50) {
784      numberBad++;
785      if (firstBad<0)
786        firstBad=i;
787    } else if (value) {
788      if (value>largestObj)
789        largestObj=value;
790      if (value<smallestObj)
791        smallestObj=value;
792    }
793    value=upper_[i]-lower_[i];
794    if (value<-primalTolerance()) {
795      numberBad++;
796      if (firstBad<0)
797        firstBad=i;
798    } else if (value<=fixTolerance) {
799      if (value) {
800        // modify
801        upper_[i] = lower_[i];
802        modifiedBounds++;
803      }
804    } else {
805      if (value<minimumGap)
806        minimumGap=value;
807    }
808    if (lower_[i]>-1.0e100&&lower_[i]) {
809      value = fabs(lower_[i]);
810      if (value>largestBound)
811        largestBound=value;
812      if (value<smallestBound)
813        smallestBound=value;
814    }
815    if (upper_[i]<1.0e100&&upper_[i]) {
816      value = fabs(upper_[i]);
817      if (value>largestBound)
818        largestBound=value;
819      if (value<smallestBound)
820        smallestBound=value;
821    }
822  }
823  if (largestBound)
824    handler_->message(CLP_RIMSTATISTICS3,messages_)
825      <<smallestBound
826      <<largestBound
827      <<minimumGap
828      <<CoinMessageEol;
829  minimumGap=1.0e100;
830  smallestBound=1.0e100;
831  largestBound=0.0;
832  for (i=0;i<numberColumns_;i++) {
833    double value;
834    value = fabs(cost_[i]);
835    if (value>1.0e50) {
836      numberBad++;
837      if (firstBad<0)
838        firstBad=i;
839    } else if (value) {
840      if (value>largestObj)
841        largestObj=value;
842      if (value<smallestObj)
843        smallestObj=value;
844    }
845    value=upper_[i]-lower_[i];
846    if (value<-primalTolerance()) {
847      numberBad++;
848      if (firstBad<0)
849        firstBad=i;
850    } else if (value<=fixTolerance) {
851      if (value) {
852        // modify
853        upper_[i] = lower_[i];
854        modifiedBounds++;
855      }
856    } else {
857      if (value<minimumGap)
858        minimumGap=value;
859    }
860    if (lower_[i]>-1.0e100&&lower_[i]) {
861      value = fabs(lower_[i]);
862      if (value>largestBound)
863        largestBound=value;
864      if (value<smallestBound)
865        smallestBound=value;
866    }
867    if (upper_[i]<1.0e100&&upper_[i]) {
868      value = fabs(upper_[i]);
869      if (value>largestBound)
870        largestBound=value;
871      if (value<smallestBound)
872        smallestBound=value;
873    }
874  }
875  char rowcol[]={'R','C'};
876  if (numberBad) {
877    handler_->message(CLP_BAD_BOUNDS,messages_)
878      <<numberBad
879      <<rowcol[isColumn(firstBad)]<<sequenceWithin(firstBad)
880      <<CoinMessageEol;
881    problemStatus_=4;
882    return false;
883  }
884  if (modifiedBounds)
885    handler_->message(CLP_MODIFIEDBOUNDS,messages_)
886      <<modifiedBounds
887      <<CoinMessageEol;
888  handler_->message(CLP_RIMSTATISTICS1,messages_)
889    <<smallestObj
890    <<largestObj
891    <<CoinMessageEol;  if (largestBound)
892    handler_->message(CLP_RIMSTATISTICS2,messages_)
893      <<smallestBound
894      <<largestBound
895      <<minimumGap
896      <<CoinMessageEol;
897  return true;
898}
899/* Loads a problem (the constraints on the
900   rows are given by lower and upper bounds). If a pointer is 0 then the
901   following values are the default:
902   <ul>
903   <li> <code>colub</code>: all columns have upper bound infinity
904   <li> <code>collb</code>: all columns have lower bound 0
905   <li> <code>rowub</code>: all rows have upper bound infinity
906   <li> <code>rowlb</code>: all rows have lower bound -infinity
907   <li> <code>obj</code>: all variables have 0 objective coefficient
908   </ul>
909*/
910void 
911ClpInterior::loadProblem (  const ClpMatrixBase& matrix,
912                    const double* collb, const double* colub,   
913                    const double* obj,
914                    const double* rowlb, const double* rowub,
915                    const double * rowObjective)
916{
917  ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
918                        rowObjective);
919}
920void 
921ClpInterior::loadProblem (  const CoinPackedMatrix& matrix,
922                    const double* collb, const double* colub,   
923                    const double* obj,
924                    const double* rowlb, const double* rowub,
925                    const double * rowObjective)
926{
927  ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
928                        rowObjective);
929}
930
931/* Just like the other loadProblem() method except that the matrix is
932   given in a standard column major ordered format (without gaps). */
933void 
934ClpInterior::loadProblem (  const int numcols, const int numrows,
935                    const CoinBigIndex* start, const int* index,
936                    const double* value,
937                    const double* collb, const double* colub,   
938                    const double* obj,
939                    const double* rowlb, const double* rowub,
940                    const double * rowObjective)
941{
942  ClpModel::loadProblem(numcols, numrows, start, index, value,
943                          collb, colub, obj, rowlb, rowub,
944                          rowObjective);
945}
946void 
947ClpInterior::loadProblem (  const int numcols, const int numrows,
948                           const CoinBigIndex* start, const int* index,
949                           const double* value,const int * length,
950                           const double* collb, const double* colub,   
951                           const double* obj,
952                           const double* rowlb, const double* rowub,
953                           const double * rowObjective)
954{
955  ClpModel::loadProblem(numcols, numrows, start, index, value, length,
956                          collb, colub, obj, rowlb, rowub,
957                          rowObjective);
958}
959// Read an mps file from the given filename
960int 
961ClpInterior::readMps(const char *filename,
962            bool keepNames,
963            bool ignoreErrors)
964{
965  int status = ClpModel::readMps(filename,keepNames,ignoreErrors);
966  return status;
967}
968#ifdef PDCO
969#include "ClpPdco.hpp"
970/* Pdco algorithm - see ClpPdco.hpp for method */
971int 
972ClpInterior::pdco()
973{
974  return ((ClpPdco *) this)->pdco();
975}
976// ** Temporary version
977int 
978ClpInterior::pdco( ClpPdcoBase * stuff, Options &options, Info &info, Outfo &outfo)
979{
980  return ((ClpPdco *) this)->pdco(stuff,options,info,outfo);
981}
982#endif
983#include "ClpPredictorCorrector.hpp"
984// Primal-Dual Predictor-Corrector barrier
985int 
986ClpInterior::primalDual()
987{ 
988  return ((ClpPredictorCorrector *) this)->solve();
989}
990
991void 
992ClpInterior::checkSolution()
993{
994  int iRow,iColumn;
995  CoinMemcpyN(cost_,numberColumns_,reducedCost_);
996  matrix_->transposeTimes(-1.0,dual_,reducedCost_);
997  // Now modify reduced costs for quadratic
998  double quadraticOffset=quadraticDjs(reducedCost_,
999                                      solution_,scaleFactor_);
1000
1001  objectiveValue_ = 0.0;
1002  // now look at solution
1003  sumPrimalInfeasibilities_=0.0;
1004  sumDualInfeasibilities_=0.0;
1005  double dualTolerance =  10.0*dblParam_[ClpDualTolerance];
1006  double primalTolerance =  dblParam_[ClpPrimalTolerance];
1007  double primalTolerance2 =  10.0*dblParam_[ClpPrimalTolerance];
1008  worstComplementarity_=0.0;
1009  complementarityGap_=0.0;
1010
1011  // Done scaled - use permanent regions for output
1012  // but internal for bounds
1013  const double * lower = lower_+numberColumns_;
1014  const double * upper = upper_+numberColumns_;
1015  for (iRow=0;iRow<numberRows_;iRow++) {
1016    double infeasibility=0.0;
1017    double distanceUp = CoinMin(upper[iRow]-
1018      rowActivity_[iRow],1.0e10);
1019    double distanceDown = CoinMin(rowActivity_[iRow] -
1020      lower[iRow],1.0e10);
1021    if (distanceUp>primalTolerance2) {
1022      double value = dual_[iRow];
1023      // should not be negative
1024      if (value<-dualTolerance) {
1025        sumDualInfeasibilities_ += -dualTolerance-value;
1026        value = - value*distanceUp;
1027        if (value>worstComplementarity_) 
1028          worstComplementarity_=value;
1029        complementarityGap_ += value;
1030      }
1031    }
1032    if (distanceDown>primalTolerance2) {
1033      double value = dual_[iRow];
1034      // should not be positive
1035      if (value>dualTolerance) {
1036        sumDualInfeasibilities_ += value-dualTolerance;
1037        value =  value*distanceDown;
1038        if (value>worstComplementarity_) 
1039          worstComplementarity_=value;
1040        complementarityGap_ += value;
1041      }
1042    }
1043    if (rowActivity_[iRow]>upper[iRow]) {
1044      infeasibility=rowActivity_[iRow]-upper[iRow];
1045    } else if (rowActivity_[iRow]<lower[iRow]) {
1046      infeasibility=lower[iRow]-rowActivity_[iRow];
1047    }
1048    if (infeasibility>primalTolerance) {
1049      sumPrimalInfeasibilities_ += infeasibility-primalTolerance;
1050    }
1051  }
1052  lower = lower_;
1053  upper = upper_;
1054  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1055    double infeasibility=0.0;
1056    objectiveValue_ += cost_[iColumn]*columnActivity_[iColumn];
1057    double distanceUp = CoinMin(upper[iColumn]-
1058      columnActivity_[iColumn],1.0e10);
1059    double distanceDown = CoinMin(columnActivity_[iColumn] -
1060      lower[iColumn],1.0e10);
1061    if (distanceUp>primalTolerance2) {
1062      double value = reducedCost_[iColumn];
1063      // should not be negative
1064      if (value<-dualTolerance) {
1065        sumDualInfeasibilities_ += -dualTolerance-value;
1066        value = - value*distanceUp;
1067        if (value>worstComplementarity_) 
1068          worstComplementarity_=value;
1069        complementarityGap_ += value;
1070      }
1071    }
1072    if (distanceDown>primalTolerance2) {
1073      double value = reducedCost_[iColumn];
1074      // should not be positive
1075      if (value>dualTolerance) {
1076        sumDualInfeasibilities_ += value-dualTolerance;
1077        value =  value*distanceDown;
1078        if (value>worstComplementarity_) 
1079          worstComplementarity_=value;
1080        complementarityGap_ += value;
1081      }
1082    }
1083    if (columnActivity_[iColumn]>upper[iColumn]) {
1084      infeasibility=columnActivity_[iColumn]-upper[iColumn];
1085    } else if (columnActivity_[iColumn]<lower[iColumn]) {
1086      infeasibility=lower[iColumn]-columnActivity_[iColumn];
1087    }
1088    if (infeasibility>primalTolerance) {
1089      sumPrimalInfeasibilities_ += infeasibility-primalTolerance;
1090    }
1091  }
1092  // add in offset
1093  objectiveValue_ += 0.5*quadraticOffset;
1094}
1095// Set cholesky (and delete present one)
1096void 
1097ClpInterior::setCholesky(ClpCholeskyBase * cholesky)
1098{
1099  delete cholesky_;
1100  cholesky_= cholesky;
1101}
1102/* Borrow model.  This is so we dont have to copy large amounts
1103   of data around.  It assumes a derived class wants to overwrite
1104   an empty model with a real one - while it does an algorithm.
1105   This is same as ClpModel one. */
1106void 
1107ClpInterior::borrowModel(ClpModel & otherModel) 
1108{
1109  ClpModel::borrowModel(otherModel);
1110}
1111/* Return model - updates any scalars */
1112void 
1113ClpInterior::returnModel(ClpModel & otherModel)
1114{
1115  ClpModel::returnModel(otherModel);
1116}
1117// Return number fixed to see if worth presolving
1118int 
1119ClpInterior::numberFixed() const
1120{
1121  int i;
1122  int nFixed=0;
1123  for (i=0;i<numberColumns_;i++) {
1124    if (columnUpper_[i]<1.0e20||columnLower_[i]>-1.0e20) { 
1125      if (columnUpper_[i]>columnLower_[i]) { 
1126        if (fixedOrFree(i))
1127          nFixed++;
1128      }
1129    }
1130  }
1131  for (i=0;i<numberRows_;i++) {
1132    if (rowUpper_[i]<1.0e20||rowLower_[i]>-1.0e20) { 
1133      if (rowUpper_[i]>rowLower_[i]) { 
1134        if (fixedOrFree(i+numberColumns_))
1135          nFixed++;
1136      }
1137    }
1138  }
1139  return nFixed;
1140}
1141// fix variables interior says should be
1142void 
1143ClpInterior::fixFixed(bool reallyFix)
1144{
1145  // Arrays for change in columns and rhs
1146  double * columnChange = new double[numberColumns_];
1147  double * rowChange = new double[numberRows_];
1148  CoinZeroN(columnChange,numberColumns_);
1149  CoinZeroN(rowChange,numberRows_);
1150  matrix_->times(1.0,columnChange,rowChange);
1151  int i;
1152  double tolerance = primalTolerance();
1153  for (i=0;i<numberColumns_;i++) {
1154    if (columnUpper_[i]<1.0e20||columnLower_[i]>-1.0e20) { 
1155      if (columnUpper_[i]>columnLower_[i]) { 
1156        if (fixedOrFree(i)) {
1157          if (columnActivity_[i]-columnLower_[i]<columnUpper_[i]-columnActivity_[i]) {
1158            double change = columnLower_[i]-columnActivity_[i]; 
1159            if (fabs(change)<tolerance) {
1160              if (reallyFix)
1161                columnUpper_[i]=columnLower_[i];
1162              columnChange[i] = change;
1163              columnActivity_[i]=columnLower_[i];
1164            }
1165          } else {
1166            double change = columnUpper_[i]-columnActivity_[i]; 
1167            if (fabs(change)<tolerance) {
1168              if (reallyFix)
1169                columnLower_[i]=columnUpper_[i];
1170              columnChange[i] = change;
1171              columnActivity_[i]=columnUpper_[i];
1172            }
1173          }
1174        }
1175      }
1176    }
1177  }
1178  CoinZeroN(rowChange,numberRows_);
1179  matrix_->times(1.0,columnChange,rowChange);
1180  // If makes mess of things then don't do
1181  double newSum=0.0;
1182  for (i=0;i<numberRows_;i++) {
1183    double value=rowActivity_[i]+rowChange[i];
1184    if (value>rowUpper_[i]+tolerance)
1185      newSum += value-rowUpper_[i]-tolerance;
1186    else if (value<rowLower_[i]-tolerance)
1187      newSum -= value-rowLower_[i]+tolerance;
1188  }
1189  if (newSum>1.0e-5+1.5*sumPrimalInfeasibilities_) {
1190    // put back and skip changes
1191    for (i=0;i<numberColumns_;i++) 
1192      columnActivity_[i] -= columnChange[i];
1193  } else {
1194    CoinZeroN(rowActivity_,numberRows_);
1195    matrix_->times(1.0,columnActivity_,rowActivity_);
1196    if (reallyFix) {
1197      for (i=0;i<numberRows_;i++) {
1198        if (rowUpper_[i]<1.0e20||rowLower_[i]>-1.0e20) { 
1199          if (rowUpper_[i]>rowLower_[i]) { 
1200            if (fixedOrFree(i+numberColumns_)) {
1201              if (rowActivity_[i]-rowLower_[i]<rowUpper_[i]-rowActivity_[i]) {
1202                double change = rowLower_[i]-rowActivity_[i]; 
1203                if (fabs(change)<tolerance) {
1204                  if (reallyFix)
1205                    rowUpper_[i]=rowLower_[i];
1206                  rowActivity_[i]=rowLower_[i];
1207                }
1208              } else {
1209                double change = rowLower_[i]-rowActivity_[i]; 
1210                if (fabs(change)<tolerance) {
1211                  if (reallyFix)
1212                    rowLower_[i]=rowUpper_[i];
1213                  rowActivity_[i]=rowUpper_[i];
1214                }
1215              }
1216            }
1217          }
1218        }
1219      }
1220    }
1221  }
1222  delete [] rowChange;
1223  delete [] columnChange;
1224}
1225/* Modifies djs to allow for quadratic.
1226   returns quadratic offset */
1227double 
1228ClpInterior::quadraticDjs(double * djRegion, const double * solution,
1229                          double scaleFactor)
1230{
1231  double quadraticOffset=0.0;
1232#ifndef NO_RTTI
1233  ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
1234#else
1235  ClpQuadraticObjective * quadraticObj = NULL;
1236  if (objective_->type()==2)
1237    quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
1238#endif
1239  if (quadraticObj) {
1240    CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
1241    const int * columnQuadratic = quadratic->getIndices();
1242    const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
1243    const int * columnQuadraticLength = quadratic->getVectorLengths();
1244    double * quadraticElement = quadratic->getMutableElements();
1245    int numberColumns = quadratic->getNumCols();
1246    for (int iColumn=0;iColumn<numberColumns;iColumn++) {
1247      double value = 0.0;
1248      for (CoinBigIndex j=columnQuadraticStart[iColumn];
1249           j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
1250        int jColumn = columnQuadratic[j];
1251        double valueJ = solution[jColumn];
1252        double elementValue = quadraticElement[j];
1253        //value += valueI*valueJ*elementValue;
1254        value += valueJ*elementValue;
1255        quadraticOffset += solution[iColumn]*valueJ*elementValue;
1256      }
1257      djRegion[iColumn] += scaleFactor*value;
1258    }
1259  }
1260  return quadraticOffset;
1261}
Note: See TracBrowser for help on using the repository browser.