Changeset 1367 for trunk/Clp


Ignore:
Timestamp:
May 23, 2009 3:47:00 AM (11 years ago)
Author:
forrest
Message:

try and improve stability - also long double interior point

Location:
trunk/Clp/src
Files:
22 edited

Legend:

Unmodified
Added
Removed
  • trunk/Clp/src/CbcOrClpParam.cpp

    r1365 r1367  
    28682868     );
    28692869  parameters[numberParameters-1].setIntValue(1);
     2870#ifdef CBC_KEEP_DEPRECATED
    28702871  parameters[numberParameters++]=
    28712872    CbcOrClpParam("strengthen","Create strengthened problem",
     
    28752876     "This creates a new problem by applying the root node cuts.  All tight constraints \
    28762877will be in resulting problem"
    2877      );
     2878     );
     2879#endif
    28782880  parameters[numberParameters++]=
    28792881    CbcOrClpParam("strong!Branching","Number of variables to look at in strong branching",
  • trunk/Clp/src/ClpCholeskyBase.cpp

    r1321 r1367  
    8686#else
    8787  // actually long double
    88   workDouble_ = (double *) ClpCopyOfArray((longWork *) rhs.workDouble_,numberRows_);
     88  workDouble_ = reinterpret_cast<double *> (ClpCopyOfArray(reinterpret_cast<CoinWorkDouble *> (rhs.workDouble_),numberRows_));
    8989#endif
    9090  link_ = ClpCopyOfArray(rhs.link_,numberRows_);
     
    171171#else
    172172    // actually long double
    173     workDouble_ = (double *) ClpCopyOfArray((longWork *) rhs.workDouble_,numberRows_);
     173    workDouble_ = reinterpret_cast<double *> (ClpCopyOfArray(reinterpret_cast<CoinWorkDouble *> (rhs.workDouble_),numberRows_));
    174174#endif
    175175    link_ = ClpCopyOfArray(rhs.link_,numberRows_);
     
    195195   region1 is rows+columns, region2 is rows */
    196196void
    197 ClpCholeskyBase::solveKKT (double * region1, double * region2, const double * diagonal,
    198                            double diagonalScaleFactor)
     197ClpCholeskyBase::solveKKT (CoinWorkDouble * region1, CoinWorkDouble * region2, const CoinWorkDouble * diagonal,
     198                           CoinWorkDouble diagonalScaleFactor)
    199199{
    200200  if (!doKKT_) {
     
    202202    int numberColumns = model_->numberColumns();
    203203    int numberTotal = numberRows_+numberColumns;
    204     double * region1Save = new double[numberTotal];
     204    CoinWorkDouble * region1Save = new CoinWorkDouble[numberTotal];
    205205    for (iColumn=0;iColumn<numberTotal;iColumn++) {
    206206      region1[iColumn] *= diagonal[iColumn];
     
    209209    multiplyAdd(region1+numberColumns,numberRows_,-1.0,region2,1.0);
    210210    model_->clpMatrix()->times(1.0,region1,region2);
    211     double maximumRHS = maximumAbsElement(region2,numberRows_);
    212     double scale=1.0;
    213     double unscale=1.0;
     211    CoinWorkDouble maximumRHS = maximumAbsElement(region2,numberRows_);
     212    CoinWorkDouble scale=1.0;
     213    CoinWorkDouble unscale=1.0;
    214214    if (maximumRHS>1.0e-30) {
    215215      if (maximumRHS<=0.5) {
    216         double factor=2.0;
     216        CoinWorkDouble factor=2.0;
    217217        while (maximumRHS<=0.5) {
    218218          maximumRHS*=factor;
     
    220220        } /* endwhile */
    221221      } else if (maximumRHS>=2.0&&maximumRHS<=COIN_DBL_MAX) {
    222         double factor=0.5;
     222        CoinWorkDouble factor=0.5;
    223223        while (maximumRHS>=2.0) {
    224224          maximumRHS*=factor;
     
    246246    int numberColumns = model_->numberColumns();
    247247    int numberTotal = numberColumns + numberRowsModel;
    248     double * array = new double [numberRows_];
     248    CoinWorkDouble * array = new CoinWorkDouble [numberRows_];
    249249    CoinMemcpyN(region1,numberTotal,array);
    250250    CoinMemcpyN(region2,numberRowsModel,array+numberTotal);
     
    253253    int iRow;
    254254    for (iRow=0;iRow<numberTotal;iRow++) {
    255       if (rowsDropped_[iRow]&&fabs(array[iRow])>1.0e-8) {
     255      if (rowsDropped_[iRow]&&CoinAbs(array[iRow])>1.0e-8) {
    256256        printf("row region1 %d dropped %g\n",iRow,array[iRow]);
    257257      }
    258258    }
    259259    for (;iRow<numberRows_;iRow++) {
    260       if (rowsDropped_[iRow]&&fabs(array[iRow])>1.0e-8) {
     260      if (rowsDropped_[iRow]&&CoinAbs(array[iRow])>1.0e-8) {
    261261        printf("row region2 %d dropped %g\n",iRow,array[iRow]);
    262262      }
     
    12371237#else
    12381238    // actually long double
    1239     workDouble_ = (double *) (new longWork[numberRows_]);
     1239    workDouble_ = reinterpret_cast<double *> (new CoinWorkDouble[numberRows_]);
    12401240#endif
    12411241    diagonal_ = new longDouble[numberRows_];
     
    15211521/* Factorize - filling in rowsDropped and returning number dropped */
    15221522int
    1523 ClpCholeskyBase::factorize(const double * diagonal, int * rowsDropped)
     1523ClpCholeskyBase::factorize(const CoinWorkDouble * diagonal, int * rowsDropped)
    15241524{
    15251525  const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
     
    15331533  int numberColumns=model_->clpMatrix()->getNumCols();
    15341534  //perturbation
    1535   longDouble perturbation=model_->diagonalPerturbation()*model_->diagonalNorm();
     1535  CoinWorkDouble perturbation=model_->diagonalPerturbation()*model_->diagonalNorm();
    15361536  //perturbation=perturbation*perturbation*100000000.0;
    15371537  if (perturbation>1.0) {
     
    15401540      std::cout <<"large perturbation "<<perturbation<<std::endl;
    15411541#endif
    1542     perturbation=sqrt(perturbation);;
     1542    perturbation=CoinSqrt(perturbation);
    15431543    perturbation=1.0;
    15441544  }
     
    15481548  CoinZeroN(work,numberRows_);
    15491549  int newDropped=0;
    1550   double largest=1.0;
    1551   double smallest=COIN_DBL_MAX;
     1550  CoinWorkDouble largest=1.0;
     1551  CoinWorkDouble smallest=COIN_DBL_MAX;
    15521552  int numberDense=0;
    15531553  if (!doKKT_) {
    1554     const double * diagonalSlack = diagonal + numberColumns;
     1554    const CoinWorkDouble * diagonalSlack = diagonal + numberColumns;
    15551555    if (dense_)
    15561556      numberDense=dense_->numberRows();
     
    15771577      }
    15781578    }
    1579     longDouble delta2 = model_->delta(); // add delta*delta to diagonal
     1579    CoinWorkDouble delta2 = model_->delta(); // add delta*delta to diagonal
    15801580    delta2 *= delta2;
    15811581    // largest in initial matrix
    1582     double largest2=1.0e-20;
     1582    CoinWorkDouble largest2=1.0e-20;
    15831583    for (iRow=0;iRow<numberRows_;iRow++) {
    15841584      longDouble * put = sparseFactor_+choleskyStart_[iRow];
     
    15971597            CoinBigIndex start=columnStart[iColumn];
    15981598            CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
    1599             longDouble multiplier = diagonal[iColumn]*elementByRow[k];
     1599            CoinWorkDouble multiplier = diagonal[iColumn]*elementByRow[k];
    16001600            for (CoinBigIndex j=start;j<end;j++) {
    16011601              int jRow=row[j];
    16021602              int jNewRow = permuteInverse_[jRow];
    16031603              if (jNewRow>=iRow&&!rowsDropped_[jRow]) {
    1604                 longDouble value=element[j]*multiplier;
     1604                CoinWorkDouble value=element[j]*multiplier;
    16051605                work[jNewRow] += value;
    16061606              }
     
    16091609        }
    16101610        diagonal_[iRow]=work[iRow];
    1611         largest2 = CoinMax(largest2,fabs(work[iRow]));
     1611        largest2 = CoinMax(largest2,CoinAbs(work[iRow]));
    16121612        work[iRow]=0.0;
    16131613        int j;
     
    16151615          int jRow = which[j];
    16161616          put[j]=work[jRow];
    1617           largest2 = CoinMax(largest2,fabs(work[jRow]));
     1617          largest2 = CoinMax(largest2,CoinAbs(work[jRow]));
    16181618          work[jRow]=0.0;
    16191619        }
     
    16291629    //check sizes
    16301630    largest2*=1.0e-20;
    1631     largest = CoinMin(largest2,1.0e-11);
     1631    largest = CoinMin(largest2,CHOL_SMALL_VALUE);
    16321632    int numberDroppedBefore=0;
    16331633    for (iRow=0;iRow<numberRows_;iRow++) {
     
    16361636      rowsDropped[iRow]=dropped;
    16371637      if (!dropped) {
    1638         longDouble diagonal = diagonal_[iRow];
     1638        CoinWorkDouble diagonal = diagonal_[iRow];
    16391639        if (diagonal>largest2) {
    16401640          diagonal_[iRow]=diagonal+perturbation;
     
    16631663      for ( i=0;i<numberRows_;i++) {
    16641664        assert (diagonal_[i]>=0.0);
    1665         diagonal_[i]=sqrt(diagonal_[i]);
     1665        diagonal_[i]=CoinSqrt(diagonal_[i]);
    16661666      }
    16671667      // Update dense columns (just L)
     
    16731673            a[j]=0.0;
    16741674        }
    1675         solveLong(a,1);
     1675        for (i=0;i<numberRows_;i++) {
     1676          int iRow = permute_[i];
     1677          workDouble_[i] = a[iRow];
     1678        }
     1679        for (i=0;i<numberRows_;i++) {
     1680          CoinWorkDouble value=workDouble_[i];
     1681          CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
     1682          CoinBigIndex j;
     1683          for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
     1684            int iRow = choleskyRow_[j+offset];
     1685            workDouble_[iRow] -= sparseFactor_[j]*value;
     1686          }
     1687        }
     1688        for (i=0;i<numberRows_;i++) {
     1689          int iRow = permute_[i];
     1690          a[iRow]=workDouble_[i]*diagonal_[i];
     1691        }
    16761692      }
    16771693      dense_->resetRowsDropped();
     
    16821698        const longDouble * a = denseColumn_+i*numberRows_;
    16831699        // do diagonal
    1684         longDouble value = denseDiagonal[i];
     1700        CoinWorkDouble value = denseDiagonal[i];
    16851701        const longDouble * b = denseColumn_+i*numberRows_;
    16861702        for (int k=0;k<numberRows_;k++)
     
    16881704        denseDiagonal[i]=value;
    16891705        for (int j=i+1;j<numberDense;j++) {
    1690           longDouble value = 0.0;
     1706          CoinWorkDouble value = 0.0;
    16911707          const longDouble * b = denseColumn_+j*numberRows_;
    16921708          for (int k=0;k<numberRows_;k++)
     
    17781794    }
    17791795    if (permuted) {
    1780       longDouble delta2 = model_->delta(); // add delta*delta to bottom
     1796      CoinWorkDouble delta2 = model_->delta(); // add delta*delta to bottom
    17811797      delta2 *= delta2;
    17821798      // Need to permute - ugly
     
    17881804          if (iOriginalRow<numberColumns) {
    17891805            iColumn=iOriginalRow;
    1790             longDouble value = diagonal[iColumn];
    1791             if (fabs(value)>1.0e-100) {
     1806            CoinWorkDouble value = diagonal[iColumn];
     1807            if (CoinAbs(value)>1.0e-100) {
    17921808              value = 1.0/value;
    1793               largest = CoinMax(largest,fabs(value));
     1809              largest = CoinMax(largest,CoinAbs(value));
    17941810              diagonal_[iRow] = -value;
    17951811              CoinBigIndex start=columnStart[iColumn];
     
    18001816                if (kRow>iRow) {
    18011817                  work[kRow]=element[j];
    1802                   largest = CoinMax(largest,fabs(element[j]));
     1818                  largest = CoinMax(largest,CoinAbs(element[j]));
    18031819                }
    18041820              }
     
    18071823            }
    18081824          } else if (iOriginalRow<numberTotal) {
    1809             longDouble value = diagonal[iOriginalRow];
    1810             if (fabs(value)>1.0e-100) {
     1825            CoinWorkDouble value = diagonal[iOriginalRow];
     1826            if (CoinAbs(value)>1.0e-100) {
    18111827              value = 1.0/value;
    1812               largest = CoinMax(largest,fabs(value));
     1828              largest = CoinMax(largest,CoinAbs(value));
    18131829            } else {
    18141830              value = 1.0e100;
     
    18281844              if (jNewRow>iRow) {
    18291845                work[jNewRow]=elementByRow[j];
    1830                 largest = CoinMax(largest,fabs(elementByRow[j]));
     1846                largest = CoinMax(largest,CoinAbs(elementByRow[j]));
    18311847              }
    18321848            }
     
    18571873            CoinBigIndex j;
    18581874            iColumn=iOriginalRow;
    1859             longDouble value = diagonal[iColumn];
    1860             if (fabs(value)>1.0e-100) {
     1875            CoinWorkDouble value = diagonal[iColumn];
     1876            if (CoinAbs(value)>1.0e-100) {
    18611877              value = 1.0/value;
    18621878              for (j=columnQuadraticStart[iColumn];
     
    18701886                }
    18711887              }
    1872               largest = CoinMax(largest,fabs(value));
     1888              largest = CoinMax(largest,CoinAbs(value));
    18731889              diagonal_[iRow] = -value;
    18741890              CoinBigIndex start=columnStart[iColumn];
     
    18791895                if (kRow>iRow) {
    18801896                  work[kRow]=element[j];
    1881                   largest = CoinMax(largest,fabs(element[j]));
     1897                  largest = CoinMax(largest,CoinAbs(element[j]));
    18821898                }
    18831899              }
     
    18861902            }
    18871903          } else if (iOriginalRow<numberTotal) {
    1888             longDouble value = diagonal[iOriginalRow];
    1889             if (fabs(value)>1.0e-100) {
     1904            CoinWorkDouble value = diagonal[iOriginalRow];
     1905            if (CoinAbs(value)>1.0e-100) {
    18901906              value = 1.0/value;
    1891               largest = CoinMax(largest,fabs(value));
     1907              largest = CoinMax(largest,CoinAbs(value));
    18921908            } else {
    18931909              value = 1.0e100;
     
    19071923              if (jNewRow>iRow) {
    19081924                work[jNewRow]=elementByRow[j];
    1909                 largest = CoinMax(largest,fabs(elementByRow[j]));
     1925                largest = CoinMax(largest,CoinAbs(elementByRow[j]));
    19101926              }
    19111927            }
     
    19311947          longDouble * put = sparseFactor_+choleskyStart_[iColumn];
    19321948          int * which = choleskyRow_+indexStart_[iColumn];
    1933           longDouble value = diagonal[iColumn];
    1934           if (fabs(value)>1.0e-100) {
     1949          CoinWorkDouble value = diagonal[iColumn];
     1950          if (CoinAbs(value)>1.0e-100) {
    19351951            value = 1.0/value;
    1936             largest = CoinMax(largest,fabs(value));
     1952            largest = CoinMax(largest,CoinAbs(value));
    19371953            diagonal_[iColumn] = -value;
    19381954            CoinBigIndex start=columnStart[iColumn];
     
    19421958              //sparseFactor_[numberElements++]=element[j];
    19431959              work[row[j]+numberTotal]=element[j];
    1944               largest = CoinMax(largest,fabs(element[j]));
     1960              largest = CoinMax(largest,CoinAbs(element[j]));
    19451961            }
    19461962          } else {
     
    19651981          int * which = choleskyRow_+indexStart_[iColumn];
    19661982          int number = choleskyStart_[iColumn+1]-choleskyStart_[iColumn];
    1967           longDouble value = diagonal[iColumn];
     1983          CoinWorkDouble value = diagonal[iColumn];
    19681984          CoinBigIndex j;
    1969           if (fabs(value)>1.0e-100) {
     1985          if (CoinAbs(value)>1.0e-100) {
    19701986            value = 1.0/value;
    19711987            for (j=columnQuadraticStart[iColumn];
     
    19781994              }
    19791995            }
    1980             largest = CoinMax(largest,fabs(value));
     1996            largest = CoinMax(largest,CoinAbs(value));
    19811997            diagonal_[iColumn] = -value;
    19821998            CoinBigIndex start=columnStart[iColumn];
     
    19842000            for (j=start;j<end;j++) {
    19852001              work[row[j]+numberTotal]=element[j];
    1986               largest = CoinMax(largest,fabs(element[j]));
     2002              largest = CoinMax(largest,CoinAbs(element[j]));
    19872003            }
    19882004            for (j=0;j<number;j++) {
     
    20052021        longDouble * put = sparseFactor_+choleskyStart_[iColumn];
    20062022        int * which = choleskyRow_+indexStart_[iColumn];
    2007         longDouble value = diagonal[iColumn];
    2008         if (fabs(value)>1.0e-100) {
     2023        CoinWorkDouble value = diagonal[iColumn];
     2024        if (CoinAbs(value)>1.0e-100) {
    20092025          value = 1.0/value;
    2010           largest = CoinMax(largest,fabs(value));
     2026          largest = CoinMax(largest,CoinAbs(value));
    20112027        } else {
    20122028          value = 1.0e100;
     
    20232039      }
    20242040      // Finish diagonal
    2025       longDouble delta2 = model_->delta(); // add delta*delta to bottom
     2041      CoinWorkDouble delta2 = model_->delta(); // add delta*delta to bottom
    20262042      delta2 *= delta2;
    20272043      for (iRow=0;iRow<numberRowsModel;iRow++) {
     
    20372053    //check sizes
    20382054    largest*=1.0e-20;
    2039     largest = CoinMin(largest,1.0e-11);
     2055    largest = CoinMin(largest,CHOL_SMALL_VALUE);
    20402056    doubleParameters_[10]=CoinMax(1.0e-20,largest);
    20412057    integerParameters_[20]=0;
     
    20562072    // Should save adjustments in ..R_
    20572073    int n1=0,n2=0;
    2058     double * primalR = model_->primalR();
    2059     double * dualR = model_->dualR();
     2074    CoinWorkDouble * primalR = model_->primalR();
     2075    CoinWorkDouble * dualR = model_->dualR();
    20602076    for (iRow=0;iRow<numberTotal;iRow++) {
    20612077      if (rowsDropped2[iRow]) {
     
    20932109ClpCholeskyBase::factorizePart2(int * rowsDropped)
    20942110{
    2095   longDouble largest=doubleParameters_[3];
    2096   longDouble smallest=doubleParameters_[4];
     2111  CoinWorkDouble largest=doubleParameters_[3];
     2112  CoinWorkDouble smallest=doubleParameters_[4];
    20972113  int numberDropped=integerParameters_[20];
    20982114  // probably done before
     
    21592175      for ( jRow=lastRow;jRow<iRow;jRow++) {
    21602176        int jCount = jRow-lastRow;
    2161         longDouble diagonalValue = diagonal_[jRow];
     2177        CoinWorkDouble diagonalValue = diagonal_[jRow];
    21622178        CoinBigIndex start=choleskyStart_[jRow];
    21632179        CoinBigIndex end=choleskyStart_[jRow+1];
     
    21652181          jCount--;
    21662182          CoinBigIndex get = choleskyStart_[kRow]+jCount;
    2167           longDouble a_jk = sparseFactor_[get];
    2168           longDouble value1 = d[kRow]*a_jk;
     2183          CoinWorkDouble a_jk = sparseFactor_[get];
     2184          CoinWorkDouble value1 = d[kRow]*a_jk;
    21692185          diagonalValue -= a_jk*value1;
    21702186          for (CoinBigIndex j=start;j<end;j++)
     
    22222238    }
    22232239    // for each column L[*,kRow] that affects L[*,iRow]
    2224     longDouble diagonalValue=diagonal_[iRow];
     2240    CoinWorkDouble diagonalValue=diagonal_[iRow];
    22252241    int nextRow = link_[iRow];
    22262242    int kRow=0;
     
    22342250      CoinBigIndex end=choleskyStart_[kRow+1];
    22352251      assert(k<end);
    2236       longDouble a_ik=sparseFactor_[k++];
    2237       longDouble value1=d[kRow]*a_ik;
     2252      CoinWorkDouble a_ik=sparseFactor_[k++];
     2253      CoinWorkDouble value1=d[kRow]*a_ik;
    22382254      // update first
    22392255      first[kRow]=k;
     
    22592275            CoinBigIndex j=first[kkRow];
    22602276            //int iiRow = choleskyRow_[j+indexStart_[kkRow]-choleskyStart_[kkRow]];
    2261             longDouble a = sparseFactor_[j];
    2262             longDouble dValue = d[kkRow]*a;
     2277            CoinWorkDouble a = sparseFactor_[j];
     2278            CoinWorkDouble dValue = d[kkRow]*a;
    22632279            diagonalValue -= a*dValue;
    22642280            work[kkRow]=dValue;
     
    22712287          for (int i=0;i<length;i++) {
    22722288            int lRow = choleskyRow_[currentIndex++];
    2273             longDouble t0 = work[lRow];
     2289            CoinWorkDouble t0 = work[lRow];
    22742290            for (int kkRow=kRow;kkRow<last;kkRow++) {
    22752291              CoinBigIndex j = first[kkRow]+i;
     
    23402356        for (CoinBigIndex j=start;j<end;j++) {
    23412357          int jRow = choleskyRow_[j+offset];
    2342           longDouble value = sparseFactor_[j]-work[jRow];
     2358          CoinWorkDouble value = sparseFactor_[j]-work[jRow];
    23432359          work[jRow]=0.0;
    23442360          sparseFactor_[j]= diagonalValue*value;
     
    23952411      CoinBigIndex offset = indexStart_[iRow]-choleskyStart_[iRow];
    23962412      if (clique_[iRow]<2) {
    2397         longDouble dValue=d[iRow];
     2413        CoinWorkDouble dValue=d[iRow];
    23982414        for (CoinBigIndex k=start;k<end;k++) {
    23992415          int kRow = choleskyRow_[k+offset];
    24002416          assert(kRow>=firstDense_);
    2401           longDouble a_ik=sparseFactor_[k];
    2402           longDouble value1=dValue*a_ik;
     2417          CoinWorkDouble a_ik=sparseFactor_[k];
     2418          CoinWorkDouble value1=dValue*a_ik;
    24032419          diagonal_[kRow] -= value1*a_ik;
    24042420          CoinBigIndex base = choleskyStart_[kRow]-kRow-1;
    24052421          for (CoinBigIndex j=k+1;j<end;j++) {
    24062422            int jRow=choleskyRow_[j+offset];
    2407             longDouble a_jk = sparseFactor_[j];
     2423            CoinWorkDouble a_jk = sparseFactor_[j];
    24082424            sparseFactor_[base+jRow] -= a_jk*value1;
    24092425          }
     
    24112427      } else if (clique_[iRow]<3) {
    24122428        // do as pair
    2413         longDouble dValue0=d[iRow];
    2414         longDouble dValue1=d[iRow+1];
     2429        CoinWorkDouble dValue0=d[iRow];
     2430        CoinWorkDouble dValue1=d[iRow+1];
    24152431        int offset1 = first[iRow+1]-start;
    24162432        // skip row
     
    24192435          int kRow = choleskyRow_[k+offset];
    24202436          assert(kRow>=firstDense_);
    2421           longDouble a_ik0=sparseFactor_[k];
    2422           longDouble value0=dValue0*a_ik0;
    2423           longDouble a_ik1=sparseFactor_[k+offset1];
    2424           longDouble value1=dValue1*a_ik1;
     2437          CoinWorkDouble a_ik0=sparseFactor_[k];
     2438          CoinWorkDouble value0=dValue0*a_ik0;
     2439          CoinWorkDouble a_ik1=sparseFactor_[k+offset1];
     2440          CoinWorkDouble value1=dValue1*a_ik1;
    24252441          diagonal_[kRow] -= value0*a_ik0 + value1*a_ik1;
    24262442          CoinBigIndex base = choleskyStart_[kRow]-kRow-1;
    24272443          for (CoinBigIndex j=k+1;j<end;j++) {
    24282444            int jRow=choleskyRow_[j+offset];
    2429             longDouble a_jk0 = sparseFactor_[j];
    2430             longDouble a_jk1 = sparseFactor_[j+offset1];
     2445            CoinWorkDouble a_jk0 = sparseFactor_[j];
     2446            CoinWorkDouble a_jk1 = sparseFactor_[j+offset1];
    24312447            sparseFactor_[base+jRow] -= a_jk0*value0+a_jk1*value1;
    24322448          }
     
    24412457        // maybe later get fancy on big cliques and do transpose copy
    24422458        // seems only worth going to 3 on Intel
    2443         longDouble dValue0=d[iRow];
    2444         longDouble dValue1=d[iRow+1];
    2445         longDouble dValue2=d[iRow+2];
     2459        CoinWorkDouble dValue0=d[iRow];
     2460        CoinWorkDouble dValue1=d[iRow+1];
     2461        CoinWorkDouble dValue2=d[iRow+2];
    24462462        // get offsets and skip rows
    24472463        int offset1 = first[++iRow]-start;
     
    24502466          int kRow = choleskyRow_[k+offset];
    24512467          assert(kRow>=firstDense_);
    2452           double diagonalValue=diagonal_[kRow];
    2453           longDouble a_ik0=sparseFactor_[k];
    2454           longDouble value0=dValue0*a_ik0;
    2455           longDouble a_ik1=sparseFactor_[k+offset1];
    2456           longDouble value1=dValue1*a_ik1;
    2457           longDouble a_ik2=sparseFactor_[k+offset2];
    2458           longDouble value2=dValue2*a_ik2;
     2468          CoinWorkDouble diagonalValue=diagonal_[kRow];
     2469          CoinWorkDouble a_ik0=sparseFactor_[k];
     2470          CoinWorkDouble value0=dValue0*a_ik0;
     2471          CoinWorkDouble a_ik1=sparseFactor_[k+offset1];
     2472          CoinWorkDouble value1=dValue1*a_ik1;
     2473          CoinWorkDouble a_ik2=sparseFactor_[k+offset2];
     2474          CoinWorkDouble value2=dValue2*a_ik2;
    24592475          CoinBigIndex base = choleskyStart_[kRow]-kRow-1;
    24602476          diagonal_[kRow] = diagonalValue - value0*a_ik0 - value1*a_ik1 - value2*a_ik2;
    24612477          for (CoinBigIndex j=k+1;j<end;j++) {
    24622478            int jRow=choleskyRow_[j+offset];
    2463             longDouble a_jk0 = sparseFactor_[j];
    2464             longDouble a_jk1 = sparseFactor_[j+offset1];
    2465             longDouble a_jk2 = sparseFactor_[j+offset2];
     2479            CoinWorkDouble a_jk0 = sparseFactor_[j];
     2480            CoinWorkDouble a_jk1 = sparseFactor_[j+offset1];
     2481            CoinWorkDouble a_jk2 = sparseFactor_[j+offset2];
    24662482            sparseFactor_[base+jRow] -= a_jk0*value0+a_jk1*value1+a_jk2*value2;
    24672483          }
     
    24722488        // maybe later get fancy on big cliques and do transpose copy
    24732489        // maybe only worth going to 3 on Intel (but may have hidden registers)
    2474         longDouble dValue0=d[iRow];
    2475         longDouble dValue1=d[iRow+1];
    2476         longDouble dValue2=d[iRow+2];
    2477         longDouble dValue3=d[iRow+3];
     2490        CoinWorkDouble dValue0=d[iRow];
     2491        CoinWorkDouble dValue1=d[iRow+1];
     2492        CoinWorkDouble dValue2=d[iRow+2];
     2493        CoinWorkDouble dValue3=d[iRow+3];
    24782494        // get offsets and skip rows
    24792495        int offset1 = first[++iRow]-start;
     
    24832499          int kRow = choleskyRow_[k+offset];
    24842500          assert(kRow>=firstDense_);
    2485           double diagonalValue=diagonal_[kRow];
    2486           longDouble a_ik0=sparseFactor_[k];
    2487           longDouble value0=dValue0*a_ik0;
    2488           longDouble a_ik1=sparseFactor_[k+offset1];
    2489           longDouble value1=dValue1*a_ik1;
    2490           longDouble a_ik2=sparseFactor_[k+offset2];
    2491           longDouble value2=dValue2*a_ik2;
    2492           longDouble a_ik3=sparseFactor_[k+offset3];
    2493           longDouble value3=dValue3*a_ik3;
     2501          CoinWorkDouble diagonalValue=diagonal_[kRow];
     2502          CoinWorkDouble a_ik0=sparseFactor_[k];
     2503          CoinWorkDouble value0=dValue0*a_ik0;
     2504          CoinWorkDouble a_ik1=sparseFactor_[k+offset1];
     2505          CoinWorkDouble value1=dValue1*a_ik1;
     2506          CoinWorkDouble a_ik2=sparseFactor_[k+offset2];
     2507          CoinWorkDouble value2=dValue2*a_ik2;
     2508          CoinWorkDouble a_ik3=sparseFactor_[k+offset3];
     2509          CoinWorkDouble value3=dValue3*a_ik3;
    24942510          CoinBigIndex base = choleskyStart_[kRow]-kRow-1;
    24952511          diagonal_[kRow] = diagonalValue - (value0*a_ik0 + value1*a_ik1 + value2*a_ik2+value3*a_ik3);
    24962512          for (CoinBigIndex j=k+1;j<end;j++) {
    24972513            int jRow=choleskyRow_[j+offset];
    2498             longDouble a_jk0 = sparseFactor_[j];
    2499             longDouble a_jk1 = sparseFactor_[j+offset1];
    2500             longDouble a_jk2 = sparseFactor_[j+offset2];
    2501             longDouble a_jk3 = sparseFactor_[j+offset3];
     2514            CoinWorkDouble a_jk0 = sparseFactor_[j];
     2515            CoinWorkDouble a_jk1 = sparseFactor_[j+offset1];
     2516            CoinWorkDouble a_jk2 = sparseFactor_[j+offset2];
     2517            CoinWorkDouble a_jk3 = sparseFactor_[j+offset3];
    25022518            sparseFactor_[base+jRow] -= a_jk0*value0+a_jk1*value1+a_jk2*value2+a_jk3*value3;
    25032519          }
     
    25202536    // do change;
    25212537    int numberDense = dense_->numberRows();
    2522     double * change = new double[numberDense];
     2538    CoinWorkDouble * change = new CoinWorkDouble[numberDense];
    25232539    for (i=0;i<numberDense;i++) {
    25242540      const longDouble * a = denseColumn_+i*numberRows_;
    2525       longDouble value =0.0;
     2541      CoinWorkDouble value =0.0;
    25262542      for (int iRow=0;iRow<numberRows_;iRow++)
    25272543        value += a[iRow]*region[iRow];
     
    25292545    }
    25302546    // solve
     2547#if CLP_LONG_CHOLESKY>0
     2548    dense_->solveLong(change);
     2549#else
    25312550    dense_->solve(change);
     2551#endif
    25322552    for (i=0;i<numberDense;i++) {
    25332553      const longDouble * a = denseColumn_+i*numberRows_;
    2534       longDouble value = change[i];
     2554      CoinWorkDouble value = change[i];
    25352555      for (int iRow=0;iRow<numberRows_;iRow++)
    25362556        region[iRow] -= value*a[iRow];
     
    25492569#ifdef CLP_DEBUG
    25502570  double * regionX=NULL;
    2551   if (sizeof(longWork)!=sizeof(double)&&type==3) {
     2571  if (sizeof(CoinWorkDouble)!=sizeof(double)&&type==3) {
    25522572    regionX=ClpCopyOfArray(region,numberRows_);
    25532573  }
    25542574#endif
    2555   longWork * work = reinterpret_cast<longWork *> (workDouble_);
     2575  CoinWorkDouble * work = reinterpret_cast<CoinWorkDouble *> (workDouble_);
    25562576  int i;
    25572577  CoinBigIndex j;
     
    25632583  case 1:
    25642584    for (i=0;i<numberRows_;i++) {
    2565       longDouble value=work[i];
     2585      CoinWorkDouble value=work[i];
    25662586      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
    25672587      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
     
    25782598    for (i=numberRows_-1;i>=0;i--) {
    25792599      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
    2580       longDouble value=work[i]*diagonal_[i];
     2600      CoinWorkDouble value=work[i]*diagonal_[i];
    25812601      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
    25822602        int iRow = choleskyRow_[j+offset];
     
    25912611    for (i=0;i<firstDense_;i++) {
    25922612      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
    2593       longDouble value=work[i];
     2613      CoinWorkDouble value=work[i];
    25942614      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
    25952615        int iRow = choleskyRow_[j+offset];
     
    26092629#endif
    26102630      for (i=numberRows_-1;i>=firstDense_;i--) {
    2611         longDouble value=work[i];
     2631        CoinWorkDouble value=work[i];
    26122632        int iRow = permute_[i];
    26132633        region[iRow]=value;
     
    26162636    for (i=firstDense_-1;i>=0;i--) {
    26172637      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
    2618       longDouble value=work[i]*diagonal_[i];
     2638      CoinWorkDouble value=work[i]*diagonal_[i];
    26192639      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
    26202640        int iRow = choleskyRow_[j+offset];
     
    26342654    double largestO=0.0;
    26352655    for (i=0;i<numberRows_;i++) {
    2636       largestO = CoinMax(largestO,fabs(regionX[i]));
     2656      largestO = CoinMax(largestO,CoinAbs(regionX[i]));
    26372657    }
    26382658    for (i=0;i<numberRows_;i++) {
     
    26422662    for (i=0;i<firstDense_;i++) {
    26432663      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
    2644       longDouble value=work[i];
     2664      CoinWorkDouble value=work[i];
    26452665      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
    26462666        int iRow = choleskyRow_[j+offset];
     
    26562676      dense.solveLong(work+firstDense_);
    26572677      for (i=numberRows_-1;i>=firstDense_;i--) {
    2658         longDouble value=work[i];
     2678        CoinWorkDouble value=work[i];
    26592679        int iRow = permute_[i];
    26602680        regionX[iRow]=value;
     
    26632683    for (i=firstDense_-1;i>=0;i--) {
    26642684      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
    2665       longDouble value=work[i]*diagonal_[i];
     2685      CoinWorkDouble value=work[i]*diagonal_[i];
    26662686      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
    26672687        int iRow = choleskyRow_[j+offset];
     
    26752695    double largestV=0.0;
    26762696    for (i=0;i<numberRows_;i++) {
    2677       largest = CoinMax(largest,fabs(region[i]-regionX[i]));
    2678       largestV = CoinMax(largestV,fabs(region[i]));
     2697      largest = CoinMax(largest,CoinAbs(region[i]-regionX[i]));
     2698      largestV = CoinMax(largestV,CoinAbs(region[i]));
    26792699    }
    26802700    printf("largest difference %g, largest %g, largest original %g\n",
     
    26842704#endif
    26852705}
    2686 /* solve - 1 just first half, 2 just second half - 3 both.
    2687    If 1 and 2 then diagonal has sqrt of inverse otherwise inverse
    2688 */
     2706#if CLP_LONG_CHOLESKY
     2707/* Uses factorization to solve. */
    26892708void
    2690 ClpCholeskyBase::solveLong(longDouble * region, int type)
     2709ClpCholeskyBase::solve (CoinWorkDouble * region)
    26912710{
     2711  assert (!whichDense_) ;
     2712  CoinWorkDouble * work = reinterpret_cast<CoinWorkDouble *> (workDouble_);
    26922713  int i;
    26932714  CoinBigIndex j;
    26942715  for (i=0;i<numberRows_;i++) {
    26952716    int iRow = permute_[i];
    2696     workDouble_[i] = region[iRow];
    2697   }
    2698   switch (type) {
    2699   case 1:
    2700     for (i=0;i<numberRows_;i++) {
    2701       longDouble value=workDouble_[i];
    2702       CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
    2703       for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
    2704         int iRow = choleskyRow_[j+offset];
    2705         workDouble_[iRow] -= sparseFactor_[j]*value;
    2706       }
    2707     }
    2708     for (i=0;i<numberRows_;i++) {
    2709       int iRow = permute_[i];
    2710       region[iRow]=workDouble_[i]*diagonal_[i];
    2711     }
    2712     break;
    2713   case 2:
    2714     for (i=numberRows_-1;i>=0;i--) {
    2715       CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
    2716       longDouble value=workDouble_[i]*diagonal_[i];
    2717       for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
    2718         int iRow = choleskyRow_[j+offset];
    2719         value -= sparseFactor_[j]*workDouble_[iRow];
    2720       }
    2721       workDouble_[i]=value;
     2717    work[i] = region[iRow];
     2718  }
     2719  for (i=0;i<firstDense_;i++) {
     2720    CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
     2721    CoinWorkDouble value=work[i];
     2722    for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
     2723      int iRow = choleskyRow_[j+offset];
     2724      work[iRow] -= sparseFactor_[j]*value;
     2725    }
     2726  }
     2727  if (firstDense_<numberRows_) {
     2728    // do dense
     2729    ClpCholeskyDense dense;
     2730    // just borrow space
     2731    int nDense = numberRows_-firstDense_;
     2732    dense.reserveSpace(this,nDense);
     2733#if CLP_LONG_CHOLESKY!=1
     2734    dense.solveLong(work+firstDense_);
     2735#else
     2736    dense.solveLongWork(work+firstDense_);
     2737#endif
     2738    for (i=numberRows_-1;i>=firstDense_;i--) {
     2739      CoinWorkDouble value=work[i];
    27222740      int iRow = permute_[i];
    27232741      region[iRow]=value;
    27242742    }
    2725     break;
    2726   case 3:
    2727     for (i=0;i<firstDense_;i++) {
    2728       CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
    2729       longDouble value=workDouble_[i];
    2730       for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
    2731         int iRow = choleskyRow_[j+offset];
    2732         workDouble_[iRow] -= sparseFactor_[j]*value;
    2733       }
    2734     }
    2735     if (firstDense_<numberRows_) {
    2736       // do dense
    2737       ClpCholeskyDense dense;
    2738       // just borrow space
    2739       int nDense = numberRows_-firstDense_;
    2740       dense.reserveSpace(this,nDense);
    2741       dense.solveLong(workDouble_+firstDense_);
    2742       for (i=numberRows_-1;i>=firstDense_;i--) {
    2743         longDouble value=workDouble_[i];
    2744         int iRow = permute_[i];
    2745         region[iRow]=value;
    2746       }
    2747     }
    2748     for (i=firstDense_-1;i>=0;i--) {
    2749       CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
    2750       longDouble value=workDouble_[i]*diagonal_[i];
    2751       for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
    2752         int iRow = choleskyRow_[j+offset];
    2753         value -= sparseFactor_[j]*workDouble_[iRow];
    2754       }
    2755       workDouble_[i]=value;
    2756       int iRow = permute_[i];
    2757       region[iRow]=value;
    2758     }
    2759     break;
     2743  }
     2744  for (i=firstDense_-1;i>=0;i--) {
     2745    CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
     2746    CoinWorkDouble value=work[i]*diagonal_[i];
     2747    for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
     2748      int iRow = choleskyRow_[j+offset];
     2749      value -= sparseFactor_[j]*work[iRow];
     2750    }
     2751    work[i]=value;
     2752    int iRow = permute_[i];
     2753    region[iRow]=value;
    27602754  }
    27612755}
    2762 
     2756#endif
  • trunk/Clp/src/ClpCholeskyBase.hpp

    r1055 r1367  
    66#include "CoinPragma.hpp"
    77#include "CoinFinite.hpp"
     8//#define CLP_LONG_CHOLESKY 0
     9#ifndef CLP_LONG_CHOLESKY
    810#define CLP_LONG_CHOLESKY 0
     11#endif
     12/* valid combinations are
     13   CLP_LONG_CHOLESKY 0 and COIN_LONG_WORK 0
     14   CLP_LONG_CHOLESKY 1 and COIN_LONG_WORK 1
     15   CLP_LONG_CHOLESKY 2 and COIN_LONG_WORK 1
     16*/
     17#if COIN_LONG_WORK==0
     18#if CLP_LONG_CHOLESKY>0
     19#define CHOLESKY_BAD_COMBINATION
     20#endif
     21#else
     22#if CLP_LONG_CHOLESKY==0
     23#define CHOLESKY_BAD_COMBINATION
     24#endif
     25#endif
     26#ifdef CHOLESKY_BAD_COMBINATION
     27#  warning("Bad combination of CLP_LONG_CHOLESKY and COIN_BIG_DOUBLE/COIN_LONG_WORK");
     28"Bad combination of CLP_LONG_CHOLESKY and COIN_LONG_WORK"
     29#endif
    930#if CLP_LONG_CHOLESKY>1
    1031typedef long double longDouble;
    11 typedef long double longWork;
     32#define CHOL_SMALL_VALUE 1.0e-15
    1233#elif CLP_LONG_CHOLESKY==1
    1334typedef double longDouble;
    14 typedef long double longWork;
     35#define CHOL_SMALL_VALUE 1.0e-11
    1536#else
    1637typedef double longDouble;
    17 typedef double longWork;
     38#define CHOL_SMALL_VALUE 1.0e-11
    1839#endif
    1940class ClpInterior;
     
    4667  /** Factorize - filling in rowsDropped and returning number dropped.
    4768      If return code negative then out of memory */
    48   virtual int factorize(const double * diagonal, int * rowsDropped) ;
     69  virtual int factorize(const CoinWorkDouble * diagonal, int * rowsDropped) ;
    4970  /** Uses factorization to solve. */
    5071  virtual void solve (double * region) ;
    5172  /** Uses factorization to solve. - given as if KKT.
    5273   region1 is rows+columns, region2 is rows */
    53   virtual void solveKKT (double * region1, double * region2, const double * diagonal,
    54                          double diagonalScaleFactor);
     74  virtual void solveKKT (CoinWorkDouble * region1, CoinWorkDouble * region2, const CoinWorkDouble * diagonal,
     75                         CoinWorkDouble diagonalScaleFactor);
     76#if CLP_LONG_CHOLESKY
     77  /** Uses factorization to solve. */
     78  virtual void solve (CoinWorkDouble * region) ;
     79#endif
    5580  //@}
    5681
     
    163188  */
    164189  void solve(double * region, int type);
    165   void solveLong(longDouble * region, int type);
    166190  /// Forms ADAT - returns nonzero if not enough memory
    167191  int preOrder(bool lowerTriangular, bool includeDiagonal, bool doKKT);
  • trunk/Clp/src/ClpCholeskyDense.cpp

    r1321 r1367  
    66#include "CoinPragma.hpp"
    77#include "CoinHelperFunctions.hpp"
     8#include "ClpHelperFunctions.hpp"
    89
    910#include "ClpInterior.hpp"
     
    165166  CoinZeroN(sparseFactor_,sizeFactor_);
    166167  //perturbation
    167   longDouble perturbation=model_->diagonalPerturbation()*model_->diagonalNorm();
     168  CoinWorkDouble perturbation=model_->diagonalPerturbation()*model_->diagonalNorm();
    168169  perturbation=perturbation*perturbation;
    169170  if (perturbation>1.0) {
     
    172173      std::cout <<"large perturbation "<<perturbation<<std::endl;
    173174#endif
    174     perturbation=sqrt(perturbation);;
     175    perturbation=CoinSqrt(perturbation);;
    175176    perturbation=1.0;
    176177  }
    177178  int iRow;
    178179  int newDropped=0;
    179   double largest=1.0;
    180   double smallest=COIN_DBL_MAX;
    181   longDouble delta2 = model_->delta(); // add delta*delta to diagonal
     180  CoinWorkDouble largest=1.0;
     181  CoinWorkDouble smallest=COIN_DBL_MAX;
     182  CoinWorkDouble delta2 = model_->delta(); // add delta*delta to diagonal
    182183  delta2 *= delta2;
    183184  if (!doKKT_) {
     
    187188    const double * diagonalSlack = diagonal + numberColumns;
    188189    // largest in initial matrix
    189     double largest2=1.0e-20;
     190    CoinWorkDouble largest2=1.0e-20;
    190191    for (iRow=0;iRow<numberRows_;iRow++) {
    191192      if (!rowsDropped_[iRow]) {
    192193        CoinBigIndex startRow=rowStart[iRow];
    193194        CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
    194         longDouble diagonalValue = diagonalSlack[iRow]+delta2;
     195        CoinWorkDouble diagonalValue = diagonalSlack[iRow]+delta2;
    195196        for (CoinBigIndex k=startRow;k<endRow;k++) {
    196197          int iColumn=column[k];
    197198          CoinBigIndex start=columnStart[iColumn];
    198199          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
    199           longDouble multiplier = diagonal[iColumn]*elementByRow[k];
     200          CoinWorkDouble multiplier = diagonal[iColumn]*elementByRow[k];
    200201          for (CoinBigIndex j=start;j<end;j++) {
    201202            int jRow=row[j];
    202203            if (!rowsDropped_[jRow]) {
    203204              if (jRow>iRow) {
    204                 longDouble value=element[j]*multiplier;
     205                CoinWorkDouble value=element[j]*multiplier;
    205206                work[jRow] += value;
    206207              } else if (jRow==iRow) {
    207                 longDouble value=element[j]*multiplier;
     208                CoinWorkDouble value=element[j]*multiplier;
    208209                diagonalValue += value;
    209210              }
     
    212213        }
    213214        for (int j=iRow+1;j<numberRows_;j++)
    214           largest2 = CoinMax(largest2,fabs(work[j]));
     215          largest2 = CoinMax(largest2,CoinAbs(work[j]));
    215216        diagonal_[iRow]=diagonalValue;
    216         largest2 = CoinMax(largest2,fabs(diagonalValue));
     217        largest2 = CoinMax(largest2,CoinAbs(diagonalValue));
    217218      } else {
    218219        // dropped
     
    224225    //check sizes
    225226    largest2*=1.0e-20;
    226     largest = CoinMin(largest2,1.0e-11);
     227    largest = CoinMin(largest2,CHOL_SMALL_VALUE);
    227228    int numberDroppedBefore=0;
    228229    for (iRow=0;iRow<numberRows_;iRow++) {
     
    231232      rowsDropped[iRow]=dropped;
    232233      if (!dropped) {
    233         longDouble diagonal = diagonal_[iRow];
     234        CoinWorkDouble diagonal = diagonal_[iRow];
    234235        if (diagonal>largest2) {
    235236          diagonal_[iRow]=diagonal+perturbation;
     
    289290    if (!quadratic) {
    290291      for (iColumn=0;iColumn<numberColumns;iColumn++) {
    291         longDouble value = diagonal[iColumn];
    292         if (fabs(value)>1.0e-100) {
     292        CoinWorkDouble value = diagonal[iColumn];
     293        if (CoinAbs(value)>1.0e-100) {
    293294          value = 1.0/value;
    294           largest = CoinMax(largest,fabs(value));
     295          largest = CoinMax(largest,CoinAbs(value));
    295296          diagonal_[iColumn] = -value;
    296297          CoinBigIndex start=columnStart[iColumn];
     
    300301            //sparseFactor_[numberElements++]=element[j];
    301302            work[row[j]+numberTotal]=element[j];
    302             largest = CoinMax(largest,fabs(element[j]));
     303            largest = CoinMax(largest,CoinAbs(element[j]));
    303304          }
    304305        } else {
     
    315316      const double * quadraticElement = quadratic->getElements();
    316317      for (iColumn=0;iColumn<numberColumns;iColumn++) {
    317         longDouble value = diagonal[iColumn];
     318        CoinWorkDouble value = diagonal[iColumn];
    318319        CoinBigIndex j;
    319         if (fabs(value)>1.0e-100) {
     320        if (CoinAbs(value)>1.0e-100) {
    320321          value = 1.0/value;
    321322          for (j=columnQuadraticStart[iColumn];
     
    328329            }
    329330          }
    330           largest = CoinMax(largest,fabs(value));
     331          largest = CoinMax(largest,CoinAbs(value));
    331332          diagonal_[iColumn] = -value;
    332333          CoinBigIndex start=columnStart[iColumn];
     
    334335          for (j=start;j<end;j++) {
    335336            work[row[j]+numberTotal]=element[j];
    336             largest = CoinMax(largest,fabs(element[j]));
     337            largest = CoinMax(largest,CoinAbs(element[j]));
    337338          }
    338339        } else {
     
    346347    // slacks
    347348    for (iColumn=numberColumns;iColumn<numberTotal;iColumn++) {
    348       longDouble value = diagonal[iColumn];
    349       if (fabs(value)>1.0e-100) {
     349      CoinWorkDouble value = diagonal[iColumn];
     350      if (CoinAbs(value)>1.0e-100) {
    350351        value = 1.0/value;
    351         largest = CoinMax(largest,fabs(value));
     352        largest = CoinMax(largest,CoinAbs(value));
    352353      } else {
    353354        value = 1.0e100;
     
    364365    //check sizes
    365366    largest*=1.0e-20;
    366     largest = CoinMin(largest,1.0e-11);
     367    largest = CoinMin(largest,CHOL_SMALL_VALUE);
    367368    doubleParameters_[10]=CoinMax(1.0e-20,largest);
    368369    integerParameters_[20]=0;
     
    387388    // Should save adjustments in ..R_
    388389    int n1=0,n2=0;
    389     double * primalR = model_->primalR();
    390     double * dualR = model_->dualR();
     390    CoinWorkDouble * primalR = model_->primalR();
     391    CoinWorkDouble * dualR = model_->dualR();
    391392    for (iRow=0;iRow<numberTotal;iRow++) {
    392393      if (rowsDropped2[iRow]) {
     
    429430  CoinMemcpyN(yy,numberRows_,diagonal_);
    430431  int numberDropped=0;
    431   double largest=0.0;
    432   double smallest=COIN_DBL_MAX;
     432  CoinWorkDouble largest=0.0;
     433  CoinWorkDouble smallest=COIN_DBL_MAX;
    433434  double dropValue = doubleParameters_[10];
    434435  int firstPositive=integerParameters_[34];
     
    441442    int addOffsetNow = numberRows_-1;;
    442443    longDouble * workNow=sparseFactor_-1+iColumn;
    443     double diagonalValue = diagonal_[iColumn];
     444    CoinWorkDouble diagonalValue = diagonal_[iColumn];
    444445    for (iRow=0;iRow<iColumn;iRow++) {
    445446      double aj = *workNow;
    446447      addOffsetNow--;
    447448      workNow += addOffsetNow;
    448       double multiplier = workDouble_[iRow];
    449       diagonalValue -= aj*aj*multiplier;
     449      diagonalValue -= aj*aj*workDouble_[iRow];
    450450    }
    451451    bool dropColumn=false;
     
    776776                longDouble * diagonal, longDouble * work, int * rowsDropped)
    777777{
    778   longDouble largest=doubleParameters_[3];
    779   longDouble smallest=doubleParameters_[4];
     778  CoinWorkDouble largest=doubleParameters_[3];
     779  CoinWorkDouble smallest=doubleParameters_[4];
    780780  double dropValue = doubleParameters_[10];
    781781  int firstPositive=integerParameters_[34];
     
    783783  int numberDropped=0;
    784784  int i, j, k;
    785   longDouble t00,temp1;
     785  CoinWorkDouble t00,temp1;
    786786  longDouble * aa;
    787787  aa = a-BLOCK;
     
    790790    t00 = aa[j];
    791791    for (k = 0; k < j; ++k) {
    792       longDouble multiplier = work[k];
     792      CoinWorkDouble multiplier = work[k];
    793793      t00 -= a[j + k * BLOCK] * a[j + k * BLOCK] * multiplier;
    794794    }
    795795    bool dropColumn=false;
    796     longDouble useT00=t00;
     796    CoinWorkDouble useT00=t00;
    797797    if (j+rowOffset<firstPositive) {
    798798      // must be negative
     
    831831        t00=aa[i];
    832832        for (k = 0; k < j; ++k) {
    833           longDouble multiplier = work[k];
     833          CoinWorkDouble multiplier = work[k];
    834834          t00 -= a[i + k * BLOCK] * a[j + k * BLOCK] * multiplier;
    835835        }
     
    874874    for (j = 0; j < BLOCK; j +=2) {
    875875      aa+=2*BLOCK;
    876       longDouble temp0 = diagonal[j];
    877       longDouble temp1 = diagonal[j+1];
     876      CoinWorkDouble temp0 = diagonal[j];
     877      CoinWorkDouble temp1 = diagonal[j+1];
    878878      for (int i=0;i<BLOCK;i+=2) {
    879         longDouble t00=aUnder[i+j*BLOCK];
    880         longDouble t10=aUnder[i+ BLOCK + j*BLOCK];
    881         longDouble t01=aUnder[i+1+j*BLOCK];
    882         longDouble t11=aUnder[i+1+ BLOCK + j*BLOCK];
     879        CoinWorkDouble t00=aUnder[i+j*BLOCK];
     880        CoinWorkDouble t10=aUnder[i+ BLOCK + j*BLOCK];
     881        CoinWorkDouble t01=aUnder[i+1+j*BLOCK];
     882        CoinWorkDouble t11=aUnder[i+1+ BLOCK + j*BLOCK];
    883883        for (int k = 0; k < j; ++k) {
    884           longDouble multiplier=work[k];
    885           longDouble au0 = aUnder[i + k * BLOCK]*multiplier;
    886           longDouble au1 = aUnder[i + 1 + k * BLOCK]*multiplier;
    887           longDouble at0 = aTri[j + k * BLOCK];
    888           longDouble at1 = aTri[j + 1 + k * BLOCK];
     884          CoinWorkDouble multiplier=work[k];
     885          CoinWorkDouble au0 = aUnder[i + k * BLOCK]*multiplier;
     886          CoinWorkDouble au1 = aUnder[i + 1 + k * BLOCK]*multiplier;
     887          CoinWorkDouble at0 = aTri[j + k * BLOCK];
     888          CoinWorkDouble at1 = aTri[j + 1 + k * BLOCK];
    889889          t00 -= au0 * at0;
    890890          t10 -= au0 * at1;
     
    893893        }
    894894        t00 *= temp0;
    895         longDouble at1 = aTri[j + 1 + j*BLOCK]*work[j];
     895        CoinWorkDouble at1 = aTri[j + 1 + j*BLOCK]*work[j];
    896896        t10 -= t00 * at1;
    897897        t01 *= temp0;
     
    908908    for (j = 0; j < BLOCK; j ++) {
    909909      aa+=BLOCK;
    910       longDouble temp1 = diagonal[j];
     910      CoinWorkDouble temp1 = diagonal[j];
    911911      for (int i=0;i<nUnder;i++) {
    912         longDouble t00=aUnder[i+j*BLOCK];
     912        CoinWorkDouble t00=aUnder[i+j*BLOCK];
    913913        for (int k = 0; k < j; ++k) {
    914           longDouble multiplier=work[k];
     914          CoinWorkDouble multiplier=work[k];
    915915          t00 -= aUnder[i + k * BLOCK] * aTri[j + k * BLOCK] * multiplier;
    916916        }
     
    937937#endif
    938938  int i, j, k;
    939   longDouble t00;
     939  CoinWorkDouble t00;
    940940  longDouble * aa;
    941941#ifdef BLOCKUNROLL
     
    944944    aa = aTri-2*BLOCK;
    945945    for (j = 0; j < BLOCK; j +=2) {
    946       longDouble t00,t01,t10,t11;
     946      CoinWorkDouble t00,t01,t10,t11;
    947947      aa+=2*BLOCK;
    948948      aUnder2+=2;
     
    951951      t10=aa[j+1+BLOCK];
    952952      for (k = 0; k < BLOCK; ++k) {
    953         longDouble multiplier = work[k];
    954         longDouble a0 = aUnder2[k * BLOCK];
    955         longDouble a1 = aUnder2[1 + k * BLOCK];
    956         longDouble x0 = a0 * multiplier;
    957         longDouble x1 = a1 * multiplier;
     953        CoinWorkDouble multiplier = work[k];
     954        CoinWorkDouble a0 = aUnder2[k * BLOCK];
     955        CoinWorkDouble a1 = aUnder2[1 + k * BLOCK];
     956        CoinWorkDouble x0 = a0 * multiplier;
     957        CoinWorkDouble x1 = a1 * multiplier;
    958958        t00 -= a0 * x0;
    959959        t01 -= a1 * x0;
     
    969969        t11=aa[i+1+BLOCK];
    970970        for (k = 0; k < BLOCK; ++k) {
    971           longDouble multiplier = work[k];
    972           longDouble a0 = aUnder2[k * BLOCK]*multiplier;
    973           longDouble a1 = aUnder2[1 + k * BLOCK]*multiplier;
     971          CoinWorkDouble multiplier = work[k];
     972          CoinWorkDouble a0 = aUnder2[k * BLOCK]*multiplier;
     973          CoinWorkDouble a1 = aUnder2[1 + k * BLOCK]*multiplier;
    974974          t00 -= aUnder[i + k * BLOCK] * a0;
    975975          t01 -= aUnder[i + k * BLOCK] * a1;
     
    991991        t00=aa[i];
    992992        for (k = 0; k < BLOCK; ++k) {
    993           longDouble multiplier = work[k];
     993          CoinWorkDouble multiplier = work[k];
    994994          t00 -= aUnder[i + k * BLOCK] * aUnder[j + k * BLOCK]*multiplier;
    995995        }
     
    10341034      aa+=2*BLOCK;
    10351035      for (i=0;i< BLOCK;i+=2) {
    1036         longDouble t00=aa[i+0*BLOCK];
    1037         longDouble t10=aa[i+1*BLOCK];
    1038         longDouble t01=aa[i+1+0*BLOCK];
    1039         longDouble t11=aa[i+1+1*BLOCK];
     1036        CoinWorkDouble t00=aa[i+0*BLOCK];
     1037        CoinWorkDouble t10=aa[i+1*BLOCK];
     1038        CoinWorkDouble t01=aa[i+1+0*BLOCK];
     1039        CoinWorkDouble t11=aa[i+1+1*BLOCK];
    10401040        for (k=0;k<BLOCK;k++) {
    1041           longDouble multiplier = work[k];
    1042           longDouble a00=aUnder[i+k*BLOCK]*multiplier;
    1043           longDouble a01=aUnder[i+1+k*BLOCK]*multiplier;
     1041          CoinWorkDouble multiplier = work[k];
     1042          CoinWorkDouble a00=aUnder[i+k*BLOCK]*multiplier;
     1043          CoinWorkDouble a01=aUnder[i+1+k*BLOCK]*multiplier;
    10441044          t00 -= a00 * above[j + 0 + k * BLOCK];
    10451045          t10 -= a00 * above[j + 1 + k * BLOCK];
     
    10571057      aa+=4*BLOCK;
    10581058      for (i=0;i< BLOCK;i+=4) {
    1059         longDouble t00=aa[i+0+0*BLOCK];
    1060         longDouble t10=aa[i+0+1*BLOCK];
    1061         longDouble t20=aa[i+0+2*BLOCK];
    1062         longDouble t30=aa[i+0+3*BLOCK];
    1063         longDouble t01=aa[i+1+0*BLOCK];
    1064         longDouble t11=aa[i+1+1*BLOCK];
    1065         longDouble t21=aa[i+1+2*BLOCK];
    1066         longDouble t31=aa[i+1+3*BLOCK];
    1067         longDouble t02=aa[i+2+0*BLOCK];
    1068         longDouble t12=aa[i+2+1*BLOCK];
    1069         longDouble t22=aa[i+2+2*BLOCK];
    1070         longDouble t32=aa[i+2+3*BLOCK];
    1071         longDouble t03=aa[i+3+0*BLOCK];
    1072         longDouble t13=aa[i+3+1*BLOCK];
    1073         longDouble t23=aa[i+3+2*BLOCK];
    1074         longDouble t33=aa[i+3+3*BLOCK];
     1059        CoinWorkDouble t00=aa[i+0+0*BLOCK];
     1060        CoinWorkDouble t10=aa[i+0+1*BLOCK];
     1061        CoinWorkDouble t20=aa[i+0+2*BLOCK];
     1062        CoinWorkDouble t30=aa[i+0+3*BLOCK];
     1063        CoinWorkDouble t01=aa[i+1+0*BLOCK];
     1064        CoinWorkDouble t11=aa[i+1+1*BLOCK];
     1065        CoinWorkDouble t21=aa[i+1+2*BLOCK];
     1066        CoinWorkDouble t31=aa[i+1+3*BLOCK];
     1067        CoinWorkDouble t02=aa[i+2+0*BLOCK];
     1068        CoinWorkDouble t12=aa[i+2+1*BLOCK];
     1069        CoinWorkDouble t22=aa[i+2+2*BLOCK];
     1070        CoinWorkDouble t32=aa[i+2+3*BLOCK];
     1071        CoinWorkDouble t03=aa[i+3+0*BLOCK];
     1072        CoinWorkDouble t13=aa[i+3+1*BLOCK];
     1073        CoinWorkDouble t23=aa[i+3+2*BLOCK];
     1074        CoinWorkDouble t33=aa[i+3+3*BLOCK];
    10751075        const longDouble * COIN_RESTRICT aUnderNow = aUnder+i;
    10761076        const longDouble * COIN_RESTRICT aboveNow = above+j;
    10771077        for (k=0;k<BLOCK;k++) {
    1078           longDouble multiplier = work[k];
    1079           longDouble a00=aUnderNow[0]*multiplier;
    1080           longDouble a01=aUnderNow[1]*multiplier;
    1081           longDouble a02=aUnderNow[2]*multiplier;
    1082           longDouble a03=aUnderNow[3]*multiplier;
     1078          CoinWorkDouble multiplier = work[k];
     1079          CoinWorkDouble a00=aUnderNow[0]*multiplier;
     1080          CoinWorkDouble a01=aUnderNow[1]*multiplier;
     1081          CoinWorkDouble a02=aUnderNow[2]*multiplier;
     1082          CoinWorkDouble a03=aUnderNow[3]*multiplier;
    10831083          t00 -= a00 * aboveNow[0];
    10841084          t10 -= a00 * aboveNow[1];
     
    11251125      aa+=4*BLOCK;
    11261126      for (i=0;i< n;i+=2) {
    1127         longDouble t00=aa[i+0*BLOCK];
    1128         longDouble t10=aa[i+1*BLOCK];
    1129         longDouble t20=aa[i+2*BLOCK];
    1130         longDouble t30=aa[i+3*BLOCK];
    1131         longDouble t01=aa[i+1+0*BLOCK];
    1132         longDouble t11=aa[i+1+1*BLOCK];
    1133         longDouble t21=aa[i+1+2*BLOCK];
    1134         longDouble t31=aa[i+1+3*BLOCK];
     1127        CoinWorkDouble t00=aa[i+0*BLOCK];
     1128        CoinWorkDouble t10=aa[i+1*BLOCK];
     1129        CoinWorkDouble t20=aa[i+2*BLOCK];
     1130        CoinWorkDouble t30=aa[i+3*BLOCK];
     1131        CoinWorkDouble t01=aa[i+1+0*BLOCK];
     1132        CoinWorkDouble t11=aa[i+1+1*BLOCK];
     1133        CoinWorkDouble t21=aa[i+1+2*BLOCK];
     1134        CoinWorkDouble t31=aa[i+1+3*BLOCK];
    11351135        const longDouble * COIN_RESTRICT aUnderNow = aUnder+i;
    11361136        const longDouble * COIN_RESTRICT aboveNow = above+j;
    11371137        for (k=0;k<BLOCK;k++) {
    1138           longDouble multiplier = work[k];
    1139           longDouble a00=aUnderNow[0]*multiplier;
    1140           longDouble a01=aUnderNow[1]*multiplier;
     1138          CoinWorkDouble multiplier = work[k];
     1139          CoinWorkDouble a00=aUnderNow[0]*multiplier;
     1140          CoinWorkDouble a01=aUnderNow[1]*multiplier;
    11411141          t00 -= a00 * aboveNow[0];
    11421142          t10 -= a00 * aboveNow[1];
     
    11601160      }
    11611161      if (odd) {
    1162         longDouble t0=aa[n+0*BLOCK];
    1163         longDouble t1=aa[n+1*BLOCK];
    1164         longDouble t2=aa[n+2*BLOCK];
    1165         longDouble t3=aa[n+3*BLOCK];
    1166         longDouble a0;
     1162        CoinWorkDouble t0=aa[n+0*BLOCK];
     1163        CoinWorkDouble t1=aa[n+1*BLOCK];
     1164        CoinWorkDouble t2=aa[n+2*BLOCK];
     1165        CoinWorkDouble t3=aa[n+3*BLOCK];
     1166        CoinWorkDouble a0;
    11671167        for (k=0;k<BLOCK;k++) {
    11681168          a0=aUnder[n+k*BLOCK]*work[k];
     
    11841184    aa+=BLOCK;
    11851185    for (i=0;i< nUnder;i++) {
    1186       longDouble t00=aa[i+0*BLOCK];
     1186      CoinWorkDouble t00=aa[i+0*BLOCK];
    11871187      for (k=0;k<BLOCK;k++) {
    1188         longDouble a00=aUnder[i+k*BLOCK]*work[k];
     1188        CoinWorkDouble a00=aUnder[i+k*BLOCK]*work[k];
    11891189        t00 -= a00 * above[j + k * BLOCK];
    11901190      }
     
    12931293  if (numberRows_<200) {
    12941294    for (int i=0;i<numberRows_;i++) {
    1295       assert(fabs(region[i]-region2[i])<1.0e-3);
     1295      assert(CoinAbs(region[i]-region2[i])<1.0e-3);
    12961296    }
    12971297    delete [] region2;
     
    13041304{
    13051305  int j, k;
    1306   longDouble t00;
     1306  CoinWorkDouble t00;
    13071307  for (j = 0; j < n; j ++) {
    13081308    t00 = region[j];
     
    13221322  if (n==BLOCK) {
    13231323    for (k = 0; k < BLOCK; k+=4) {
    1324       longDouble t0 = region2[0];
    1325       longDouble t1 = region2[1];
    1326       longDouble t2 = region2[2];
    1327       longDouble t3 = region2[3];
     1324      CoinWorkDouble t0 = region2[0];
     1325      CoinWorkDouble t1 = region2[1];
     1326      CoinWorkDouble t2 = region2[2];
     1327      CoinWorkDouble t3 = region2[3];
    13281328      t0 -= region[0] * a[0 + 0 * BLOCK];
    13291329      t1 -= region[0] * a[1 + 0 * BLOCK];
     
    14161416#endif
    14171417    for (k = 0; k < n; ++k) {
    1418       longDouble t00 = region2[k];
     1418      CoinWorkDouble t00 = region2[k];
    14191419      for (j = 0; j < BLOCK; j ++) {
    14201420        t00 -= region[j] * a[k + j * BLOCK];
     
    14311431{
    14321432  int j, k;
    1433   longDouble t00;
     1433  CoinWorkDouble t00;
    14341434  for (j = n-1; j >=0; j --) {
    14351435    t00 = region[j];
     
    14491449  if (n==BLOCK) {
    14501450    for (j = 0; j < BLOCK; j +=4) {
    1451       longDouble t0 = region[0];
    1452       longDouble t1 = region[1];
    1453       longDouble t2 = region[2];
    1454       longDouble t3 = region[3];
     1451      CoinWorkDouble t0 = region[0];
     1452      CoinWorkDouble t1 = region[1];
     1453      CoinWorkDouble t2 = region[2];
     1454      CoinWorkDouble t3 = region[3];
    14551455      t0 -= region2[0] * a[0 + 0*BLOCK];
    14561456      t1 -= region2[0] * a[0 + 1*BLOCK];
     
    15441544#endif
    15451545    for (j = 0; j < BLOCK; j ++) {
    1546       longDouble t00 = region[j];
     1546      CoinWorkDouble t00 = region[j];
    15471547      for (k = 0; k < n; ++k) {
    15481548        t00 -= region2[k] * a[k + j * BLOCK];
     
    15561556/* Uses factorization to solve. */
    15571557void
    1558 ClpCholeskyDense::solveLong (longDouble * region)
     1558ClpCholeskyDense::solveLong (CoinWorkDouble * region)
    15591559{
    15601560  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
     
    16191619// Forward part of solve 1
    16201620void
    1621 ClpCholeskyDense::solveF1Long(longDouble * a,int n,longDouble * region)
     1621ClpCholeskyDense::solveF1Long(longDouble * a,int n,CoinWorkDouble * region)
    16221622{
    16231623  int j, k;
    1624   longDouble t00;
     1624  CoinWorkDouble t00;
    16251625  for (j = 0; j < n; j ++) {
    16261626    t00 = region[j];
     
    16341634// Forward part of solve 2
    16351635void
    1636 ClpCholeskyDense::solveF2Long(longDouble * a,int n,longDouble * region, longDouble * region2)
     1636ClpCholeskyDense::solveF2Long(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2)
    16371637{
    16381638  int j, k;
     
    16401640  if (n==BLOCK) {
    16411641    for (k = 0; k < BLOCK; k+=4) {
    1642       longDouble t0 = region2[0];
    1643       longDouble t1 = region2[1];
    1644       longDouble t2 = region2[2];
    1645       longDouble t3 = region2[3];
     1642      CoinWorkDouble t0 = region2[0];
     1643      CoinWorkDouble t1 = region2[1];
     1644      CoinWorkDouble t2 = region2[2];
     1645      CoinWorkDouble t3 = region2[3];
    16461646      t0 -= region[0] * a[0 + 0 * BLOCK];
    16471647      t1 -= region[0] * a[1 + 0 * BLOCK];
     
    17341734#endif
    17351735    for (k = 0; k < n; ++k) {
    1736       longDouble t00 = region2[k];
     1736      CoinWorkDouble t00 = region2[k];
    17371737      for (j = 0; j < BLOCK; j ++) {
    17381738        t00 -= region[j] * a[k + j * BLOCK];
     
    17461746// Backward part of solve 1
    17471747void
    1748 ClpCholeskyDense::solveB1Long(longDouble * a,int n,longDouble * region)
     1748ClpCholeskyDense::solveB1Long(longDouble * a,int n,CoinWorkDouble * region)
    17491749{
    17501750  int j, k;
    1751   longDouble t00;
     1751  CoinWorkDouble t00;
    17521752  for (j = n-1; j >=0; j --) {
    17531753    t00 = region[j];
     
    17611761// Backward part of solve 2
    17621762void
    1763 ClpCholeskyDense::solveB2Long(longDouble * a,int n,longDouble * region, longDouble * region2)
     1763ClpCholeskyDense::solveB2Long(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2)
    17641764{
    17651765  int j, k;
     
    17671767  if (n==BLOCK) {
    17681768    for (j = 0; j < BLOCK; j +=4) {
    1769       longDouble t0 = region[0];
    1770       longDouble t1 = region[1];
    1771       longDouble t2 = region[2];
    1772       longDouble t3 = region[3];
     1769      CoinWorkDouble t0 = region[0];
     1770      CoinWorkDouble t1 = region[1];
     1771      CoinWorkDouble t2 = region[2];
     1772      CoinWorkDouble t3 = region[3];
    17731773      t0 -= region2[0] * a[0 + 0*BLOCK];
    17741774      t1 -= region2[0] * a[0 + 1*BLOCK];
     
    18621862#endif
    18631863    for (j = 0; j < BLOCK; j ++) {
    1864       longDouble t00 = region[j];
     1864      CoinWorkDouble t00 = region[j];
    18651865      for (k = 0; k < n; ++k) {
    18661866        t00 -= region2[k] * a[k + j * BLOCK];
     
    18741874/* Uses factorization to solve. */
    18751875void
    1876 ClpCholeskyDense::solveLongWork (longWork * region)
     1876ClpCholeskyDense::solveLongWork (CoinWorkDouble * region)
    18771877{
    18781878  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
     
    19371937// Forward part of solve 1
    19381938void
    1939 ClpCholeskyDense::solveF1LongWork(longDouble * a,int n,longWork * region)
     1939ClpCholeskyDense::solveF1LongWork(longDouble * a,int n,CoinWorkDouble * region)
    19401940{
    19411941  int j, k;
    1942   longWork t00;
     1942  CoinWorkDouble t00;
    19431943  for (j = 0; j < n; j ++) {
    19441944    t00 = region[j];
     
    19521952// Forward part of solve 2
    19531953void
    1954 ClpCholeskyDense::solveF2LongWork(longDouble * a,int n,longWork * region, longWork * region2)
     1954ClpCholeskyDense::solveF2LongWork(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2)
    19551955{
    19561956  int j, k;
     
    19581958  if (n==BLOCK) {
    19591959    for (k = 0; k < BLOCK; k+=4) {
    1960       longWork t0 = region2[0];
    1961       longWork t1 = region2[1];
    1962       longWork t2 = region2[2];
    1963       longWork t3 = region2[3];
     1960      CoinWorkDouble t0 = region2[0];
     1961      CoinWorkDouble t1 = region2[1];
     1962      CoinWorkDouble t2 = region2[2];
     1963      CoinWorkDouble t3 = region2[3];
    19641964      t0 -= region[0] * a[0 + 0 * BLOCK];
    19651965      t1 -= region[0] * a[1 + 0 * BLOCK];
     
    20522052#endif
    20532053    for (k = 0; k < n; ++k) {
    2054       longWork t00 = region2[k];
     2054      CoinWorkDouble t00 = region2[k];
    20552055      for (j = 0; j < BLOCK; j ++) {
    20562056        t00 -= region[j] * a[k + j * BLOCK];
     
    20642064// Backward part of solve 1
    20652065void
    2066 ClpCholeskyDense::solveB1LongWork(longDouble * a,int n,longWork * region)
     2066ClpCholeskyDense::solveB1LongWork(longDouble * a,int n,CoinWorkDouble * region)
    20672067{
    20682068  int j, k;
    2069   longWork t00;
     2069  CoinWorkDouble t00;
    20702070  for (j = n-1; j >=0; j --) {
    20712071    t00 = region[j];
     
    20792079// Backward part of solve 2
    20802080void
    2081 ClpCholeskyDense::solveB2LongWork(longDouble * a,int n,longWork * region, longWork * region2)
     2081ClpCholeskyDense::solveB2LongWork(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2)
    20822082{
    20832083  int j, k;
     
    20852085  if (n==BLOCK) {
    20862086    for (j = 0; j < BLOCK; j +=4) {
    2087       longWork t0 = region[0];
    2088       longWork t1 = region[1];
    2089       longWork t2 = region[2];
    2090       longWork t3 = region[3];
     2087      CoinWorkDouble t0 = region[0];
     2088      CoinWorkDouble t1 = region[1];
     2089      CoinWorkDouble t2 = region[2];
     2090      CoinWorkDouble t3 = region[3];
    20912091      t0 -= region2[0] * a[0 + 0*BLOCK];
    20922092      t1 -= region2[0] * a[0 + 1*BLOCK];
     
    21802180#endif
    21812181    for (j = 0; j < BLOCK; j ++) {
    2182       longWork t00 = region[j];
     2182      CoinWorkDouble t00 = region[j];
    21832183      for (k = 0; k < n; ++k) {
    21842184        t00 -= region2[k] * a[k + j * BLOCK];
  • trunk/Clp/src/ClpCholeskyDense.hpp

    r1321 r1367  
    8989  void solveB2(longDouble * a,int n,double * region,double * region2);
    9090  /** Uses factorization to solve. */
    91   void solveLong (longDouble * region) ;
     91  void solveLong (CoinWorkDouble * region) ;
    9292  /// Forward part of solve
    93   void solveF1Long(longDouble * a,int n,longDouble * region);
    94   void solveF2Long(longDouble * a,int n,longDouble * region,longDouble * region2);
     93  void solveF1Long(longDouble * a,int n,CoinWorkDouble * region);
     94  void solveF2Long(longDouble * a,int n,CoinWorkDouble * region,CoinWorkDouble * region2);
    9595  /// Backward part of solve
    96   void solveB1Long(longDouble * a,int n,longDouble * region);
    97   void solveB2Long(longDouble * a,int n,longDouble * region,longDouble * region2);
     96  void solveB1Long(longDouble * a,int n,CoinWorkDouble * region);
     97  void solveB2Long(longDouble * a,int n,CoinWorkDouble * region,CoinWorkDouble * region2);
    9898  /** Uses factorization to solve. */
    99   void solveLongWork (longWork * region) ;
     99  void solveLongWork (CoinWorkDouble * region) ;
    100100  /// Forward part of solve
    101   void solveF1LongWork(longDouble * a,int n,longWork * region);
    102   void solveF2LongWork(longDouble * a,int n,longWork * region,longWork * region2);
     101  void solveF1LongWork(longDouble * a,int n,CoinWorkDouble * region);
     102  void solveF2LongWork(longDouble * a,int n,CoinWorkDouble * region,CoinWorkDouble * region2);
    103103  /// Backward part of solve
    104   void solveB1LongWork(longDouble * a,int n,longWork * region);
    105   void solveB2LongWork(longDouble * a,int n,longWork * region,longWork * region2);
     104  void solveB1LongWork(longDouble * a,int n,CoinWorkDouble * region);
     105  void solveB2LongWork(longDouble * a,int n,CoinWorkDouble * region,CoinWorkDouble * region2);
    106106  int bNumber(const longDouble * array,int &, int&);
    107107  /// A
  • trunk/Clp/src/ClpHelperFunctions.cpp

    r1185 r1367  
    112112  }
    113113}
     114#if COIN_LONG_WORK
     115// For long double versions
     116CoinWorkDouble
     117maximumAbsElement(const CoinWorkDouble * region, int size)
     118{
     119  int i;
     120  CoinWorkDouble maxValue=0.0;
     121  for (i=0;i<size;i++)
     122    maxValue = CoinMax(maxValue,CoinAbs(region[i]));
     123  return maxValue;
     124}
     125void
     126setElements(CoinWorkDouble * region, int size, CoinWorkDouble value)
     127{
     128  int i;
     129  for (i=0;i<size;i++)
     130    region[i]=value;
     131}
     132void
     133multiplyAdd(const CoinWorkDouble * region1, int size, CoinWorkDouble multiplier1,
     134                 CoinWorkDouble * region2, CoinWorkDouble multiplier2)
     135{
     136  int i;
     137  if (multiplier1==1.0) {
     138    if (multiplier2==1.0) {
     139      for (i=0;i<size;i++)
     140        region2[i] = region1[i] + region2[i];
     141    } else if (multiplier2==-1.0) {
     142      for (i=0;i<size;i++)
     143        region2[i] = region1[i] - region2[i];
     144    } else if (multiplier2==0.0) {
     145      for (i=0;i<size;i++)
     146        region2[i] = region1[i] ;
     147    } else {
     148      for (i=0;i<size;i++)
     149        region2[i] = region1[i] + multiplier2*region2[i];
     150    }
     151  } else if (multiplier1==-1.0) {
     152    if (multiplier2==1.0) {
     153      for (i=0;i<size;i++)
     154        region2[i] = -region1[i] + region2[i];
     155    } else if (multiplier2==-1.0) {
     156      for (i=0;i<size;i++)
     157        region2[i] = -region1[i] - region2[i];
     158    } else if (multiplier2==0.0) {
     159      for (i=0;i<size;i++)
     160        region2[i] = -region1[i] ;
     161    } else {
     162      for (i=0;i<size;i++)
     163        region2[i] = -region1[i] + multiplier2*region2[i];
     164    }
     165  } else if (multiplier1==0.0) {
     166    if (multiplier2==1.0) {
     167      // nothing to do
     168    } else if (multiplier2==-1.0) {
     169      for (i=0;i<size;i++)
     170        region2[i] =  -region2[i];
     171    } else if (multiplier2==0.0) {
     172      for (i=0;i<size;i++)
     173        region2[i] =  0.0;
     174    } else {
     175      for (i=0;i<size;i++)
     176        region2[i] =  multiplier2*region2[i];
     177    }
     178  } else {
     179    if (multiplier2==1.0) {
     180      for (i=0;i<size;i++)
     181        region2[i] = multiplier1*region1[i] + region2[i];
     182    } else if (multiplier2==-1.0) {
     183      for (i=0;i<size;i++)
     184        region2[i] = multiplier1*region1[i] - region2[i];
     185    } else if (multiplier2==0.0) {
     186      for (i=0;i<size;i++)
     187        region2[i] = multiplier1*region1[i] ;
     188    } else {
     189      for (i=0;i<size;i++)
     190        region2[i] = multiplier1*region1[i] + multiplier2*region2[i];
     191    }
     192  }
     193}
     194CoinWorkDouble
     195innerProduct(const CoinWorkDouble * region1, int size, const CoinWorkDouble * region2)
     196{
     197  int i;
     198  CoinWorkDouble value=0.0;
     199  for (i=0;i<size;i++)
     200    value += region1[i]*region2[i];
     201  return value;
     202}
     203void
     204getNorms(const CoinWorkDouble * region, int size, CoinWorkDouble & norm1, CoinWorkDouble & norm2)
     205{
     206  norm1 = 0.0;
     207  norm2 = 0.0;
     208  int i;
     209  for (i=0;i<size;i++) {
     210    norm2 += region[i]*region[i];
     211    norm1 = CoinMax(norm1,CoinAbs(region[i]));
     212  }
     213}
     214#endif
    114215#ifdef DEBUG_MEMORY
    115216#include <malloc.h>
  • trunk/Clp/src/ClpHelperFunctions.hpp

    r754 r1367  
    1818double innerProduct(const double * region1, int size, const double * region2);
    1919void getNorms(const double * region, int size, double & norm1, double & norm2);
     20#if COIN_LONG_WORK
     21  // For long double versions
     22CoinWorkDouble maximumAbsElement(const CoinWorkDouble * region, int size);
     23void setElements(CoinWorkDouble * region, int size, CoinWorkDouble value);
     24void multiplyAdd(const CoinWorkDouble * region1, int size, CoinWorkDouble multiplier1,
     25                 CoinWorkDouble * region2, CoinWorkDouble multiplier2);
     26CoinWorkDouble innerProduct(const CoinWorkDouble * region1, int size, const CoinWorkDouble * region2);
     27void getNorms(const CoinWorkDouble * region, int size, CoinWorkDouble & norm1, CoinWorkDouble & norm2);
     28inline void
     29CoinMemcpyN(const double * from, const int size, CoinWorkDouble * to)
     30{
     31  for (int i=0;i<size;i++)
     32    to[i]=from[i];
     33}
     34inline void
     35CoinMemcpyN(const CoinWorkDouble * from, const int size, double * to)
     36{
     37  for (int i=0;i<size;i++)
     38    to[i]=from[i];
     39}
     40inline CoinWorkDouble
     41CoinMax(const CoinWorkDouble x1, const double x2)
     42{
     43    return (x1 > x2) ? x1 : x2;
     44}
     45inline CoinWorkDouble
     46CoinMax(double x1, const CoinWorkDouble x2)
     47{
     48    return (x1 > x2) ? x1 : x2;
     49}
     50inline CoinWorkDouble
     51CoinMin(const CoinWorkDouble x1, const double x2)
     52{
     53    return (x1 < x2) ? x1 : x2;
     54}
     55inline CoinWorkDouble
     56CoinMin(double x1, const CoinWorkDouble x2)
     57{
     58    return (x1 < x2) ? x1 : x2;
     59}
     60inline CoinWorkDouble CoinSqrt(CoinWorkDouble x)
     61{
     62  return sqrtl(x);
     63}
     64#else
     65inline double CoinSqrt(double x)
     66{
     67  return sqrt(x);
     68}
     69#endif
    2070
    2171/// Following only included if ClpPdco defined
  • trunk/Clp/src/ClpInterior.cpp

    r1321 r1367  
    109109  algorithm_(-1)
    110110{
    111   memset(historyInfeasibility_,0,LENGTH_HISTORY*sizeof(double));
     111  memset(historyInfeasibility_,0,LENGTH_HISTORY*sizeof(CoinWorkDouble));
    112112  solveType_=3; // say interior based life form
    113113  cholesky_ = new ClpCholeskyDense(); // put in placeholder
     
    198198    algorithm_(-1)
    199199{
    200   memset(historyInfeasibility_,0,LENGTH_HISTORY*sizeof(double));
     200  memset(historyInfeasibility_,0,LENGTH_HISTORY*sizeof(CoinWorkDouble));
    201201  solveType_=3; // say interior based life form
    202202  cholesky_= new ClpCholeskyDense();
     
    351351  algorithm_(-1)
    352352{
    353   memset(historyInfeasibility_,0,LENGTH_HISTORY*sizeof(double));
     353  memset(historyInfeasibility_,0,LENGTH_HISTORY*sizeof(CoinWorkDouble));
    354354  solveType_=3; // say interior based life form
    355355  cholesky_= new ClpCholeskyDense();
     
    537537  int nTotal = numberRows_+numberColumns_;
    538538  delete [] solution_;
    539   solution_ = new double[nTotal];
     539  solution_ = new CoinWorkDouble[nTotal];
    540540  CoinMemcpyN(columnActivity_,  numberColumns_,solution_);
    541541  CoinMemcpyN(rowActivity_,     numberRows_,solution_+numberColumns_);
    542542  delete [] cost_;
    543   cost_ = new double[nTotal];
     543  cost_ = new CoinWorkDouble[nTotal];
    544544  int i;
    545   double direction = optimizationDirection_*objectiveScale_;
     545  CoinWorkDouble direction = optimizationDirection_*objectiveScale_;
    546546  // direction is actually scale out not scale in
    547547  if (direction)
     
    550550  for (i=0;i<numberColumns_;i++)
    551551    cost_[i] = direction*obj[i];
    552   memset(cost_+numberColumns_,0,numberRows_*sizeof(double));
     552  memset(cost_+numberColumns_,0,numberRows_*sizeof(CoinWorkDouble));
    553553  // do scaling if needed
    554554  if (scalingFlag_>0&&!rowScale_) {
     
    558558  delete [] lower_;
    559559  delete [] upper_;
    560   lower_ = new double[nTotal];
    561   upper_ = new double[nTotal];
     560  lower_ = new CoinWorkDouble[nTotal];
     561  upper_ = new CoinWorkDouble[nTotal];
    562562  rowLowerWork_ = lower_+numberColumns_;
    563563  columnLowerWork_ = lower_;
     
    587587  if(rowScale_) {
    588588    for (i=0;i<numberColumns_;i++) {
    589       double multiplier = rhsScale_/columnScale_[i];
     589      CoinWorkDouble multiplier = rhsScale_/columnScale_[i];
    590590      cost_[i] *= columnScale_[i];
    591591      if (columnLowerWork_[i]>-1.0e50)
     
    596596    }
    597597    for (i=0;i<numberRows_;i++) {
    598       double multiplier = rhsScale_*rowScale_[i];
     598      CoinWorkDouble multiplier = rhsScale_*rowScale_[i];
    599599      if (rowLowerWork_[i]>-1.0e50)
    600600        rowLowerWork_[i] *= multiplier;
     
    612612  }
    613613  assert (!errorRegion_);
    614   errorRegion_ = new double [numberRows_];
     614  errorRegion_ = new CoinWorkDouble [numberRows_];
    615615  assert (!rhsFixRegion_);
    616   rhsFixRegion_ = new double [numberRows_];
     616  rhsFixRegion_ = new CoinWorkDouble [numberRows_];
    617617  assert (!deltaY_);
    618   deltaY_ = new double [numberRows_];
     618  deltaY_ = new CoinWorkDouble [numberRows_];
    619619  CoinZeroN(deltaY_,numberRows_);
    620620  assert (!upperSlack_);
    621   upperSlack_ = new double [nTotal];
     621  upperSlack_ = new CoinWorkDouble [nTotal];
    622622  assert (!lowerSlack_);
    623   lowerSlack_ = new double [nTotal];
     623  lowerSlack_ = new CoinWorkDouble [nTotal];
    624624  assert (!diagonal_);
    625   diagonal_ = new double [nTotal];
     625  diagonal_ = new CoinWorkDouble [nTotal];
    626626  assert (!deltaX_);
    627   deltaX_ = new double [nTotal];
     627  deltaX_ = new CoinWorkDouble [nTotal];
    628628  CoinZeroN(deltaX_,nTotal);
    629629  assert (!deltaZ_);
    630   deltaZ_ = new double [nTotal];
     630  deltaZ_ = new CoinWorkDouble [nTotal];
    631631  CoinZeroN(deltaZ_,nTotal);
    632632  assert (!deltaW_);
    633   deltaW_ = new double [nTotal];
     633  deltaW_ = new CoinWorkDouble [nTotal];
    634634  CoinZeroN(deltaW_,nTotal);
    635635  assert (!deltaSU_);
    636   deltaSU_ = new double [nTotal];
     636  deltaSU_ = new CoinWorkDouble [nTotal];
    637637  CoinZeroN(deltaSU_,nTotal);
    638638  assert (!deltaSL_);
    639   deltaSL_ = new double [nTotal];
     639  deltaSL_ = new CoinWorkDouble [nTotal];
    640640  CoinZeroN(deltaSL_,nTotal);
    641641  assert (!primalR_);
     
    643643  // create arrays if we are doing KKT
    644644  if (cholesky_->type()>=20) {
    645     primalR_ = new double [nTotal];
     645    primalR_ = new CoinWorkDouble [nTotal];
    646646    CoinZeroN(primalR_,nTotal);
    647     dualR_ = new double [numberRows_];
     647    dualR_ = new CoinWorkDouble [numberRows_];
    648648    CoinZeroN(dualR_,numberRows_);
    649649  }
    650650  assert (!rhsB_);
    651   rhsB_ = new double [numberRows_];
     651  rhsB_ = new CoinWorkDouble [numberRows_];
    652652  CoinZeroN(rhsB_,numberRows_);
    653653  assert (!rhsU_);
    654   rhsU_ = new double [nTotal];
     654  rhsU_ = new CoinWorkDouble [nTotal];
    655655  CoinZeroN(rhsU_,nTotal);
    656656  assert (!rhsL_);
    657   rhsL_ = new double [nTotal];
     657  rhsL_ = new CoinWorkDouble [nTotal];
    658658  CoinZeroN(rhsL_,nTotal);
    659659  assert (!rhsZ_);
    660   rhsZ_ = new double [nTotal];
     660  rhsZ_ = new CoinWorkDouble [nTotal];
    661661  CoinZeroN(rhsZ_,nTotal);
    662662  assert (!rhsW_);
    663   rhsW_ = new double [nTotal];
     663  rhsW_ = new CoinWorkDouble [nTotal];
    664664  CoinZeroN(rhsW_,nTotal);
    665665  assert (!rhsC_);
    666   rhsC_ = new double [nTotal];
     666  rhsC_ = new CoinWorkDouble [nTotal];
    667667  CoinZeroN(rhsC_,nTotal);
    668668  assert (!workArray_);
    669   workArray_ = new double [nTotal];
     669  workArray_ = new CoinWorkDouble [nTotal];
    670670  CoinZeroN(workArray_,nTotal);
    671671  assert (!zVec_);
    672   zVec_ = new double [nTotal];
     672  zVec_ = new CoinWorkDouble [nTotal];
    673673  CoinZeroN(zVec_,nTotal);
    674674  assert (!wVec_);
    675   wVec_ = new double [nTotal];
     675  wVec_ = new CoinWorkDouble [nTotal];
    676676  CoinZeroN(wVec_,nTotal);
    677677  assert (!dj_);
    678   dj_ = new double [nTotal];
     678  dj_ = new CoinWorkDouble [nTotal];
    679679  if (!status_)
    680680    status_ = new unsigned char [numberRows_+numberColumns_];
     
    687687  int i;
    688688  if (optimizationDirection_!=1.0||objectiveScale_!=1.0) {
    689     double scaleC = optimizationDirection_/objectiveScale_;
     689    CoinWorkDouble scaleC = optimizationDirection_/objectiveScale_;
    690690    // and modify all dual signs
    691691    for (i=0;i<numberColumns_;i++)
     
    695695  }
    696696  if (rowScale_) {
    697     double scaleR = 1.0/rhsScale_;
     697    CoinWorkDouble scaleR = 1.0/rhsScale_;
    698698    for (i=0;i<numberColumns_;i++) {
    699       double scaleFactor = columnScale_[i];
    700       double valueScaled = columnActivity_[i];
     699      CoinWorkDouble scaleFactor = columnScale_[i];
     700      CoinWorkDouble valueScaled = columnActivity_[i];
    701701      columnActivity_[i] = valueScaled*scaleFactor*scaleR;
    702       double valueScaledDual = reducedCost_[i];
     702      CoinWorkDouble valueScaledDual = reducedCost_[i];
    703703      reducedCost_[i] = valueScaledDual/scaleFactor;
    704704    }
    705705    for (i=0;i<numberRows_;i++) {
    706       double scaleFactor = rowScale_[i];
    707       double valueScaled = rowActivity_[i];
     706      CoinWorkDouble scaleFactor = rowScale_[i];
     707      CoinWorkDouble valueScaled = rowActivity_[i];
    708708      rowActivity_[i] = (valueScaled*scaleR)/scaleFactor;
    709       double valueScaledDual = dual_[i];
     709      CoinWorkDouble valueScaledDual = dual_[i];
    710710      dual_[i] = valueScaledDual*scaleFactor;
    711711    }
    712712  } else if (rhsScale_!=1.0) {
    713     double scaleR = 1.0/rhsScale_;
     713    CoinWorkDouble scaleR = 1.0/rhsScale_;
    714714    for (i=0;i<numberColumns_;i++) {
    715       double valueScaled = columnActivity_[i];
     715      CoinWorkDouble valueScaled = columnActivity_[i];
    716716      columnActivity_[i] = valueScaled*scaleR;
    717717    }
    718718    for (i=0;i<numberRows_;i++) {
    719       double valueScaled = rowActivity_[i];
     719      CoinWorkDouble valueScaled = rowActivity_[i];
    720720      rowActivity_[i] = valueScaled*scaleR;
    721721    }
     
    762762  }
    763763  int numberBad ;
    764   double largestBound, smallestBound, minimumGap;
    765   double smallestObj, largestObj;
     764  CoinWorkDouble largestBound, smallestBound, minimumGap;
     765  CoinWorkDouble smallestObj, largestObj;
    766766  int firstBad;
    767767  int modifiedBounds=0;
     
    775775  largestObj=0.0;
    776776  // If bounds are too close - fix
    777   double fixTolerance = 1.1*primalTolerance();
     777  CoinWorkDouble fixTolerance = 1.1*primalTolerance();
    778778  for (i=numberColumns_;i<numberColumns_+numberRows_;i++) {
    779     double value;
    780     value = fabs(cost_[i]);
     779    CoinWorkDouble value;
     780    value = CoinAbs(cost_[i]);
    781781    if (value>1.0e50) {
    782782      numberBad++;
     
    805805    }
    806806    if (lower_[i]>-1.0e100&&lower_[i]) {
    807       value = fabs(lower_[i]);
     807      value = CoinAbs(lower_[i]);
    808808      if (value>largestBound)
    809809        largestBound=value;
     
    812812    }
    813813    if (upper_[i]<1.0e100&&upper_[i]) {
    814       value = fabs(upper_[i]);
     814      value = CoinAbs(upper_[i]);
    815815      if (value>largestBound)
    816816        largestBound=value;
     
    821821  if (largestBound)
    822822    handler_->message(CLP_RIMSTATISTICS3,messages_)
    823       <<smallestBound
    824       <<largestBound
    825       <<minimumGap
     823      <<static_cast<double>(smallestBound)
     824      <<static_cast<double>(largestBound)
     825      <<static_cast<double>(minimumGap)
    826826      <<CoinMessageEol;
    827827  minimumGap=1.0e100;
     
    829829  largestBound=0.0;
    830830  for (i=0;i<numberColumns_;i++) {
    831     double value;
    832     value = fabs(cost_[i]);
     831    CoinWorkDouble value;
     832    value = CoinAbs(cost_[i]);
    833833    if (value>1.0e50) {
    834834      numberBad++;
     
    857857    }
    858858    if (lower_[i]>-1.0e100&&lower_[i]) {
    859       value = fabs(lower_[i]);
     859      value = CoinAbs(lower_[i]);
    860860      if (value>largestBound)
    861861        largestBound=value;
     
    864864    }
    865865    if (upper_[i]<1.0e100&&upper_[i]) {
    866       value = fabs(upper_[i]);
     866      value = CoinAbs(upper_[i]);
    867867      if (value>largestBound)
    868868        largestBound=value;
     
    885885      <<CoinMessageEol;
    886886  handler_->message(CLP_RIMSTATISTICS1,messages_)
    887     <<smallestObj
    888     <<largestObj
     887    <<static_cast<double>(smallestObj)
     888    <<static_cast<double>(largestObj)
    889889    <<CoinMessageEol;  if (largestBound)
    890890    handler_->message(CLP_RIMSTATISTICS2,messages_)
    891       <<smallestBound
    892       <<largestBound
    893       <<minimumGap
     891      <<static_cast<double>(smallestBound)
     892      <<static_cast<double>(largestBound)
     893      <<static_cast<double>(minimumGap)
    894894      <<CoinMessageEol;
    895895  return true;
     
    991991{
    992992  int iRow,iColumn;
    993   CoinMemcpyN(cost_,numberColumns_,reducedCost_);
    994   matrix_->transposeTimes(-1.0,dual_,reducedCost_);
     993  CoinWorkDouble * reducedCost = reinterpret_cast<CoinWorkDouble *>(reducedCost_);
     994  CoinWorkDouble * dual = reinterpret_cast<CoinWorkDouble *>(dual_);
     995  CoinMemcpyN(cost_,numberColumns_,reducedCost);
     996  matrix_->transposeTimes(-1.0,dual,reducedCost);
    995997  // Now modify reduced costs for quadratic
    996   double quadraticOffset=quadraticDjs(reducedCost_,
     998  CoinWorkDouble quadraticOffset=quadraticDjs(reducedCost,
    997999                                      solution_,scaleFactor_);
    9981000
     
    10011003  sumPrimalInfeasibilities_=0.0;
    10021004  sumDualInfeasibilities_=0.0;
    1003   double dualTolerance =  10.0*dblParam_[ClpDualTolerance];
    1004   double primalTolerance =  dblParam_[ClpPrimalTolerance];
    1005   double primalTolerance2 =  10.0*dblParam_[ClpPrimalTolerance];
     1005  CoinWorkDouble dualTolerance =  10.0*dblParam_[ClpDualTolerance];
     1006  CoinWorkDouble primalTolerance =  dblParam_[ClpPrimalTolerance];
     1007  CoinWorkDouble primalTolerance2 =  10.0*dblParam_[ClpPrimalTolerance];
    10061008  worstComplementarity_=0.0;
    10071009  complementarityGap_=0.0;
     
    10091011  // Done scaled - use permanent regions for output
    10101012  // but internal for bounds
    1011   const double * lower = lower_+numberColumns_;
    1012   const double * upper = upper_+numberColumns_;
     1013  const CoinWorkDouble * lower = lower_+numberColumns_;
     1014  const CoinWorkDouble * upper = upper_+numberColumns_;
    10131015  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);
     1016    CoinWorkDouble infeasibility=0.0;
     1017    CoinWorkDouble distanceUp = CoinMin(upper[iRow]-
     1018                                        rowActivity_[iRow],static_cast<CoinWorkDouble>(1.0e10));
     1019    CoinWorkDouble distanceDown = CoinMin(rowActivity_[iRow] -
     1020                                          lower[iRow],static_cast<CoinWorkDouble>(1.0e10));
    10191021    if (distanceUp>primalTolerance2) {
    1020       double value = dual_[iRow];
     1022      CoinWorkDouble value = dual[iRow];
    10211023      // should not be negative
    10221024      if (value<-dualTolerance) {
     
    10291031    }
    10301032    if (distanceDown>primalTolerance2) {
    1031       double value = dual_[iRow];
     1033      CoinWorkDouble value = dual[iRow];
    10321034      // should not be positive
    10331035      if (value>dualTolerance) {
     
    10511053  upper = upper_;
    10521054  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    1053     double infeasibility=0.0;
     1055    CoinWorkDouble infeasibility=0.0;
    10541056    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);
     1057    CoinWorkDouble distanceUp = CoinMin(upper[iColumn]-
     1058                                        columnActivity_[iColumn],static_cast<CoinWorkDouble>(1.0e10));
     1059    CoinWorkDouble distanceDown = CoinMin(columnActivity_[iColumn] -
     1060                                          lower[iColumn],static_cast<CoinWorkDouble>(1.0e10));
    10591061    if (distanceUp>primalTolerance2) {
    1060       double value = reducedCost_[iColumn];
     1062      CoinWorkDouble value = reducedCost[iColumn];
    10611063      // should not be negative
    10621064      if (value<-dualTolerance) {
     
    10691071    }
    10701072    if (distanceDown>primalTolerance2) {
    1071       double value = reducedCost_[iColumn];
     1073      CoinWorkDouble value = reducedCost[iColumn];
    10721074      // should not be positive
    10731075      if (value>dualTolerance) {
     
    10881090    }
    10891091  }
     1092#if COIN_LONG_WORK
     1093  // ok as packs down
     1094  CoinMemcpyN(reducedCost,numberColumns_,reducedCost_);
     1095#endif
    10901096  // add in offset
    10911097  objectiveValue_ += 0.5*quadraticOffset;
     
    11421148{
    11431149  // Arrays for change in columns and rhs
    1144   double * columnChange = new double[numberColumns_];
    1145   double * rowChange = new double[numberRows_];
     1150  CoinWorkDouble * columnChange = new CoinWorkDouble[numberColumns_];
     1151  CoinWorkDouble * rowChange = new CoinWorkDouble[numberRows_];
    11461152  CoinZeroN(columnChange,numberColumns_);
    11471153  CoinZeroN(rowChange,numberRows_);
    11481154  matrix_->times(1.0,columnChange,rowChange);
    11491155  int i;
    1150   double tolerance = primalTolerance();
     1156  CoinWorkDouble tolerance = primalTolerance();
    11511157  for (i=0;i<numberColumns_;i++) {
    11521158    if (columnUpper_[i]<1.0e20||columnLower_[i]>-1.0e20) {
     
    11541160        if (fixedOrFree(i)) {
    11551161          if (columnActivity_[i]-columnLower_[i]<columnUpper_[i]-columnActivity_[i]) {
    1156             double change = columnLower_[i]-columnActivity_[i];
    1157             if (fabs(change)<tolerance) {
     1162            CoinWorkDouble change = columnLower_[i]-columnActivity_[i];
     1163            if (CoinAbs(change)<tolerance) {
    11581164              if (reallyFix)
    11591165                columnUpper_[i]=columnLower_[i];
     
    11621168            }
    11631169          } else {
    1164             double change = columnUpper_[i]-columnActivity_[i];
    1165             if (fabs(change)<tolerance) {
     1170            CoinWorkDouble change = columnUpper_[i]-columnActivity_[i];
     1171            if (CoinAbs(change)<tolerance) {
    11661172              if (reallyFix)
    11671173                columnLower_[i]=columnUpper_[i];
     
    11771183  matrix_->times(1.0,columnChange,rowChange);
    11781184  // If makes mess of things then don't do
    1179   double newSum=0.0;
     1185  CoinWorkDouble newSum=0.0;
    11801186  for (i=0;i<numberRows_;i++) {
    1181     double value=rowActivity_[i]+rowChange[i];
     1187    CoinWorkDouble value=rowActivity_[i]+rowChange[i];
    11821188    if (value>rowUpper_[i]+tolerance)
    11831189      newSum += value-rowUpper_[i]-tolerance;
     
    11981204            if (fixedOrFree(i+numberColumns_)) {
    11991205              if (rowActivity_[i]-rowLower_[i]<rowUpper_[i]-rowActivity_[i]) {
    1200                 double change = rowLower_[i]-rowActivity_[i];
    1201                 if (fabs(change)<tolerance) {
     1206                CoinWorkDouble change = rowLower_[i]-rowActivity_[i];
     1207                if (CoinAbs(change)<tolerance) {
    12021208                  if (reallyFix)
    12031209                    rowUpper_[i]=rowLower_[i];
     
    12051211                }
    12061212              } else {
    1207                 double change = rowLower_[i]-rowActivity_[i];
    1208                 if (fabs(change)<tolerance) {
     1213                CoinWorkDouble change = rowLower_[i]-rowActivity_[i];
     1214                if (CoinAbs(change)<tolerance) {
    12091215                  if (reallyFix)
    12101216                    rowLower_[i]=rowUpper_[i];
     
    12231229/* Modifies djs to allow for quadratic.
    12241230   returns quadratic offset */
    1225 double
    1226 ClpInterior::quadraticDjs(double * djRegion, const double * solution,
    1227                           double scaleFactor)
    1228 {
    1229   double quadraticOffset=0.0;
     1231CoinWorkDouble
     1232ClpInterior::quadraticDjs(CoinWorkDouble * djRegion, const CoinWorkDouble * solution,
     1233                          CoinWorkDouble scaleFactor)
     1234{
     1235  CoinWorkDouble quadraticOffset=0.0;
    12301236#ifndef NO_RTTI
    12311237  ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
     
    12431249    int numberColumns = quadratic->getNumCols();
    12441250    for (int iColumn=0;iColumn<numberColumns;iColumn++) {
    1245       double value = 0.0;
     1251      CoinWorkDouble value = 0.0;
    12461252      for (CoinBigIndex j=columnQuadraticStart[iColumn];
    12471253           j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
    12481254        int jColumn = columnQuadratic[j];
    1249         double valueJ = solution[jColumn];
    1250         double elementValue = quadraticElement[j];
     1255        CoinWorkDouble valueJ = solution[jColumn];
     1256        CoinWorkDouble elementValue = quadraticElement[j];
    12511257        //value += valueI*valueJ*elementValue;
    12521258        value += valueJ*elementValue;
  • trunk/Clp/src/ClpInterior.hpp

    r1264 r1367  
    174174  {algorithm_=value; }
    175175  /// Sum of dual infeasibilities
    176   inline double sumDualInfeasibilities() const
     176  inline CoinWorkDouble sumDualInfeasibilities() const
    177177          { return sumDualInfeasibilities_;}
    178178  /// Sum of primal infeasibilities
    179   inline double sumPrimalInfeasibilities() const
     179  inline CoinWorkDouble sumPrimalInfeasibilities() const
    180180          { return sumPrimalInfeasibilities_;}
    181181  /// dualObjective.
    182   inline double dualObjective() const
     182  inline CoinWorkDouble dualObjective() const
    183183  { return dualObjective_;}
    184184  /// primalObjective.
    185   inline double primalObjective() const
     185  inline CoinWorkDouble primalObjective() const
    186186  { return primalObjective_;}
    187187  /// diagonalNorm
    188   inline double diagonalNorm() const
     188  inline CoinWorkDouble diagonalNorm() const
    189189  { return diagonalNorm_;}
    190190  /// linearPerturbation
    191   inline double linearPerturbation() const
     191  inline CoinWorkDouble linearPerturbation() const
    192192  { return linearPerturbation_;}
    193   inline void setLinearPerturbation(double value)
     193  inline void setLinearPerturbation(CoinWorkDouble value)
    194194  { linearPerturbation_=value;}
     195  /// projectionTolerance
     196  inline CoinWorkDouble projectionTolerance() const
     197  { return projectionTolerance_;}
     198  inline void setProjectionTolerance(CoinWorkDouble value)
     199  { projectionTolerance_=value;}
    195200  /// diagonalPerturbation
    196   inline double diagonalPerturbation() const
     201  inline CoinWorkDouble diagonalPerturbation() const
    197202  { return diagonalPerturbation_;}
    198   inline void setDiagonalPerturbation(double value)
     203  inline void setDiagonalPerturbation(CoinWorkDouble value)
    199204  { diagonalPerturbation_=value;}
    200205  /// gamma
    201   inline double gamma() const
     206  inline CoinWorkDouble gamma() const
    202207  { return gamma_;}
    203   inline void setGamma(double value)
     208  inline void setGamma(CoinWorkDouble value)
    204209  { gamma_=value;}
    205210  /// delta
    206   inline double delta() const
     211  inline CoinWorkDouble delta() const
    207212  { return delta_;}
    208   inline void setDelta(double value)
     213  inline void setDelta(CoinWorkDouble value)
    209214  { delta_=value;}
    210215  /// ComplementarityGap
    211   inline double complementarityGap() const
     216  inline CoinWorkDouble complementarityGap() const
    212217          { return complementarityGap_;}
    213218  //@}
     
    216221  //@{
    217222  /// Largest error on Ax-b
    218   inline double largestPrimalError() const
     223  inline CoinWorkDouble largestPrimalError() const
    219224          { return largestPrimalError_;}
    220225  /// Largest error on basic duals
    221   inline double largestDualError() const
     226  inline CoinWorkDouble largestDualError() const
    222227          { return largestDualError_;}
    223228  /// Maximum iterations
     
    234239  void fixFixed(bool reallyFix=true);
    235240  /// Primal erturbation vector
    236   inline double * primalR() const
     241  inline CoinWorkDouble * primalR() const
    237242  { return primalR_;}
    238243  /// Dual erturbation vector
    239   inline double * dualR() const
     244  inline CoinWorkDouble * dualR() const
    240245  { return dualR_;}
    241246  //@}
     
    260265  //@{
    261266  /// Raw objective value (so always minimize)
    262   inline double rawObjectiveValue() const
     267  inline CoinWorkDouble rawObjectiveValue() const
    263268  { return objectiveValue_;}
    264269  /// Returns 1 if sequence indicates column
     
    272277  /** Modifies djs to allow for quadratic.
    273278      returns quadratic offset */
    274   double quadraticDjs(double * djRegion, const double * solution,
    275                       double scaleFactor);
     279  CoinWorkDouble quadraticDjs(CoinWorkDouble * djRegion, const CoinWorkDouble * solution,
     280                      CoinWorkDouble scaleFactor);
    276281
    277282  /// To say a variable is fixed
     
    370375  //@{
    371376  /// Largest error on Ax-b
    372   double largestPrimalError_;
     377  CoinWorkDouble largestPrimalError_;
    373378  /// Largest error on basic duals
    374   double largestDualError_;
     379  CoinWorkDouble largestDualError_;
    375380  /// Sum of dual infeasibilities
    376   double sumDualInfeasibilities_;
     381  CoinWorkDouble sumDualInfeasibilities_;
    377382  /// Sum of primal infeasibilities
    378   double sumPrimalInfeasibilities_;
     383  CoinWorkDouble sumPrimalInfeasibilities_;
    379384  /// Worst complementarity
    380   double worstComplementarity_;
     385  CoinWorkDouble worstComplementarity_;
    381386  ///
    382387public:
    383   double xsize_;
    384   double zsize_;
     388  CoinWorkDouble xsize_;
     389  CoinWorkDouble zsize_;
    385390protected:
    386391  /// Working copy of lower bounds (Owner of arrays below)
    387   double * lower_;
     392  CoinWorkDouble * lower_;
    388393  /// Row lower bounds - working copy
    389   double * rowLowerWork_;
     394  CoinWorkDouble * rowLowerWork_;
    390395  /// Column lower bounds - working copy
    391   double * columnLowerWork_;
     396  CoinWorkDouble * columnLowerWork_;
    392397  /// Working copy of upper bounds (Owner of arrays below)
    393   double * upper_;
     398  CoinWorkDouble * upper_;
    394399  /// Row upper bounds - working copy
    395   double * rowUpperWork_;
     400  CoinWorkDouble * rowUpperWork_;
    396401  /// Column upper bounds - working copy
    397   double * columnUpperWork_;
     402  CoinWorkDouble * columnUpperWork_;
    398403  /// Working copy of objective
    399   double * cost_;
     404  CoinWorkDouble * cost_;
    400405public:
    401406  /// Rhs
    402   double * rhs_;
    403   double * x_;
    404   double * y_;
    405   double * dj_;
     407  CoinWorkDouble * rhs_;
     408  CoinWorkDouble * x_;
     409  CoinWorkDouble * y_;
     410  CoinWorkDouble * dj_;
    406411protected:
    407412  /// Pointer to Lsqr object
     
    411416  /// Below here is standard barrier stuff
    412417  /// mu.
    413   double mu_;
     418  CoinWorkDouble mu_;
    414419  /// objectiveNorm.
    415   double objectiveNorm_;
     420  CoinWorkDouble objectiveNorm_;
    416421  /// rhsNorm.
    417   double rhsNorm_;
     422  CoinWorkDouble rhsNorm_;
    418423  /// solutionNorm.
    419   double solutionNorm_;
     424  CoinWorkDouble solutionNorm_;
    420425  /// dualObjective.
    421   double dualObjective_;
     426  CoinWorkDouble dualObjective_;
    422427  /// primalObjective.
    423   double primalObjective_;
     428  CoinWorkDouble primalObjective_;
    424429  /// diagonalNorm.
    425   double diagonalNorm_;
     430  CoinWorkDouble diagonalNorm_;
    426431  /// stepLength
    427   double stepLength_;
     432  CoinWorkDouble stepLength_;
    428433  /// linearPerturbation
    429   double linearPerturbation_;
     434  CoinWorkDouble linearPerturbation_;
    430435  /// diagonalPerturbation
    431   double diagonalPerturbation_;
     436  CoinWorkDouble diagonalPerturbation_;
    432437  // gamma from Saunders and Tomlin regularized
    433   double gamma_;
     438  CoinWorkDouble gamma_;
    434439  // delta from Saunders and Tomlin regularized
    435   double delta_;
     440  CoinWorkDouble delta_;
    436441  /// targetGap
    437   double targetGap_;
     442  CoinWorkDouble targetGap_;
    438443  /// projectionTolerance
    439   double projectionTolerance_;
     444  CoinWorkDouble projectionTolerance_;
    440445  /// maximumRHSError.  maximum Ax
    441   double maximumRHSError_;
     446  CoinWorkDouble maximumRHSError_;
    442447  /// maximumBoundInfeasibility.
    443   double maximumBoundInfeasibility_;
     448  CoinWorkDouble maximumBoundInfeasibility_;
    444449  /// maximumDualError.
    445   double maximumDualError_;
     450  CoinWorkDouble maximumDualError_;
    446451  /// diagonalScaleFactor.
    447   double diagonalScaleFactor_;
     452  CoinWorkDouble diagonalScaleFactor_;
    448453  /// scaleFactor.  For scaling objective
    449   double scaleFactor_;
     454  CoinWorkDouble scaleFactor_;
    450455  /// actualPrimalStep
    451   double actualPrimalStep_;
     456  CoinWorkDouble actualPrimalStep_;
    452457  /// actualDualStep
    453   double actualDualStep_;
     458  CoinWorkDouble actualDualStep_;
    454459  /// smallestInfeasibility
    455   double smallestInfeasibility_;
     460  CoinWorkDouble smallestInfeasibility_;
    456461  /// historyInfeasibility.
    457462#define LENGTH_HISTORY 5
    458   double historyInfeasibility_[LENGTH_HISTORY];
     463  CoinWorkDouble historyInfeasibility_[LENGTH_HISTORY];
    459464  /// complementarityGap.
    460   double complementarityGap_;
     465  CoinWorkDouble complementarityGap_;
    461466  /// baseObjectiveNorm
    462   double baseObjectiveNorm_;
     467  CoinWorkDouble baseObjectiveNorm_;
    463468  /// worstDirectionAccuracy
    464   double worstDirectionAccuracy_;
     469  CoinWorkDouble worstDirectionAccuracy_;
    465470  /// maximumRHSChange
    466   double maximumRHSChange_;
     471  CoinWorkDouble maximumRHSChange_;
    467472  /// errorRegion. i.e. Ax
    468   double * errorRegion_;
     473  CoinWorkDouble * errorRegion_;
    469474  /// rhsFixRegion.
    470   double * rhsFixRegion_;
     475  CoinWorkDouble * rhsFixRegion_;
    471476  /// upperSlack
    472   double * upperSlack_;
     477  CoinWorkDouble * upperSlack_;
    473478  /// lowerSlack
    474   double * lowerSlack_;
     479  CoinWorkDouble * lowerSlack_;
    475480  /// diagonal
    476   double * diagonal_;
     481  CoinWorkDouble * diagonal_;
    477482  /// solution
    478   double * solution_;
     483  CoinWorkDouble * solution_;
    479484  /// work array
    480   double * workArray_;
     485  CoinWorkDouble * workArray_;
    481486  /// delta X
    482   double * deltaX_;
     487  CoinWorkDouble * deltaX_;
    483488  /// delta Y
    484   double * deltaY_;
     489  CoinWorkDouble * deltaY_;
    485490  /// deltaZ.
    486   double * deltaZ_;
     491  CoinWorkDouble * deltaZ_;
    487492  /// deltaW.
    488   double * deltaW_;
     493  CoinWorkDouble * deltaW_;
    489494  /// deltaS.
    490   double * deltaSU_;
    491   double * deltaSL_;
     495  CoinWorkDouble * deltaSU_;
     496  CoinWorkDouble * deltaSL_;
    492497  /// Primal regularization array
    493   double * primalR_;
     498  CoinWorkDouble * primalR_;
    494499  /// Dual regularization array
    495   double * dualR_;
     500  CoinWorkDouble * dualR_;
    496501  /// rhs B
    497   double * rhsB_;
     502  CoinWorkDouble * rhsB_;
    498503  /// rhsU.
    499   double * rhsU_;
     504  CoinWorkDouble * rhsU_;
    500505  /// rhsL.
    501   double * rhsL_;
     506  CoinWorkDouble * rhsL_;
    502507  /// rhsZ.
    503   double * rhsZ_;
     508  CoinWorkDouble * rhsZ_;
    504509  /// rhsW.
    505   double * rhsW_;
     510  CoinWorkDouble * rhsW_;
    506511  /// rhs C
    507   double * rhsC_;
     512  CoinWorkDouble * rhsC_;
    508513  /// zVec
    509   double * zVec_;
     514  CoinWorkDouble * zVec_;
    510515  /// wVec
    511   double * wVec_;
     516  CoinWorkDouble * wVec_;
    512517  /// cholesky.
    513518  ClpCholeskyBase * cholesky_;
  • trunk/Clp/src/ClpMatrixBase.cpp

    r1034 r1367  
    605605  abort();
    606606}
     607#if COIN_LONG_WORK
     608// For long double versions (aborts if not supported)
     609void
     610ClpMatrixBase::times(CoinWorkDouble scalar,
     611      const CoinWorkDouble * x, CoinWorkDouble * y) const
     612{
     613  std::cerr<<"long times not supported - ClpMatrixBase"<<std::endl;
     614  abort();
     615}
     616void
     617ClpMatrixBase::transposeTimes(CoinWorkDouble scalar,
     618                              const CoinWorkDouble * x, CoinWorkDouble * y) const
     619{
     620  std::cerr<<"long transposeTimes not supported - ClpMatrixBase"<<std::endl;
     621  abort();
     622}
     623#endif
  • trunk/Clp/src/ClpMatrixBase.hpp

    r1302 r1367  
    55
    66#include "CoinPragma.hpp"
     7#include "CoinFinite.hpp"
    78
    89#include "CoinPackedMatrix.hpp"
     
    277278                              const double * columnScale,
    278279                              double * spare=NULL) const;
     280#if COIN_LONG_WORK
     281  // For long double versions (aborts if not supported)
     282  virtual void times(CoinWorkDouble scalar,
     283                     const CoinWorkDouble * x, CoinWorkDouble * y) const ;
     284  virtual void transposeTimes(CoinWorkDouble scalar,
     285                              const CoinWorkDouble * x, CoinWorkDouble * y) const ;
     286#endif
    279287  /** Return <code>x * scalar *A + y</code> in <code>z</code>.
    280288      Can use y as temporary array (will be empty at end)
  • trunk/Clp/src/ClpMessage.cpp

    r1034 r1367  
    8282  {CLP_BARRIER_INFO,45,3,"Detail - %s"},
    8383  {CLP_BARRIER_END,46,1,"At end primal/dual infeasibilities %g/%g, complementarity gap %g, objective %g"},
    84   {CLP_BARRIER_ACCURACY,47,2,"Relative error in phase %d, refinement %d is %g"},
     84  {CLP_BARRIER_ACCURACY,47,2,"Relative error in phase %d, %d refinements reduced from %g to %g"},
    8585  {CLP_BARRIER_SAFE,48,2,"Initial safe primal value %g, objective norm %g"},
    8686  {CLP_BARRIER_NEGATIVE_GAPS,49,3,"%d negative gaps summing to %g"},
  • trunk/Clp/src/ClpPackedMatrix.cpp

    r1344 r1367  
    59615961  output->setPackedMode(true);
    59625962}
     5963#if COIN_LONG_WORK
     5964// For long double versions
     5965void
     5966ClpPackedMatrix::times(CoinWorkDouble scalar,
     5967                   const CoinWorkDouble * x, CoinWorkDouble * y) const
     5968{
     5969  int iRow,iColumn;
     5970  // get matrix data pointers
     5971  const int * row = matrix_->getIndices();
     5972  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
     5973  const double * elementByColumn = matrix_->getElements();
     5974  //memset(y,0,matrix_->getNumRows()*sizeof(double));
     5975  assert (((flags_&2)!=0)==matrix_->hasGaps());
     5976  if (!(flags_&2)) {
     5977    for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
     5978      CoinBigIndex j;
     5979      CoinWorkDouble value = x[iColumn];
     5980      if (value) {
     5981        CoinBigIndex start = columnStart[iColumn];
     5982        CoinBigIndex end = columnStart[iColumn+1];
     5983        value *= scalar;
     5984        for (j=start; j<end; j++) {
     5985          iRow=row[j];
     5986          y[iRow] += value*elementByColumn[j];
     5987        }
     5988      }
     5989    }
     5990  } else {
     5991    const int * columnLength = matrix_->getVectorLengths();
     5992    for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
     5993      CoinBigIndex j;
     5994      CoinWorkDouble value = x[iColumn];
     5995      if (value) {
     5996        CoinBigIndex start = columnStart[iColumn];
     5997        CoinBigIndex end = start + columnLength[iColumn];
     5998        value *= scalar;
     5999        for (j=start; j<end; j++) {
     6000          iRow=row[j];
     6001          y[iRow] += value*elementByColumn[j];
     6002        }
     6003      }
     6004    }
     6005  }
     6006}
     6007void
     6008ClpPackedMatrix::transposeTimes(CoinWorkDouble scalar,
     6009                                const CoinWorkDouble * x, CoinWorkDouble * y) const
     6010{
     6011  int iColumn;
     6012  // get matrix data pointers
     6013  const int * row = matrix_->getIndices();
     6014  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
     6015  const double * elementByColumn = matrix_->getElements();
     6016  if (!(flags_&2)) {
     6017    if (scalar==-1.0) {
     6018      CoinBigIndex start=columnStart[0];
     6019      for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
     6020        CoinBigIndex j;
     6021        CoinBigIndex next=columnStart[iColumn+1];
     6022        CoinWorkDouble value=y[iColumn];
     6023        for (j=start;j<next;j++) {
     6024          int jRow=row[j];
     6025          value -= x[jRow]*elementByColumn[j];
     6026        }
     6027        start=next;
     6028        y[iColumn] = value;
     6029      }
     6030    } else {
     6031      CoinBigIndex start=columnStart[0];
     6032      for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
     6033        CoinBigIndex j;
     6034        CoinBigIndex next=columnStart[iColumn+1];
     6035        CoinWorkDouble value=0.0;
     6036        for (j=start;j<next;j++) {
     6037          int jRow=row[j];
     6038          value += x[jRow]*elementByColumn[j];
     6039        }
     6040        start=next;
     6041        y[iColumn] += value*scalar;
     6042      }
     6043    }
     6044  } else {
     6045    const int * columnLength = matrix_->getVectorLengths();
     6046    for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
     6047      CoinBigIndex j;
     6048      CoinWorkDouble value=0.0;
     6049      CoinBigIndex start = columnStart[iColumn];
     6050      CoinBigIndex end = start + columnLength[iColumn];
     6051      for (j=start; j<end; j++) {
     6052        int jRow=row[j];
     6053        value += x[jRow]*elementByColumn[j];
     6054      }
     6055      y[iColumn] += value*scalar;
     6056    }
     6057  }
     6058}
     6059#endif
    59636060#ifdef CLP_ALL_ONE_FILE
    59646061#undef reference
  • trunk/Clp/src/ClpPackedMatrix.hpp

    r1344 r1367  
    264264  /// Sets up an effective RHS
    265265  void useEffectiveRhs(ClpSimplex * model);
     266#if COIN_LONG_WORK
     267  // For long double versions
     268  virtual void times(CoinWorkDouble scalar,
     269                     const CoinWorkDouble * x, CoinWorkDouble * y) const ;
     270  virtual void transposeTimes(CoinWorkDouble scalar,
     271                              const CoinWorkDouble * x, CoinWorkDouble * y) const ;
     272#endif
    266273//@}
    267274
  • trunk/Clp/src/ClpPredictorCorrector.cpp

    r1366 r1367  
    3333    fwrite(&numberRows,sizeof(int),1,fp);
    3434    fwrite(&numberColumns,sizeof(int),1,fp);
    35     double dsave[20];
     35    CoinWorkDouble dsave[20];
    3636    memset(dsave,0,sizeof(dsave));
    37     fwrite(dsave,sizeof(double),20,fp);
     37    fwrite(dsave,sizeof(CoinWorkDouble),20,fp);
    3838    int msave[20];
    3939    memset(msave,0,sizeof(msave));
    4040    msave[0]=numberIterations_;
    4141    fwrite(msave,sizeof(int),20,fp);
    42     fwrite(dual_,sizeof(double),numberRows,fp);
    43     fwrite(errorRegion_,sizeof(double),numberRows,fp);
    44     fwrite(rhsFixRegion_,sizeof(double),numberRows,fp);
    45     fwrite(solution_,sizeof(double),numberColumns,fp);
    46     fwrite(solution_+numberColumns,sizeof(double),numberRows,fp);
    47     fwrite(diagonal_,sizeof(double),numberColumns,fp);
    48     fwrite(diagonal_+numberColumns,sizeof(double),numberRows,fp);
    49     fwrite(wVec_,sizeof(double),numberColumns,fp);
    50     fwrite(wVec_+numberColumns,sizeof(double),numberRows,fp);
    51     fwrite(zVec_,sizeof(double),numberColumns,fp);
    52     fwrite(zVec_+numberColumns,sizeof(double),numberRows,fp);
    53     fwrite(upperSlack_,sizeof(double),numberColumns,fp);
    54     fwrite(upperSlack_+numberColumns,sizeof(double),numberRows,fp);
    55     fwrite(lowerSlack_,sizeof(double),numberColumns,fp);
    56     fwrite(lowerSlack_+numberColumns,sizeof(double),numberRows,fp);
     42    fwrite(dual_,sizeof(CoinWorkDouble),numberRows,fp);
     43    fwrite(errorRegion_,sizeof(CoinWorkDouble),numberRows,fp);
     44    fwrite(rhsFixRegion_,sizeof(CoinWorkDouble),numberRows,fp);
     45    fwrite(solution_,sizeof(CoinWorkDouble),numberColumns,fp);
     46    fwrite(solution_+numberColumns,sizeof(CoinWorkDouble),numberRows,fp);
     47    fwrite(diagonal_,sizeof(CoinWorkDouble),numberColumns,fp);
     48    fwrite(diagonal_+numberColumns,sizeof(CoinWorkDouble),numberRows,fp);
     49    fwrite(wVec_,sizeof(CoinWorkDouble),numberColumns,fp);
     50    fwrite(wVec_+numberColumns,sizeof(CoinWorkDouble),numberRows,fp);
     51    fwrite(zVec_,sizeof(CoinWorkDouble),numberColumns,fp);
     52    fwrite(zVec_+numberColumns,sizeof(CoinWorkDouble),numberRows,fp);
     53    fwrite(upperSlack_,sizeof(CoinWorkDouble),numberColumns,fp);
     54    fwrite(upperSlack_+numberColumns,sizeof(CoinWorkDouble),numberRows,fp);
     55    fwrite(lowerSlack_,sizeof(CoinWorkDouble),numberColumns,fp);
     56    fwrite(lowerSlack_+numberColumns,sizeof(CoinWorkDouble),numberRows,fp);
    5757    fclose(fp);
    5858  } else {
     
    6161}
    6262#endif
    63 static double eScale=1.0e27;
    64 static double eBaseCaution=1.0e-12;
    65 static double eBase=1.0e-12;
    66 static double eRatio=1.0e40;
    67 static double eRatioCaution=1.0e25;
    68 static double eDiagonal=1.0e25;
    69 static double eDiagonalCaution=1.0e18;
    70 static double eExtra=1.0e-12;
     63// Could change on CLP_LONG_CHOLESKY or COIN_LONG_WORK?
     64static CoinWorkDouble eScale=1.0e27;
     65static CoinWorkDouble eBaseCaution=1.0e-12;
     66static CoinWorkDouble eBase=1.0e-12;
     67static CoinWorkDouble eRatio=1.0e40;
     68static CoinWorkDouble eRatioCaution=1.0e25;
     69static CoinWorkDouble eDiagonal=1.0e25;
     70static CoinWorkDouble eDiagonalCaution=1.0e18;
     71static CoinWorkDouble eExtra=1.0e-12;
    7172
    7273// main function
     
    8182    return 2;
    8283  }
     84#if COIN_LONG_WORK
     85  // reallocate some regions
     86  double * dualSave = dual_;
     87  dual_ = reinterpret_cast<double *>(new CoinWorkDouble[numberRows_]);
     88  double * reducedCostSave = reducedCost_;
     89  reducedCost_ = reinterpret_cast<double *>(new CoinWorkDouble[numberColumns_]);
     90#endif 
    8391  //diagonalPerturbation_=1.0e-25;
    8492  ClpMatrixBase * saveMatrix = NULL;
     
    165173    return -1;
    166174  }
     175  CoinWorkDouble * dualArray = reinterpret_cast<CoinWorkDouble *>(dual_);
    167176  // Could try centering steps without any original step i.e. just center
    168177  //firstFactorization(false);
    169   CoinZeroN(dual_,numberRows_);
     178  CoinZeroN(dualArray,numberRows_);
    170179  multiplyAdd(solution_+numberColumns_,numberRows_,-1.0,errorRegion_,0.0);
    171180  matrix_->times(1.0,solution_,errorRegion_);
    172181  maximumRHSError_=maximumAbsElement(errorRegion_,numberRows_);
    173182  maximumBoundInfeasibility_=maximumRHSError_;
    174   //double maximumDualError_=COIN_DBL_MAX;
     183  //CoinWorkDouble maximumDualError_=COIN_DBL_MAX;
    175184  //initialize
    176185  actualDualStep_=0.0;
     
    184193  int numberFixedTotal=numberFixed;
    185194  //int numberRows_DroppedBefore=0;
    186   //double extra=eExtra;
    187   //double maximumPerturbation=COIN_DBL_MAX;
     195  //CoinWorkDouble extra=eExtra;
     196  //CoinWorkDouble maximumPerturbation=COIN_DBL_MAX;
    188197  //constants for infeas interior point
    189   const double beta2 = 0.99995;
    190   const double tau   = 0.00002;
    191   double lastComplementarityGap=COIN_DBL_MAX;
    192   double lastStep=1.0;
     198  const CoinWorkDouble beta2 = 0.99995;
     199  const CoinWorkDouble tau   = 0.00002;
     200  CoinWorkDouble lastComplementarityGap=COIN_DBL_MAX;
     201  CoinWorkDouble lastStep=1.0;
    193202  // use to see if to take affine
    194   double checkGap = COIN_DBL_MAX;
     203  CoinWorkDouble checkGap = COIN_DBL_MAX;
    195204  int lastGoodIteration=0;
    196   double bestObjectiveGap=COIN_DBL_MAX;
    197   double bestObjective=COIN_DBL_MAX;
     205  CoinWorkDouble bestObjectiveGap=COIN_DBL_MAX;
     206  CoinWorkDouble bestObjective=COIN_DBL_MAX;
     207  int bestKilled=-1;
    198208  int saveIteration=-1;
     209  int saveIteration2=-1;
    199210  bool sloppyOptimal=false;
    200   double * savePi=NULL;
    201   double * savePrimal=NULL;
     211  CoinWorkDouble * savePi=NULL;
     212  CoinWorkDouble * savePrimal=NULL;
     213  CoinWorkDouble * savePi2=NULL;
     214  CoinWorkDouble * savePrimal2=NULL;
    202215  // Extra regions for centering
    203   double * saveX = new double[numberTotal];
    204   double * saveY = new double[numberRows_];
    205   double * saveZ = new double[numberTotal];
    206   double * saveW = new double[numberTotal];
    207   double * saveSL = new double[numberTotal];
    208   double * saveSU = new double[numberTotal];
     216  CoinWorkDouble * saveX = new CoinWorkDouble[numberTotal];
     217  CoinWorkDouble * saveY = new CoinWorkDouble[numberRows_];
     218  CoinWorkDouble * saveZ = new CoinWorkDouble[numberTotal];
     219  CoinWorkDouble * saveW = new CoinWorkDouble[numberTotal];
     220  CoinWorkDouble * saveSL = new CoinWorkDouble[numberTotal];
     221  CoinWorkDouble * saveSU = new CoinWorkDouble[numberTotal];
    209222  // Save smallest mu used in primal dual moves
    210   double smallestPrimalDualMu=COIN_DBL_MAX;
    211   double objScale = optimizationDirection_/
     223  CoinWorkDouble smallestPrimalDualMu=COIN_DBL_MAX;
     224  CoinWorkDouble objScale = optimizationDirection_/
    212225    (rhsScale_*objectiveScale_);
    213226  while (problemStatus_<0) {
     
    231244    handler_->message(CLP_BARRIER_ITERATION,messages_)
    232245      <<numberIterations_
    233       <<primalObjective_*objScale- dblParam_[ClpObjOffset]
    234       << dualObjective_*objScale- dblParam_[ClpObjOffset]
    235       <<complementarityGap_
     246      <<static_cast<double>(primalObjective_*objScale- dblParam_[ClpObjOffset])
     247      << static_cast<double>(dualObjective_*objScale- dblParam_[ClpObjOffset])
     248      <<static_cast<double>(complementarityGap_)
    236249      <<numberFixedTotal
    237250      <<cholesky_->rank()
     
    253266    //saveIteration=-1;
    254267    lastStep = CoinMin(actualPrimalStep_,actualDualStep_);
    255     double goodGapChange;
    256     //#define KEEP_GOING_IF_FIXED 0
     268    CoinWorkDouble goodGapChange;
     269    //#define KEEP_GOING_IF_FIXED 5
    257270#ifndef KEEP_GOING_IF_FIXED
    258271#define KEEP_GOING_IF_FIXED 10000
     
    265278        goodGapChange=0.99; // make more likely to carry on
    266279    }
    267     double gapO;
    268     double lastGood=bestObjectiveGap;
     280    CoinWorkDouble gapO;
     281    CoinWorkDouble lastGood=bestObjectiveGap;
    269282    if (gonePrimalFeasible_&&goneDualFeasible_) {
    270       double largestObjective;
    271       if (fabs(primalObjective_)>fabs(dualObjective_)) {
    272         largestObjective = fabs(primalObjective_);
     283      CoinWorkDouble largestObjective;
     284      if (CoinAbs(primalObjective_)>CoinAbs(dualObjective_)) {
     285        largestObjective = CoinAbs(primalObjective_);
    273286      } else {
    274         largestObjective = fabs(dualObjective_);
     287        largestObjective = CoinAbs(dualObjective_);
    275288      }
    276289      if (largestObjective<1.0) {
    277290        largestObjective=1.0;
    278291      }
    279       gapO=fabs(primalObjective_-dualObjective_)/largestObjective;
     292      gapO=CoinAbs(primalObjective_-dualObjective_)/largestObjective;
    280293      handler_->message(CLP_BARRIER_OBJECTIVE_GAP,messages_)
    281         <<gapO
     294        <<static_cast<double>(gapO)
    282295        <<CoinMessageEol;
    283296      //start saving best
    284       if (gapO<bestObjectiveGap)
     297      bool saveIt=false;
     298      if (gapO<bestObjectiveGap) {
    285299        bestObjectiveGap=gapO;
     300#ifndef SAVE_ON_OBJ
     301        saveIt=true;
     302#endif
     303      }
    286304      if (primalObjective_<bestObjective) {
     305        bestObjective=primalObjective_;
     306#ifdef SAVE_ON_OBJ
     307        saveIt=true;
     308#endif
     309      }
     310      if (numberFixedTotal>bestKilled) {
     311        bestKilled=numberFixedTotal;
     312#if KEEP_GOING_IF_FIXED<10
     313        saveIt=true;
     314#endif
     315      }
     316      if (saveIt) {
     317#if KEEP_GOING_IF_FIXED<10
     318        printf("saving\n");
     319#endif
    287320        saveIteration=numberIterations_;
    288         bestObjective=primalObjective_;
    289321        if (!savePi) {
    290           savePi=new double[numberRows_];
    291           savePrimal = new double [numberTotal];
    292         }
    293         CoinMemcpyN(dual_,numberRows_,savePi);
     322          savePi=new CoinWorkDouble[numberRows_];
     323          savePrimal = new CoinWorkDouble [numberTotal];
     324        }
     325        CoinMemcpyN(dualArray,numberRows_,savePi);
    294326        CoinMemcpyN(solution_,numberTotal,savePrimal);
    295327      } else if(gapO>2.0*bestObjectiveGap) {
     
    300332      //std::cout <<"could stop"<<std::endl;
    301333      //gapO=0.0;
    302       if (fabs(primalObjective_-dualObjective_)<dualTolerance()) {
     334      if (CoinAbs(primalObjective_-dualObjective_)<dualTolerance()) {
    303335        gapO=0.0;
    304336      }
     
    308340        handler_->message(CLP_BARRIER_GONE_INFEASIBLE,messages_)
    309341          <<CoinMessageEol;
     342        CoinWorkDouble scaledRHSError=maximumRHSError_/(solutionNorm_+10.0);
     343        // save alternate
     344        if (numberFixedTotal>bestKilled&&
     345            maximumBoundInfeasibility_<1.0e-6&&
     346            scaledRHSError<1.0e-2) {
     347          bestKilled=numberFixedTotal;
     348#if KEEP_GOING_IF_FIXED<10
     349          printf("saving alternate\n");
     350#endif
     351          saveIteration2=numberIterations_;
     352          if (!savePi2) {
     353            savePi2=new CoinWorkDouble[numberRows_];
     354            savePrimal2 = new CoinWorkDouble [numberTotal];
     355          }
     356          CoinMemcpyN(dualArray,numberRows_,savePi2);
     357          CoinMemcpyN(solution_,numberTotal,savePrimal2);
     358        }
    310359        if (sloppyOptimal) {
    311360          // vaguely optimal
    312           double scaledRHSError=maximumRHSError_/solutionNorm_;
    313361          if (maximumBoundInfeasibility_>1.0e-2||
    314362              scaledRHSError>1.0e-2||
     
    322370        } else {
    323371          // not close to optimal but check if getting bad
    324           double scaledRHSError=maximumRHSError_/solutionNorm_;
     372          CoinWorkDouble scaledRHSError=maximumRHSError_/(solutionNorm_+10.0);
    325373          if ((maximumBoundInfeasibility_>1.0e-1||
    326374              scaledRHSError>1.0e-1||
     
    354402      sloppyOptimal=true;
    355403      handler_->message(CLP_BARRIER_CLOSE_TO_OPTIMAL,messages_)
    356         <<numberIterations_<<complementarityGap_
     404        <<numberIterations_<<static_cast<double>(complementarityGap_)
    357405        <<CoinMessageEol;
    358406    }
     
    363411    if (complementarityGap_>=1.05*lastComplementarityGap) {
    364412      handler_->message(CLP_BARRIER_COMPLEMENTARITY,messages_)
    365         <<complementarityGap_<<"increasing"
     413        <<static_cast<double>(complementarityGap_)<<"increasing"
    366414        <<CoinMessageEol;
    367415      if (saveIteration>=0&&sloppyOptimal) {
     
    380428               complementarityGap_<1.0e-3) {
    381429      handler_->message(CLP_BARRIER_COMPLEMENTARITY,messages_)
    382         <<complementarityGap_<<"not decreasing"
     430        <<static_cast<double>(complementarityGap_)<<"not decreasing"
    383431        <<CoinMessageEol;
    384432      if (gapO>0.75*lastGood&&numberFixed<KEEP_GOING_IF_FIXED) {
     
    388436               complementarityGap_<1.0e-6) {
    389437      handler_->message(CLP_BARRIER_COMPLEMENTARITY,messages_)
    390         <<complementarityGap_<<"not decreasing"
     438        <<static_cast<double>(complementarityGap_)<<"not decreasing"
    391439        <<CoinMessageEol;
    392440      break;
     
    421469    }
    422470    if (gapO<1.0e-9) {
    423       double value=gapO*complementarityGap_;
     471      CoinWorkDouble value=gapO*complementarityGap_;
    424472      value*=actualPrimalStep_;
    425473      value*=actualDualStep_;
     
    433481      }
    434482    }
    435     double nextGap=COIN_DBL_MAX;
     483    CoinWorkDouble nextGap=COIN_DBL_MAX;
    436484    int nextNumber=0;
    437485    int nextNumberItems=0;
     
    441489    //prepare for cholesky.  Set up scaled diagonal in deltaX
    442490    //  ** for efficiency may be better if scale factor known before
    443     double norm2=0.0;
    444     double maximumValue;
     491    CoinWorkDouble norm2=0.0;
     492    CoinWorkDouble maximumValue;
    445493    getNorms(diagonal_,numberTotal,maximumValue,norm2);
    446     diagonalNorm_ = sqrt(norm2/numberComplementarityPairs_);
     494    diagonalNorm_ = CoinSqrt(norm2/numberComplementarityPairs_);
    447495    diagonalScaleFactor_=1.0;
    448     double maximumAllowable=eScale;
     496    CoinWorkDouble maximumAllowable=eScale;
    449497    //scale so largest is less than allowable ? could do better
    450     double factor=0.5;
     498    CoinWorkDouble factor=0.5;
    451499    while (maximumValue>maximumAllowable) {
    452500      diagonalScaleFactor_*=factor;
     
    455503    if (diagonalScaleFactor_!=1.0) {
    456504      handler_->message(CLP_BARRIER_SCALING,messages_)
    457         <<"diagonal"<<diagonalScaleFactor_
     505        <<"diagonal"<<static_cast<double>(diagonalScaleFactor_)
    458506        <<CoinMessageEol;
    459507      diagonalNorm_*=diagonalScaleFactor_;
     
    496544    }
    497545    int phase=0; // predictor, corrector , primal dual
    498     double directionAccuracy=0.0;
     546    CoinWorkDouble directionAccuracy=0.0;
    499547    bool doCorrector=true;
    500548    bool goodMove=true;
     
    517565      debugMove(0,1.0e-2,1.0e-2);
    518566    }
    519     double affineGap=nextGap;
     567    CoinWorkDouble affineGap=nextGap;
    520568    int bestPhase=0;
    521     double bestNextGap=nextGap;
     569    CoinWorkDouble bestNextGap=nextGap;
    522570    // ?
    523571    bestNextGap=CoinMax(nextGap,0.8*complementarityGap_);
     
    526574    if (complementarityGap_>1.0e-4*numberComplementarityPairs_) {
    527575      //std::cout <<"predicted duality gap "<<nextGap<<std::endl;
    528       double part1=nextGap/numberComplementarityPairs_;
     576      CoinWorkDouble part1=nextGap/numberComplementarityPairs_;
    529577      part1=nextGap/numberComplementarityItems_;
    530       double part2=nextGap/complementarityGap_;
     578      CoinWorkDouble part2=nextGap/complementarityGap_;
    531579      mu_=part1*part2*part2;
    532580#if 0
    533       double papermu =complementarityGap_/numberComplementarityPairs_;
    534       double affmu = nextGap/nextNumber;
    535       double sigma = pow(affmu/papermu,3);
     581      CoinWorkDouble papermu =complementarityGap_/numberComplementarityPairs_;
     582      CoinWorkDouble affmu = nextGap/nextNumber;
     583      CoinWorkDouble sigma = pow(affmu/papermu,3);
    536584      printf("mu %g, papermu %g, affmu %g, sigma %g sigmamu %g\n",
    537585             mu_,papermu,affmu,sigma,sigma*papermu);
    538586#endif
    539587      //printf("paper mu %g\n",(nextGap*nextGap*nextGap)/(complementarityGap_*complementarityGap_*
    540       //                                            (double) numberComplementarityPairs_));
     588      //                                            (CoinWorkDouble) numberComplementarityPairs_));
    541589    } else {
    542       double phi;
     590      CoinWorkDouble phi;
    543591      if (numberComplementarityPairs_<=5000) {
    544         phi=pow(static_cast<double> (numberComplementarityPairs_),2.0);
     592        phi=pow(static_cast<CoinWorkDouble> (numberComplementarityPairs_),2.0);
    545593      } else {
    546         phi=pow(static_cast<double> (numberComplementarityPairs_),1.5);
     594        phi=pow(static_cast<CoinWorkDouble> (numberComplementarityPairs_),1.5);
    547595        if (phi<500.0*500.0) {
    548596          phi=500.0*500.0;
     
    552600    }
    553601    //save information
    554     double product=affineProduct();
     602    CoinWorkDouble product=affineProduct();
    555603#if 0
    556604    //can we do corrector step?
    557     double xx= complementarityGap_*(beta2-tau) +product;
     605    CoinWorkDouble xx= complementarityGap_*(beta2-tau) +product;
    558606    if (xx>0.0) {
    559       double saveMu = mu_;
    560       double mu2=numberComplementarityPairs_;
     607      CoinWorkDouble saveMu = mu_;
     608      CoinWorkDouble mu2=numberComplementarityPairs_;
    561609      mu2=xx/mu2;
    562610      if (mu2>mu_) {
     
    593641        bestNextGap=COIN_DBL_MAX;
    594642      }
    595       //double floatNumber = 2.0*numberComplementarityPairs_;
     643      //CoinWorkDouble floatNumber = 2.0*numberComplementarityPairs_;
    596644      //floatNumber = 1.0*numberComplementarityItems_;
    597645      //mu_=nextGap/floatNumber;
     
    615663      CoinMemcpyN(deltaSU_,numberTotal,saveSU);
    616664#ifdef HALVE
    617       double savePrimalStep = actualPrimalStep_;
    618       double saveDualStep = actualDualStep_;
    619       double saveMu = mu_;
     665      CoinWorkDouble savePrimalStep = actualPrimalStep_;
     666      CoinWorkDouble saveDualStep = actualDualStep_;
     667      CoinWorkDouble saveMu = mu_;
    620668#endif
    621669      //set up for next step
    622670      setupForSolve(phase);
    623       double directionAccuracy2=findDirectionVector(phase);
     671      CoinWorkDouble directionAccuracy2=findDirectionVector(phase);
    624672      if (directionAccuracy2>worstDirectionAccuracy_) {
    625673        worstDirectionAccuracy_=directionAccuracy2;
    626674      }
    627       double testValue=1.0e2*directionAccuracy;
     675      CoinWorkDouble testValue=1.0e2*directionAccuracy;
    628676      if (1.0e2*projectionTolerance_>testValue) {
    629677        testValue=1.0e2*projectionTolerance_;
     
    644692      if (goodMove) {
    645693        phase=1;
    646         double norm = findStepLength(phase);
     694        CoinWorkDouble norm = findStepLength(phase);
    647695        nextGap = complementarityGap(nextNumber,nextNumberItems,1);
    648696        debugMove(1,actualPrimalStep_,actualDualStep_);
     
    668716        //printf("halve %d\n",nHalve);
    669717        nHalve++;
    670         const double lambda=0.5;
     718        const CoinWorkDouble lambda=0.5;
    671719        for (i=0;i<numberRows_;i++)
    672720          deltaY_[i] = lambda*deltaY_[i]+(1.0-lambda)*saveY[i];
     
    699747    if (!goodMove) {
    700748      // Just primal dual step
    701       double floatNumber;
     749      CoinWorkDouble floatNumber;
    702750      floatNumber = 2.0*numberComplementarityPairs_;
    703751      //floatNumber = numberComplementarityItems_;
    704       double saveMu=mu_; // use one from predictor corrector
     752      CoinWorkDouble saveMu=mu_; // use one from predictor corrector
    705753      mu_=complementarityGap_/floatNumber;
    706754      // If going well try small mu
    707       mu_ *= sqrt((1.0-lastStep)/(1.0+10.0*lastStep));
    708       double mu1=mu_;
    709       double phi;
     755      mu_ *= CoinSqrt((1.0-lastStep)/(1.0+10.0*lastStep));
     756      CoinWorkDouble mu1=mu_;
     757      CoinWorkDouble phi;
    710758      if (numberComplementarityPairs_<=500) {
    711         phi=pow(static_cast<double> (numberComplementarityPairs_),2.0);
     759        phi=pow(static_cast<CoinWorkDouble> (numberComplementarityPairs_),2.0);
    712760      } else {
    713         phi=pow(static_cast<double> (numberComplementarityPairs_),1.5);
     761        phi=pow(static_cast<CoinWorkDouble> (numberComplementarityPairs_),1.5);
    714762        if (phi<500.0*500.0) {
    715763          phi=500.0*500.0;
     
    719767      //printf("pd mu %g, alternate %g, smallest %g\n",
    720768      //     mu_,mu1,smallestPrimalDualMu);
    721       mu_ = sqrt(mu_*mu1);
     769      mu_ = CoinSqrt(mu_*mu1);
    722770      mu_=mu1;
    723771      if ((numberIterations_&1)==0||numberIterations_<10)
     
    737785      setupForSolve(2);
    738786      findDirectionVector(2);
    739       double norm=findStepLength(2);
     787      CoinWorkDouble norm=findStepLength(2);
    740788      // just for debug
    741789      bestNextGap = complementarityGap_*1.0005;
     
    760808      smallestPrimalDualMu=mu_;
    761809    if (!goodMove)
    762       mu_=nextGap / (static_cast<double> (nextNumber)*1.1);
     810      mu_=nextGap / (static_cast<CoinWorkDouble> (nextNumber)*1.1);
    763811    //if (quadraticObj)
    764812    //goodMove=true;
     
    766814    // Do centering steps
    767815    int numberTries=0;
    768     double nextCenterGap=0.0;
     816    CoinWorkDouble nextCenterGap=0.0;
    769817    int numberGoodTries=0;
    770     double originalDualStep=actualDualStep_;
    771     double originalPrimalStep=actualPrimalStep_;
     818    CoinWorkDouble originalDualStep=actualDualStep_;
     819    CoinWorkDouble originalPrimalStep=actualPrimalStep_;
    772820    if (actualDualStep_>0.9&&actualPrimalStep_>0.9)
    773821      goodMove=false; // don't bother
     
    781829      CoinMemcpyN(deltaZ_,numberTotal,saveZ);
    782830      CoinMemcpyN(deltaW_,numberTotal,saveW);
    783       double savePrimalStep = actualPrimalStep_;
    784       double saveDualStep = actualDualStep_;
    785       double saveMu = mu_;
     831      CoinWorkDouble savePrimalStep = actualPrimalStep_;
     832      CoinWorkDouble saveDualStep = actualDualStep_;
     833      CoinWorkDouble saveMu = mu_;
    786834      setupForSolve(3);
    787835      findDirectionVector(3);
     
    789837      debugMove(3,actualPrimalStep_,actualDualStep_);
    790838      //debugMove(3,1.0e-7,1.0e-7);
    791       double xGap = complementarityGap(nextNumber,nextNumberItems,3);
     839      CoinWorkDouble xGap = complementarityGap(nextNumber,nextNumberItems,3);
    792840      // If one small then that's the one that counts
    793       double checkDual=saveDualStep;
    794       double checkPrimal=savePrimalStep;
     841      CoinWorkDouble checkDual=saveDualStep;
     842      CoinWorkDouble checkPrimal=savePrimalStep;
    795843      if (checkDual>5.0*checkPrimal) {
    796844        checkDual=2.0*checkPrimal;
     
    848896  delete [] saveSU;
    849897  if (savePi) {
    850     std::cout<<"Restoring from iteration "<<saveIteration<<std::endl;
    851     CoinMemcpyN(savePi,numberRows_,dual_);
    852     CoinMemcpyN(savePrimal,numberTotal,solution_);
     898    if (numberIterations_-saveIteration>20&&
     899        numberIterations_-saveIteration2<5) {
     900#if KEEP_GOING_IF_FIXED<10
     901      std::cout<<"Restoring2 from iteration "<<saveIteration2<<std::endl;
     902#endif
     903      CoinMemcpyN(savePi2,numberRows_,dualArray);
     904      CoinMemcpyN(savePrimal2,numberTotal,solution_);
     905    } else {
     906#if KEEP_GOING_IF_FIXED<10
     907      std::cout<<"Restoring from iteration "<<saveIteration<<std::endl;
     908#endif
     909      CoinMemcpyN(savePi,numberRows_,dualArray);
     910      CoinMemcpyN(savePrimal,numberTotal,solution_);
     911    }
    853912    delete [] savePi;
    854913    delete [] savePrimal;
    855914  }
     915  delete [] savePi2;
     916  delete [] savePrimal2;
    856917  //recompute slacks
    857918  // Split out solution
     
    861922  //unscale objective
    862923  multiplyAdd(NULL,numberTotal,0.0,cost_,scaleFactor_);
    863   multiplyAdd(NULL,numberRows_,0,dual_,scaleFactor_);
     924  multiplyAdd(NULL,numberRows_,0,dualArray,scaleFactor_);
    864925  checkSolution();
    865   CoinMemcpyN(reducedCost_,numberColumns_,dj_);
     926  //CoinMemcpyN(reducedCost_,numberColumns_,dj_);
    866927  // If quadratic use last solution
    867928  // Restore quadratic objective if necessary
     
    872933  }
    873934  handler_->message(CLP_BARRIER_END,messages_)
    874     <<sumPrimalInfeasibilities_
    875     <<sumDualInfeasibilities_
    876     <<complementarityGap_
    877     <<objectiveValue()
     935    <<static_cast<double>(sumPrimalInfeasibilities_)
     936    <<static_cast<double>(sumDualInfeasibilities_)
     937    <<static_cast<double>(complementarityGap_)
     938    <<static_cast<double>(objectiveValue())
    878939    <<CoinMessageEol;
    879940  //#ifdef SOME_DEBUG
     
    886947  //std::cout<<"Primal objective "<<objectiveValue()<<std::endl;
    887948  //std::cout<<"maximum complementarity "<<worstComplementarity_<<std::endl;
     949#if COIN_LONG_WORK
     950  // put back dual
     951  delete [] dual_;
     952  delete [] reducedCost_;
     953  dual_ = dualSave;
     954  reducedCost_=reducedCostSave;
     955#endif
    888956  //delete all temporary regions
    889957  deleteWorkingData();
     958#if KEEP_GOING_IF_FIXED<10
     959#ifndef NDEBUG
     960  {
     961    static int kk=0;
     962    char name[20];
     963    sprintf(name,"save.sol.%d",kk);
     964    kk++;
     965    printf("saving to file %s\n",name);
     966    FILE * fp = fopen(name,"wb");
     967    int numberWritten;
     968    numberWritten = fwrite(&numberColumns_,sizeof(int),1,fp);
     969    assert (numberWritten==1);
     970    numberWritten = fwrite(columnActivity_,sizeof(double),numberColumns_,fp);
     971    assert (numberWritten==numberColumns_);
     972    fclose(fp);
     973  }
     974#endif
     975#endif
    890976  if (saveMatrix) {
    891977    // restore normal copy
     
    899985//         1 corrector
    900986//         2 primal dual
    901 double ClpPredictorCorrector::findStepLength( int phase)
     987CoinWorkDouble ClpPredictorCorrector::findStepLength( int phase)
    902988{
    903   double directionNorm=0.0;
    904   double maximumPrimalStep=COIN_DBL_MAX;
    905   double maximumDualStep=COIN_DBL_MAX;
     989  CoinWorkDouble directionNorm=0.0;
     990  CoinWorkDouble maximumPrimalStep=COIN_DBL_MAX;
     991  CoinWorkDouble maximumDualStep=COIN_DBL_MAX;
    906992  int numberTotal = numberRows_+numberColumns_;
    907   double tolerance = 1.0e-12;
     993#if CLP_LONG_CHOLESKY<2
     994  CoinWorkDouble tolerance = 1.0e-12;
     995#else
     996  CoinWorkDouble tolerance = 1.0e-14;
     997#endif
    908998  int chosenPrimalSequence=-1;
    909999  int chosenDualSequence=-1;
     
    9111001  bool lowDual=false;
    9121002  // If done many iterations then allow to hit boundary
    913   double hitTolerance;
     1003  CoinWorkDouble hitTolerance;
    9141004  //printf("objective norm %g\n",objectiveNorm_);
    9151005  if (numberIterations_<80||!gonePrimalFeasible_)
     
    9221012  for (iColumn=0;iColumn<numberTotal;iColumn++) {
    9231013    if (!flagged(iColumn)) {
    924       double directionElement=deltaX_[iColumn];
    925       if (directionNorm<fabs(directionElement)) {
    926         directionNorm=fabs(directionElement);
    927       }
    928       if (0)
    929       printf("%d %g %g %g %g %g %g %g %g %g %g %g\n",
    930              iColumn,solution_[iColumn],deltaX_[iColumn],
    931              lowerSlack_[iColumn],deltaSL_[iColumn],
    932              upperSlack_[iColumn],deltaSU_[iColumn],
    933              dj_[iColumn],
    934              zVec_[iColumn],deltaZ_[iColumn],
    935              wVec_[iColumn],deltaW_[iColumn]);
     1014      CoinWorkDouble directionElement=deltaX_[iColumn];
     1015      if (directionNorm<CoinAbs(directionElement)) {
     1016        directionNorm=CoinAbs(directionElement);
     1017      }
    9361018      if (lowerBound(iColumn)) {
    937         double delta = - deltaSL_[iColumn];
    938         double z1 = deltaZ_[iColumn];
    939         double newZ = zVec_[iColumn]+z1;
     1019        CoinWorkDouble delta = - deltaSL_[iColumn];
     1020        CoinWorkDouble z1 = deltaZ_[iColumn];
     1021        CoinWorkDouble newZ = zVec_[iColumn]+z1;
    9401022        if (zVec_[iColumn]>tolerance) {
    9411023          if (zVec_[iColumn]<-z1*maximumDualStep) {
     
    9461028        }
    9471029        if (lowerSlack_[iColumn]<maximumPrimalStep*delta) {
    948           double newStep=lowerSlack_[iColumn]/delta;
     1030          CoinWorkDouble newStep=lowerSlack_[iColumn]/delta;
    9491031          if (newStep>0.2||newZ<hitTolerance||delta>1.0e3||delta<=1.0e-6||dj_[iColumn]<hitTolerance) {
    9501032            maximumPrimalStep = newStep;
     
    9571039      }
    9581040      if (upperBound(iColumn)) {
    959         double delta = - deltaSU_[iColumn];;
    960         double w1 = deltaW_[iColumn];
    961         double newT = wVec_[iColumn]+w1;
     1041        CoinWorkDouble delta = - deltaSU_[iColumn];;
     1042        CoinWorkDouble w1 = deltaW_[iColumn];
     1043        CoinWorkDouble newT = wVec_[iColumn]+w1;
    9621044        if (wVec_[iColumn]>tolerance) {
    9631045          if (wVec_[iColumn]<-w1*maximumDualStep) {
     
    9681050        }
    9691051        if (upperSlack_[iColumn]<maximumPrimalStep*delta) {
    970           double newStep=upperSlack_[iColumn]/delta;
     1052          CoinWorkDouble newStep=upperSlack_[iColumn]/delta;
    9711053          if (newStep>0.2||newT<hitTolerance||delta>1.0e3||delta<=1.0e-6||dj_[iColumn]>-hitTolerance) {
    9721054            maximumPrimalStep = newStep;
     
    10261108  if (quadraticObj) {
    10271109    // Use smaller unless very small
    1028     double smallerStep=CoinMin(actualDualStep_,actualPrimalStep_);
     1110    CoinWorkDouble smallerStep=CoinMin(actualDualStep_,actualPrimalStep_);
    10291111    if (smallerStep>0.0001) {
    10301112      actualDualStep_=smallerStep;
     
    10421124      // at first - just check better - if not
    10431125      // Complementarity gap will be a*change*change + b*change +c
    1044       double a=0.0;
    1045       double b=0.0;
    1046       double c=0.0;
     1126      CoinWorkDouble a=0.0;
     1127      CoinWorkDouble b=0.0;
     1128      CoinWorkDouble c=0.0;
    10471129      /* SQUARE of dual infeasibility will be:
    10481130         square of dj - ......
    10491131      */
    1050       double aq=0.0;
    1051       double bq=0.0;
    1052       double cq=0.0;
    1053       double gamma2 = gamma_*gamma_; // gamma*gamma will be added to diagonal
    1054       double * linearDjChange = new double[numberTotal];
     1132      CoinWorkDouble aq=0.0;
     1133      CoinWorkDouble bq=0.0;
     1134      CoinWorkDouble cq=0.0;
     1135      CoinWorkDouble gamma2 = gamma_*gamma_; // gamma*gamma will be added to diagonal
     1136      CoinWorkDouble * linearDjChange = new CoinWorkDouble[numberTotal];
    10551137      CoinZeroN(linearDjChange,numberColumns_);
    10561138      multiplyAdd(deltaY_,numberRows_,1.0,linearDjChange+numberColumns_,0.0);
     
    10601142      const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
    10611143      const int * columnQuadraticLength = quadratic->getVectorLengths();
    1062       double * quadraticElement = quadratic->getMutableElements();
     1144      CoinWorkDouble * quadraticElement = quadratic->getMutableElements();
    10631145      for (iColumn=0;iColumn<numberTotal;iColumn++) {
    1064         double oldPrimal = solution_[iColumn];
     1146        CoinWorkDouble oldPrimal = solution_[iColumn];
    10651147        if (!flagged(iColumn)) {
    10661148          if (lowerBound(iColumn)) {
    1067             double change =oldPrimal+deltaX_[iColumn]-lowerSlack_[iColumn]-lower_[iColumn];
     1149            CoinWorkDouble change =oldPrimal+deltaX_[iColumn]-lowerSlack_[iColumn]-lower_[iColumn];
    10681150            c += lowerSlack_[iColumn]*zVec_[iColumn];
    10691151            b += lowerSlack_[iColumn]*deltaZ_[iColumn]+zVec_[iColumn]*change;
     
    10711153          }
    10721154          if (upperBound(iColumn)) {
    1073             double change =upper_[iColumn]-oldPrimal-deltaX_[iColumn]-upperSlack_[iColumn];
     1155            CoinWorkDouble change =upper_[iColumn]-oldPrimal-deltaX_[iColumn]-upperSlack_[iColumn];
    10741156            c += upperSlack_[iColumn]*wVec_[iColumn];
    10751157            b += upperSlack_[iColumn]*deltaW_[iColumn]+wVec_[iColumn]*change;
     
    10771159          }
    10781160          // new djs are dj_ + change*value
    1079           double djChange = linearDjChange[iColumn];
     1161          CoinWorkDouble djChange = linearDjChange[iColumn];
    10801162          if (iColumn<numberColumns_) {
    10811163            for (CoinBigIndex j=columnQuadraticStart[iColumn];
    10821164                 j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
    10831165              int jColumn = columnQuadratic[j];
    1084               double changeJ = deltaX_[jColumn];
    1085               double elementValue = quadraticElement[j];
     1166              CoinWorkDouble changeJ = deltaX_[jColumn];
     1167              CoinWorkDouble elementValue = quadraticElement[j];
    10861168              djChange += changeJ*elementValue;
    10871169            }
    10881170          }
    1089           double gammaTerm = gamma2;
     1171          CoinWorkDouble gammaTerm = gamma2;
    10901172          if (primalR_) {
    10911173            gammaTerm += primalR_[iColumn];
     
    10931175          djChange += gammaTerm;
    10941176          // and dual infeasibility
    1095           double oldInf = dj_[iColumn]-zVec_[iColumn]+wVec_[iColumn]+
     1177          CoinWorkDouble oldInf = dj_[iColumn]-zVec_[iColumn]+wVec_[iColumn]+
    10961178            gammaTerm*solution_[iColumn];
    1097           double changeInf = djChange-deltaZ_[iColumn]+deltaW_[iColumn];
     1179          CoinWorkDouble changeInf = djChange-deltaZ_[iColumn]+deltaW_[iColumn];
    10981180          cq += oldInf*oldInf;
    10991181          bq += 2.0*oldInf*changeInf;
     
    11081190          }
    11091191          // new djs are dj_ + change*value
    1110           double djChange = linearDjChange[iColumn];
     1192          CoinWorkDouble djChange = linearDjChange[iColumn];
    11111193          if (iColumn<numberColumns_) {
    11121194            for (CoinBigIndex j=columnQuadraticStart[iColumn];
    11131195                 j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
    11141196              int jColumn = columnQuadratic[j];
    1115               double changeJ = deltaX_[jColumn];
    1116               double elementValue = quadraticElement[j];
     1197              CoinWorkDouble changeJ = deltaX_[jColumn];
     1198              CoinWorkDouble elementValue = quadraticElement[j];
    11171199              djChange += changeJ*elementValue;
    11181200            }
    11191201          }
    1120           double gammaTerm = gamma2;
     1202          CoinWorkDouble gammaTerm = gamma2;
    11211203          if (primalR_) {
    11221204            gammaTerm += primalR_[iColumn];
     
    11241206          djChange += gammaTerm;
    11251207          // and dual infeasibility
    1126           double oldInf = dj_[iColumn]-zVec_[iColumn]+wVec_[iColumn]+
     1208          CoinWorkDouble oldInf = dj_[iColumn]-zVec_[iColumn]+wVec_[iColumn]+
    11271209            gammaTerm*solution_[iColumn];
    1128           double changeInf = djChange-deltaZ_[iColumn]+deltaW_[iColumn];
     1210          CoinWorkDouble changeInf = djChange-deltaZ_[iColumn]+deltaW_[iColumn];
    11291211          cq += oldInf*oldInf;
    11301212          bq += 2.0*oldInf*changeInf;
     
    11361218      // maybe use inf and do line search
    11371219      // To check see if matches at current step
    1138       double step=actualPrimalStep_;
    1139       //Current gap + solutionNorm_ * sqrt (sum square inf)
    1140       double multiplier = solutionNorm_;
     1220      CoinWorkDouble step=actualPrimalStep_;
     1221      //Current gap + solutionNorm_ * CoinSqrt (sum square inf)
     1222      CoinWorkDouble multiplier = solutionNorm_;
    11411223      multiplier *= 0.01;
    11421224      multiplier=1.0;
    1143       double currentInf =  multiplier*sqrt(cq);
    1144       double nextInf =  multiplier*sqrt(CoinMax(cq+step*bq+step*step*aq,0.0));
    1145       double allowedIncrease=1.4;
     1225      CoinWorkDouble currentInf =  multiplier*CoinSqrt(cq);
     1226      CoinWorkDouble nextInf =  multiplier*CoinSqrt(CoinMax(cq+step*bq+step*step*aq,0.0));
     1227      CoinWorkDouble allowedIncrease=1.4;
    11461228#ifdef SOME_DEBUG
    11471229      printf("lin %g %g %g -> %g\n",a,b,c,
     
    11611243        //cq = CoinMax(cq,10.0);
    11621244        // convert to (x+q)*(x+q) = w
    1163         double q = bq/(1.0*aq);
    1164         double w = CoinMax(q*q + (cq/aq)*(allowedIncrease-1.0),0.0);
    1165         w = sqrt(w);
    1166         double stepX = w-q;
     1245        CoinWorkDouble q = bq/(1.0*aq);
     1246        CoinWorkDouble w = CoinMax(q*q + (cq/aq)*(allowedIncrease-1.0),0.0);
     1247        w = CoinSqrt(w);
     1248        CoinWorkDouble stepX = w-q;
    11671249        step=stepX;
    11681250        nextInf =
    1169           multiplier*sqrt(CoinMax(cq+step*bq+step*step*aq,0.0));
     1251          multiplier*CoinSqrt(CoinMax(cq+step*bq+step*step*aq,0.0));
    11701252#ifdef SOME_DEBUG
    11711253        printf ("with step of %g dualInf is %g\n",
     
    11801262    // minimize complementarity
    11811263    // Complementarity gap will be a*change*change + b*change +c
    1182     double a=0.0;
    1183     double b=0.0;
    1184     double c=0.0;
     1264    CoinWorkDouble a=0.0;
     1265    CoinWorkDouble b=0.0;
     1266    CoinWorkDouble c=0.0;
    11851267    for (iColumn=0;iColumn<numberTotal;iColumn++) {
    1186       double oldPrimal = solution_[iColumn];
     1268      CoinWorkDouble oldPrimal = solution_[iColumn];
    11871269      if (!flagged(iColumn)) {
    11881270        if (lowerBound(iColumn)) {
    1189           double change =oldPrimal+deltaX_[iColumn]-lowerSlack_[iColumn]-lower_[iColumn];
     1271          CoinWorkDouble change =oldPrimal+deltaX_[iColumn]-lowerSlack_[iColumn]-lower_[iColumn];
    11901272          c += lowerSlack_[iColumn]*zVec_[iColumn];
    11911273          b += lowerSlack_[iColumn]*deltaZ_[iColumn]+zVec_[iColumn]*change;
     
    11931275        }
    11941276        if (upperBound(iColumn)) {
    1195           double change =upper_[iColumn]-oldPrimal-deltaX_[iColumn]-upperSlack_[iColumn];
     1277          CoinWorkDouble change =upper_[iColumn]-oldPrimal-deltaX_[iColumn]-upperSlack_[iColumn];
    11961278          c += upperSlack_[iColumn]*wVec_[iColumn];
    11971279          b += upperSlack_[iColumn]*deltaW_[iColumn]+wVec_[iColumn]*change;
     
    12111293    // maybe use inf and do line search
    12121294    // To check see if matches at current step
    1213     double step=CoinMin(actualPrimalStep_,actualDualStep_);
    1214     double next = c+b*step+a*step*step;
     1295    CoinWorkDouble step=CoinMin(actualPrimalStep_,actualDualStep_);
     1296    CoinWorkDouble next = c+b*step+a*step*step;
    12151297#ifdef SOME_DEBUG
    12161298    printf("lin %g %g %g -> %g\n",a,b,c,
     
    12301312    a=0.0;
    12311313    b=0.0;
    1232     double ratio = actualDualStep_/actualPrimalStep_;
     1314    CoinWorkDouble ratio = actualDualStep_/actualPrimalStep_;
    12331315    for (iColumn=0;iColumn<numberTotal;iColumn++) {
    1234       double oldPrimal = solution_[iColumn];
     1316      CoinWorkDouble oldPrimal = solution_[iColumn];
    12351317      if (!flagged(iColumn)) {
    12361318        if (lowerBound(iColumn)) {
    1237           double change =oldPrimal+deltaX_[iColumn]-lowerSlack_[iColumn]-lower_[iColumn];
     1319          CoinWorkDouble change =oldPrimal+deltaX_[iColumn]-lowerSlack_[iColumn]-lower_[iColumn];
    12381320          b += lowerSlack_[iColumn]*deltaZ_[iColumn]*ratio+zVec_[iColumn]*change;
    12391321          a += deltaZ_[iColumn]*change*ratio;
    12401322        }
    12411323        if (upperBound(iColumn)) {
    1242           double change =upper_[iColumn]-oldPrimal-deltaX_[iColumn]-upperSlack_[iColumn];
     1324          CoinWorkDouble change =upper_[iColumn]-oldPrimal-deltaX_[iColumn]-upperSlack_[iColumn];
    12431325          b += upperSlack_[iColumn]*deltaW_[iColumn]*ratio+wVec_[iColumn]*change;
    12441326          a += deltaW_[iColumn]*change*ratio;
     
    12501332    // To check see if matches at current step
    12511333    step=actualPrimalStep_;
    1252     double next2 = c+b*step+a*step*step;
     1334    CoinWorkDouble next2 = c+b*step+a*step*step;
    12531335    if (next2>next) {
    12541336      actualPrimalStep_=CoinMin(actualPrimalStep_,actualDualStep_);
     
    12751357#ifdef FULL_DEBUG
    12761358  if (phase==3){
    1277     double minBeta = 0.1*mu_;
    1278     double maxBeta = 10.0*mu_;
     1359    CoinWorkDouble minBeta = 0.1*mu_;
     1360    CoinWorkDouble maxBeta = 10.0*mu_;
    12791361    for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
    12801362      if (!flagged(iColumn)) {
    12811363        if (lowerBound(iColumn)) {
    1282           double change = -rhsL_[iColumn] + deltaX_[iColumn];
    1283           double dualValue=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];
    1284           double primalValue=lowerSlack_[iColumn]+actualPrimalStep_*change;
    1285           double gapProduct=dualValue*primalValue;
     1364          CoinWorkDouble change = -rhsL_[iColumn] + deltaX_[iColumn];
     1365          CoinWorkDouble dualValue=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];
     1366          CoinWorkDouble primalValue=lowerSlack_[iColumn]+actualPrimalStep_*change;
     1367          CoinWorkDouble gapProduct=dualValue*primalValue;
    12861368          if (delta2Z_[iColumn]<minBeta||delta2Z_[iColumn]>maxBeta)
    12871369            printf("3lower %d primal %g, dual %g, gap %g, old gap %g\n",
     
    12891371        } 
    12901372        if (upperBound(iColumn)) {
    1291           double change = rhsU_[iColumn]-deltaX_[iColumn];
    1292           double dualValue=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn];
    1293           double primalValue=upperSlack_[iColumn]+actualPrimalStep_*change;
    1294           double gapProduct=dualValue*primalValue;
     1373          CoinWorkDouble change = rhsU_[iColumn]-deltaX_[iColumn];
     1374          CoinWorkDouble dualValue=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn];
     1375          CoinWorkDouble primalValue=upperSlack_[iColumn]+actualPrimalStep_*change;
     1376          CoinWorkDouble gapProduct=dualValue*primalValue;
    12951377          if (delta2W_[iColumn]<minBeta||delta2W_[iColumn]>maxBeta)
    12961378            printf("3upper %d primal %g, dual %g, gap %g, old gap %g\n",
     
    13031385#ifdef SOME_DEBUG_not
    13041386  {
    1305     double largestL=0.0;
    1306     double smallestL=COIN_DBL_MAX;
    1307     double largestU=0.0;
    1308     double smallestU=COIN_DBL_MAX;
    1309     double sumL=0.0;
    1310     double sumU=0.0;
     1387    CoinWorkDouble largestL=0.0;
     1388    CoinWorkDouble smallestL=COIN_DBL_MAX;
     1389    CoinWorkDouble largestU=0.0;
     1390    CoinWorkDouble smallestU=COIN_DBL_MAX;
     1391    CoinWorkDouble sumL=0.0;
     1392    CoinWorkDouble sumU=0.0;
    13111393    int nL=0;
    13121394    int nU=0;
     
    13141396      if (!flagged(iColumn)) {
    13151397        if (lowerBound(iColumn)) {
    1316           double change = -rhsL_[iColumn] + deltaX_[iColumn];
    1317           double dualValue=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];
    1318           double primalValue=lowerSlack_[iColumn]+actualPrimalStep_*change;
    1319           double gapProduct=dualValue*primalValue;
     1398          CoinWorkDouble change = -rhsL_[iColumn] + deltaX_[iColumn];
     1399          CoinWorkDouble dualValue=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];
     1400          CoinWorkDouble primalValue=lowerSlack_[iColumn]+actualPrimalStep_*change;
     1401          CoinWorkDouble gapProduct=dualValue*primalValue;
    13201402          largestL = CoinMax(largestL,gapProduct);
    13211403          smallestL = CoinMin(smallestL,gapProduct);
     
    13241406        } 
    13251407        if (upperBound(iColumn)) {
    1326           double change = rhsU_[iColumn]-deltaX_[iColumn];
    1327           double dualValue=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn];
    1328           double primalValue=upperSlack_[iColumn]+actualPrimalStep_*change;
    1329           double gapProduct=dualValue*primalValue;
     1408          CoinWorkDouble change = rhsU_[iColumn]-deltaX_[iColumn];
     1409          CoinWorkDouble dualValue=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn];
     1410          CoinWorkDouble primalValue=upperSlack_[iColumn]+actualPrimalStep_*change;
     1411          CoinWorkDouble gapProduct=dualValue*primalValue;
    13301412          largestU = CoinMax(largestU,gapProduct);
    13311413          smallestU = CoinMin(smallestU,gapProduct);
     
    13351417      }
    13361418    }
    1337     double mu = (sumL+sumU)/(static_cast<double> (nL+nU));
     1419    CoinWorkDouble mu = (sumL+sumU)/(static_cast<CoinWorkDouble> (nL+nU));
    13381420
    1339     double minBeta = 0.1*mu;
    1340     double maxBeta = 10.0*mu;
     1421    CoinWorkDouble minBeta = 0.1*mu;
     1422    CoinWorkDouble maxBeta = 10.0*mu;
    13411423    int nBL=0;
    13421424    int nAL=0;
     
    13461428      if (!flagged(iColumn)) {
    13471429        if (lowerBound(iColumn)) {
    1348           double change = -rhsL_[iColumn] + deltaX_[iColumn];
    1349           double dualValue=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];
    1350           double primalValue=lowerSlack_[iColumn]+actualPrimalStep_*change;
    1351           double gapProduct=dualValue*primalValue;
     1430          CoinWorkDouble change = -rhsL_[iColumn] + deltaX_[iColumn];
     1431          CoinWorkDouble dualValue=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];
     1432          CoinWorkDouble primalValue=lowerSlack_[iColumn]+actualPrimalStep_*change;
     1433          CoinWorkDouble gapProduct=dualValue*primalValue;
    13521434          if (gapProduct<minBeta)
    13531435            nBL++;
     
    13591441        } 
    13601442        if (upperBound(iColumn)) {
    1361           double change = rhsU_[iColumn]-deltaX_[iColumn];
    1362           double dualValue=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn];
    1363           double primalValue=upperSlack_[iColumn]+actualPrimalStep_*change;
    1364           double gapProduct=dualValue*primalValue;
     1443          CoinWorkDouble change = rhsU_[iColumn]-deltaX_[iColumn];
     1444          CoinWorkDouble dualValue=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn];
     1445          CoinWorkDouble primalValue=upperSlack_[iColumn]+actualPrimalStep_*change;
     1446          CoinWorkDouble gapProduct=dualValue*primalValue;
    13651447          if (gapProduct<minBeta)
    13661448            nBU++;
     
    13841466/* Does solve. region1 is for deltaX (columns+rows), region2 for deltaPi (rows) */
    13851467void
    1386 ClpPredictorCorrector::solveSystem(double * region1, double * region2,
    1387                                    const double * region1In, const double * region2In,
    1388                                    const double * saveRegion1, const double * saveRegion2,
     1468ClpPredictorCorrector::solveSystem(CoinWorkDouble * region1, CoinWorkDouble * region2,
     1469                                   const CoinWorkDouble * region1In, const CoinWorkDouble * region2In,
     1470                                   const CoinWorkDouble * saveRegion1, const CoinWorkDouble * saveRegion2,
    13891471                                   bool gentleRefine)
    13901472{
     
    14061488    multiplyAdd(region1+numberColumns_,numberRows_,-1.0,region2,1.0);
    14071489    matrix_->times(1.0,region1,region2);
    1408     double maximumRHS = maximumAbsElement(region2,numberRows_);
    1409     double scale=1.0;
    1410     double unscale=1.0;
     1490    CoinWorkDouble maximumRHS = maximumAbsElement(region2,numberRows_);
     1491    CoinWorkDouble scale=1.0;
     1492    CoinWorkDouble unscale=1.0;
    14111493    if (maximumRHS>1.0e-30) {
    14121494      if (maximumRHS<=0.5) {
    1413         double factor=2.0;
     1495        CoinWorkDouble factor=2.0;
    14141496        while (maximumRHS<=0.5) {
    14151497          maximumRHS*=factor;
     
    14171499        } /* endwhile */
    14181500      } else if (maximumRHS>=2.0&&maximumRHS<=COIN_DBL_MAX) {
    1419         double factor=0.5;
     1501        CoinWorkDouble factor=0.5;
    14201502        while (maximumRHS>=2.0) {
    14211503          maximumRHS*=factor;
     
    14441526  if (saveRegion2) {
    14451527    //refine?
    1446     double scaleX=1.0;
     1528    CoinWorkDouble scaleX=1.0;
    14471529    if (gentleRefine)
    14481530      scaleX=0.8;
     
    14531535}
    14541536// findDirectionVector.
    1455 double ClpPredictorCorrector::findDirectionVector(const int phase)
     1537CoinWorkDouble ClpPredictorCorrector::findDirectionVector(const int phase)
    14561538{
    1457   double projectionTolerance=projectionTolerance_;
     1539  CoinWorkDouble projectionTolerance=projectionTolerance_;
    14581540  //temporary
    14591541  //projectionTolerance=1.0e-15;
    1460   double errorCheck=0.9*maximumRHSError_/solutionNorm_;
     1542  CoinWorkDouble errorCheck=0.9*maximumRHSError_/solutionNorm_;
    14611543  if (errorCheck>primalTolerance()) {
    14621544    if (errorCheck<projectionTolerance) {
     
    14681550    }
    14691551  }
    1470   double * newError = new double [numberRows_];
     1552  CoinWorkDouble * newError = new CoinWorkDouble [numberRows_];
    14711553  int numberTotal = numberRows_+numberColumns_;
    14721554  //if flagged then entries zero so can do
    14731555  // For KKT separate out
    1474   double * region1Save=NULL;//for refinement
     1556  CoinWorkDouble * region1Save=NULL;//for refinement
    14751557  int iColumn;
    14761558  if (cholesky_->type()<20) {
     
    14931575  }
    14941576  bool goodSolve=false;
    1495   double * regionSave=NULL;//for refinement
     1577  CoinWorkDouble * regionSave=NULL;//for refinement
    14961578  int numberTries=0;
    1497   double relativeError=COIN_DBL_MAX;
    1498   double tryError=1.0e31;
     1579  CoinWorkDouble relativeError=COIN_DBL_MAX;
     1580  CoinWorkDouble tryError=1.0e31;
     1581  CoinWorkDouble saveMaximum=0.0;
     1582  double firstError=0.0;
     1583  double lastError=0.0;
    14991584  while (!goodSolve&&numberTries<30) {
    1500     double lastError=relativeError;
     1585    CoinWorkDouble lastError=relativeError;
    15011586    goodSolve=true;
    1502     double maximumRHS;
    1503     double saveMaximum;
     1587    CoinWorkDouble maximumRHS;
    15041588    maximumRHS = CoinMax(maximumAbsElement(deltaY_,numberRows_),1.0e-12);
    1505     saveMaximum = maximumRHS;
     1589    if (!numberTries)
     1590      saveMaximum = maximumRHS;
    15061591    if (cholesky_->type()<20) {
    15071592      // no kkt
    1508       double scale=1.0;
    1509       double unscale=1.0;
     1593      CoinWorkDouble scale=1.0;
     1594      CoinWorkDouble unscale=1.0;
    15101595      if (maximumRHS>1.0e-30) {
    15111596        if (maximumRHS<=0.5) {
    1512           double factor=2.0;
     1597          CoinWorkDouble factor=2.0;
    15131598          while (maximumRHS<=0.5) {
    15141599            maximumRHS*=factor;
     
    15161601          } /* endwhile */
    15171602        } else if (maximumRHS>=2.0&&maximumRHS<=COIN_DBL_MAX) {
    1518           double factor=0.5;
     1603          CoinWorkDouble factor=0.5;
    15191604          while (maximumRHS>=2.0) {
    15201605            maximumRHS*=factor;
     
    15441629      if (numberTries) {
    15451630        //refine?
    1546         double scaleX=1.0;
     1631        CoinWorkDouble scaleX=1.0;
    15471632        if (lastError>1.0e-5)
    15481633          scaleX=0.8;
     
    15671652   
    15681653    //now add in old Ax - doing extra checking
    1569     double maximumRHSError=0.0;
    1570     double maximumRHSChange=0.0;
     1654    CoinWorkDouble maximumRHSError=0.0;
     1655    CoinWorkDouble maximumRHSChange=0.0;
    15711656    int iRow;
    15721657    char * dropped = cholesky_->rowsDropped();
    15731658    for (iRow=0;iRow<numberRows_;iRow++) {
    15741659      if (!dropped[iRow]) {
    1575         double newValue=newError[iRow];
    1576         double oldValue=errorRegion_[iRow];
     1660        CoinWorkDouble newValue=newError[iRow];
     1661        CoinWorkDouble oldValue=errorRegion_[iRow];
    15771662        //severity of errors depend on signs
    15781663        //**later                                                             */
    1579         if (fabs(newValue)>maximumRHSChange) {
    1580           maximumRHSChange=fabs(newValue);
     1664        if (CoinAbs(newValue)>maximumRHSChange) {
     1665          maximumRHSChange=CoinAbs(newValue);
    15811666        }
    1582         double result=newValue+oldValue;
    1583         if (fabs(result)>maximumRHSError) {
    1584           maximumRHSError=fabs(result);
     1667        CoinWorkDouble result=newValue+oldValue;
     1668        if (CoinAbs(result)>maximumRHSError) {
     1669          maximumRHSError=CoinAbs(result);
    15851670        }
    15861671        newError[iRow]=result;
    15871672      } else {
    1588         double newValue=newError[iRow];
    1589         double oldValue=errorRegion_[iRow];
    1590         if (fabs(newValue)>maximumRHSChange) {
    1591           maximumRHSChange=fabs(newValue);
     1673        CoinWorkDouble newValue=newError[iRow];
     1674        CoinWorkDouble oldValue=errorRegion_[iRow];
     1675        if (CoinAbs(newValue)>maximumRHSChange) {
     1676          maximumRHSChange=CoinAbs(newValue);
    15921677        }
    1593         double result=newValue+oldValue;
     1678        CoinWorkDouble result=newValue+oldValue;
    15941679        newError[iRow]=result;
    15951680        //newError[iRow]=0.0;
     
    16021687    if (relativeError>tryError)
    16031688      relativeError=tryError;
     1689    if (numberTries==1)
     1690      firstError=relativeError;
    16041691    if (relativeError<lastError) {
     1692      lastError=relativeError;
    16051693      maximumRHSChange_= maximumRHSChange;
    1606       if (relativeError>1.0e-9
    1607           ||numberTries>1) {
    1608         handler_->message(CLP_BARRIER_ACCURACY,messages_)
    1609           <<phase<<numberTries<<relativeError
    1610           <<CoinMessageEol;
    1611       }
    16121694      if (relativeError>projectionTolerance&&numberTries<=3) {
    16131695        //try and refine
     
    16171699      if (!goodSolve) {
    16181700        if (!regionSave) {
    1619           regionSave = new double [numberRows_];
     1701          regionSave = new CoinWorkDouble [numberRows_];
    16201702          if (cholesky_->type()>=20)
    1621             region1Save = new double [numberTotal];
     1703            region1Save = new CoinWorkDouble [numberTotal];
    16221704        }
    16231705        CoinMemcpyN(deltaY_,numberRows_,regionSave);
     
    16531735      } else {
    16541736        // disaster
    1655         CoinFillN(deltaX_,numberTotal,1.0);
    1656         CoinFillN(deltaY_,numberRows_,1.0);
     1737        CoinFillN(deltaX_,numberTotal,static_cast<CoinWorkDouble>(1.0));
     1738        CoinFillN(deltaY_,numberRows_,static_cast<CoinWorkDouble>(1.0));
    16571739        printf("bad cholesky\n");
    16581740      }
    16591741    }
    16601742  } /* endwhile */
     1743  if (firstError>1.0e-9||numberTries>1) {
     1744    handler_->message(CLP_BARRIER_ACCURACY,messages_)
     1745      <<phase<<numberTries<<static_cast<double>(firstError)
     1746      <<static_cast<double>(lastError)
     1747      <<CoinMessageEol;
     1748  }
    16611749  delete [] regionSave;
    16621750  delete [] region1Save;
    16631751  delete [] newError;
    16641752  // now rest
    1665   double extra=eExtra;
     1753  CoinWorkDouble extra=eExtra;
    16661754  //multiplyAdd(deltaY_,numberRows_,1.0,deltaW_+numberColumns_,0.0);
    16671755  //CoinZeroN(deltaW_,numberColumns_);
     
    16721760    deltaSL_[iColumn]=0.0;
    16731761    deltaZ_[iColumn]=0.0;
    1674     double dd=deltaW_[iColumn];
     1762    CoinWorkDouble dd=deltaW_[iColumn];
    16751763    deltaW_[iColumn]=0.0;
    16761764    if (!flagged(iColumn)) {
    1677       double deltaX = deltaX_[iColumn];
     1765      CoinWorkDouble deltaX = deltaX_[iColumn];
    16781766      if (lowerBound(iColumn)) {
    1679         double zValue = rhsZ_[iColumn];
    1680         double gHat = zValue + zVec_[iColumn]*rhsL_[iColumn];
    1681         double slack = lowerSlack_[iColumn]+extra;
     1767        CoinWorkDouble zValue = rhsZ_[iColumn];
     1768        CoinWorkDouble gHat = zValue + zVec_[iColumn]*rhsL_[iColumn];
     1769        CoinWorkDouble slack = lowerSlack_[iColumn]+extra;
    16821770        deltaSL_[iColumn] = -rhsL_[iColumn]+deltaX;
    16831771        deltaZ_[iColumn]=(gHat-zVec_[iColumn]*deltaX)/slack;
    16841772      }
    16851773      if (upperBound(iColumn)) {
    1686         double wValue = rhsW_[iColumn];
    1687         double hHat = wValue - wVec_[iColumn]*rhsU_[iColumn];
    1688         double slack = upperSlack_[iColumn]+extra;
     1774        CoinWorkDouble wValue = rhsW_[iColumn];
     1775        CoinWorkDouble hHat = wValue - wVec_[iColumn]*rhsU_[iColumn];
     1776        CoinWorkDouble slack = upperSlack_[iColumn]+extra;
    16891777        deltaSU_[iColumn] = rhsU_[iColumn]-deltaX;
    16901778        deltaW_[iColumn]=(hHat+wVec_[iColumn]*deltaX)/slack;
     
    16921780      if (0) {
    16931781        // different way of calculating
    1694         double gamma2 = gamma_*gamma_;
    1695         double dZ=0.0;
    1696         double dW=0.0;
    1697         double zValue = rhsZ_[iColumn];
    1698         double gHat = zValue + zVec_[iColumn]*rhsL_[iColumn];
    1699         double slackL = lowerSlack_[iColumn]+extra;
    1700         double wValue = rhsW_[iColumn];
    1701         double hHat = wValue - wVec_[iColumn]*rhsU_[iColumn];
    1702         double slackU = upperSlack_[iColumn]+extra;
    1703         double q = rhsC_[iColumn]+gamma2 * deltaX +dd;
     1782        CoinWorkDouble gamma2 = gamma_*gamma_;
     1783        CoinWorkDouble dZ=0.0;
     1784        CoinWorkDouble dW=0.0;
     1785        CoinWorkDouble zValue = rhsZ_[iColumn];
     1786        CoinWorkDouble gHat = zValue + zVec_[iColumn]*rhsL_[iColumn];
     1787        CoinWorkDouble slackL = lowerSlack_[iColumn]+extra;
     1788        CoinWorkDouble wValue = rhsW_[iColumn];
     1789        CoinWorkDouble hHat = wValue - wVec_[iColumn]*rhsU_[iColumn];
     1790        CoinWorkDouble slackU = upperSlack_[iColumn]+extra;
     1791        CoinWorkDouble q = rhsC_[iColumn]+gamma2 * deltaX +dd;
    17041792        if (primalR_)
    17051793          q += deltaX*primalR_[iColumn];
     
    17281816  }
    17291817#if 0
    1730   double * check = new double[numberTotal]; 
     1818  CoinWorkDouble * check = new CoinWorkDouble[numberTotal]; 
    17311819  // Check out rhsC_
    17321820  multiplyAdd(deltaY_,numberRows_,-1.0,check+numberColumns_,0.0);
     
    17361824  for (iColumn=0;iColumn<numberTotal;iColumn++) {
    17371825    check[iColumn] += deltaZ_[iColumn]-deltaW_[iColumn];
    1738     if (fabs(check[iColumn]-rhsC_[iColumn])>1.0e-3)
     1826    if (CoinAbs(check[iColumn]-rhsC_[iColumn])>1.0e-3)
    17391827      printf("rhsC %d %g %g\n",iColumn,check[iColumn],rhsC_[iColumn]);
    17401828  }
     
    17431831    check[iColumn] += lowerSlack_[iColumn]*deltaZ_[iColumn]+
    17441832      zVec_[iColumn]*deltaSL_[iColumn];
    1745     if (fabs(check[iColumn]-rhsZ_[iColumn])>1.0e-3)
     1833    if (CoinAbs(check[iColumn]-rhsZ_[iColumn])>1.0e-3)
    17461834      printf("rhsZ %d %g %g\n",iColumn,check[iColumn],rhsZ_[iColumn]);
    17471835  }
     
    17501838    check[iColumn] += upperSlack_[iColumn]*deltaW_[iColumn]+
    17511839      wVec_[iColumn]*deltaSU_[iColumn];
    1752     if (fabs(check[iColumn]-rhsW_[iColumn])>1.0e-3)
     1840    if (CoinAbs(check[iColumn]-rhsW_[iColumn])>1.0e-3)
    17531841      printf("rhsW %d %g %g\n",iColumn,check[iColumn],rhsW_[iColumn]);
    17541842  }
     
    17621850  int numberTotal = numberRows_+numberColumns_;
    17631851  int iColumn;
    1764   double tolerance = primalTolerance();
     1852  CoinWorkDouble tolerance = primalTolerance();
    17651853  // See if quadratic objective
    17661854#ifndef NO_RTTI
     
    17831871      clearFixed(iColumn);
    17841872  }
    1785   double initialValue =1.0e-12;
    17861873
    1787   double maximumObjective=0.0;
    1788   double objectiveNorm2=0.0;
     1874  CoinWorkDouble maximumObjective=0.0;
     1875  CoinWorkDouble objectiveNorm2=0.0;
    17891876  getNorms(cost_,numberTotal,maximumObjective,objectiveNorm2);
    17901877  if (!maximumObjective) {
    17911878    maximumObjective=1.0; // objective all zero
    17921879  }
    1793   objectiveNorm2=sqrt(objectiveNorm2)/static_cast<double> (numberTotal);
     1880  objectiveNorm2=CoinSqrt(objectiveNorm2)/static_cast<CoinWorkDouble> (numberTotal);
    17941881  objectiveNorm_=maximumObjective;
    17951882  scaleFactor_=1.0;
     
    18091896  if (quadraticObj) {
    18101897    // If scaled then really scale matrix
    1811     double scaleFactor =
     1898    CoinWorkDouble scaleFactor =
    18121899      scaleFactor_*optimizationDirection_*objectiveScale_*
    18131900      rhsScale_;
     
    18191906      double * quadraticElement = quadratic->getMutableElements();
    18201907      int numberColumns = quadratic->getNumCols();
    1821       double scale = 1.0/scaleFactor;
     1908      CoinWorkDouble scale = 1.0/scaleFactor;
    18221909      if (scalingFlag_>0&&rowScale_) {
    18231910        for (int iColumn=0;iColumn<numberColumns;iColumn++) {
    1824           double scaleI = columnScale_[iColumn]*scale;
     1911          CoinWorkDouble scaleI = columnScale_[iColumn]*scale;
    18251912          for (CoinBigIndex j=columnQuadraticStart[iColumn];
    18261913               j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
    18271914            int jColumn = columnQuadratic[j];
    1828             double scaleJ = columnScale_[jColumn];
     1915            CoinWorkDouble scaleJ = columnScale_[jColumn];
    18291916            quadraticElement[j] *= scaleI*scaleJ;
    1830             objectiveNorm_ = CoinMax(objectiveNorm_,fabs(quadraticElement[j]));
     1917            objectiveNorm_ = CoinMax(objectiveNorm_,CoinAbs(quadraticElement[j]));
    18311918          }
    18321919        }
     
    18371924               j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
    18381925            quadraticElement[j] *= scale;
    1839             objectiveNorm_ = CoinMax(objectiveNorm_,fabs(quadraticElement[j]));
     1926            objectiveNorm_ = CoinMax(objectiveNorm_,CoinAbs(quadraticElement[j]));
    18401927          }
    18411928        }
     
    18481935  //DZ in lowerDual
    18491936  //DW in upperDual
    1850   double infiniteCheck=1.0e40;
    1851   //double     fakeCheck=1.0e10;
     1937  CoinWorkDouble infiniteCheck=1.0e40;
     1938  //CoinWorkDouble     fakeCheck=1.0e10;
    18521939  //use deltaX region for work region
    18531940  for (iColumn=0;iColumn<numberTotal;iColumn++) {
    1854     double primalValue = solution_[iColumn];
     1941    CoinWorkDouble primalValue = solution_[iColumn];
    18551942    clearFlagged(iColumn);
    18561943    clearFixedOrFree(iColumn);
     
    18631950      diagonal_[iColumn]=1.0;
    18641951      deltaX_[iColumn]=1.0;
    1865       double lowerValue=lower_[iColumn];
    1866       double upperValue=upper_[iColumn];
     1952      CoinWorkDouble lowerValue=lower_[iColumn];
     1953      CoinWorkDouble upperValue=upper_[iColumn];
    18671954      if (lowerValue>-infiniteCheck) {
    18681955        if (upperValue<infiniteCheck) {
     
    19532040    solveSystem(deltaX_,errorRegion_,solution_,NULL,NULL,NULL,false);
    19542041  }
    1955   initialValue=1.0e2;
     2042  CoinWorkDouble initialValue = 1.0e2;
    19562043  if (rhsNorm_*1.0e-2>initialValue) {
    19572044    initialValue=rhsNorm_*1.0e-2;
    19582045  }
    19592046  //initialValue = CoinMax(1.0,rhsNorm_);
    1960   double smallestBoundDifference=COIN_DBL_MAX;
    1961   double * fakeSolution = deltaX_;
     2047  CoinWorkDouble smallestBoundDifference=COIN_DBL_MAX;
     2048  CoinWorkDouble * fakeSolution = deltaX_;
    19622049  for ( iColumn=0;iColumn<numberTotal;iColumn++) {
    19632050    if (!flagged(iColumn)) {
     
    19752062  solutionNorm_=1.0e-12;
    19762063  handler_->message(CLP_BARRIER_SAFE,messages_)
    1977     <<initialValue<<objectiveNorm_
     2064    <<static_cast<double>(initialValue)<<static_cast<double>(objectiveNorm_)
    19782065    <<CoinMessageEol;
    1979   double extra=1.0e-10;
    1980   double largeGap=1.0e15;
    1981   //double safeObjectiveValue=2.0*objectiveNorm_;
    1982   double safeObjectiveValue=objectiveNorm_+1.0;
    1983   double safeFree=1.0e-1*initialValue;
     2066  CoinWorkDouble extra=1.0e-10;
     2067  CoinWorkDouble largeGap=1.0e15;
     2068  //CoinWorkDouble safeObjectiveValue=2.0*objectiveNorm_;
     2069  CoinWorkDouble safeObjectiveValue=objectiveNorm_+1.0;
     2070  CoinWorkDouble safeFree=1.0e-1*initialValue;
    19842071  //printf("normal safe dual value of %g, primal value of %g\n",
    19852072  // safeObjectiveValue,initialValue);
     
    19882075  //printf("temp safe dual value of %g, primal value of %g\n",
    19892076  // safeObjectiveValue,initialValue);
    1990   double zwLarge=1.0e2*initialValue;
     2077  CoinWorkDouble zwLarge=1.0e2*initialValue;
    19912078  //zwLarge=1.0e40;
    19922079  if (cholesky_->choleskyCondition()<0.0&&cholesky_->type()<20) {
     
    19962083    safeFree *= 10.0;
    19972084  }
    1998   double gamma2 = gamma_*gamma_; // gamma*gamma will be added to diagonal
    1999 #if 0
    2000   fakeSolution[0 ] =   0.072310129 ;
    2001   fakeSolution[1 ] =   0.053083871;
    2002   fakeSolution[2 ] =      0.178127;
    2003   fakeSolution[3 ] =    0.13215151;
    2004   fakeSolution[4 ] =   0.072715642;
    2005   fakeSolution[5 ] =    0.15680727;
    2006   fakeSolution[6 ] =    0.16841689;
    2007   fakeSolution[7 ] =   0.093612798 ;
    2008   fakeSolution[8 ] =   0.072774891 ;
    2009   fakeSolution[9]=1.0;
    2010   initialValue=1.0e-5;
    2011   safeObjectiveValue=1.0e-5;
    2012 #endif
     2085  CoinWorkDouble gamma2 = gamma_*gamma_; // gamma*gamma will be added to diagonal
    20132086  // First do primal side
    20142087  for ( iColumn=0;iColumn<numberTotal;iColumn++) {
    20152088    if (!flagged(iColumn)) {
    2016       double lowerValue=lower_[iColumn];
    2017       double upperValue=upper_[iColumn];
    2018       double newValue;
    2019       double setPrimal=initialValue;
     2089      CoinWorkDouble lowerValue=lower_[iColumn];
     2090      CoinWorkDouble upperValue=upper_[iColumn];
     2091      CoinWorkDouble newValue;
     2092      CoinWorkDouble setPrimal=initialValue;
    20202093      if (quadraticObj) {
    20212094        // perturb primal solution a bit
     
    20262099          //upper and lower bounds
    20272100          if (upperValue-lowerValue>2.0*setPrimal) {
    2028             double fakeValue=fakeSolution[iColumn];
     2101            CoinWorkDouble fakeValue=fakeSolution[iColumn];
    20292102            if (fakeValue<lowerValue+setPrimal) {
    20302103              fakeValue=lowerValue+setPrimal;
     
    20392112        } else {
    20402113          //just lower bound
    2041           double fakeValue=fakeSolution[iColumn];
     2114          CoinWorkDouble fakeValue=fakeSolution[iColumn];
    20422115          if (fakeValue<lowerValue+setPrimal) {
    20432116            fakeValue=lowerValue+setPrimal;
     
    20482121        if (upperBound(iColumn)) {
    20492122          //just upper bound
    2050           double fakeValue=fakeSolution[iColumn];
     2123          CoinWorkDouble fakeValue=fakeSolution[iColumn];
    20512124          if (fakeValue>upperValue-setPrimal) {
    20522125            fakeValue=upperValue-setPrimal;
     
    20932166    quadraticElement = quadratic->getElements();
    20942167  }
    2095   double quadraticNorm=0.0;
     2168  CoinWorkDouble quadraticNorm=0.0;
    20962169  for ( iColumn=0;iColumn<numberTotal;iColumn++) {
    20972170    if (!flagged(iColumn)) {
    2098       double primalValue=solution_[iColumn];
    2099       double lowerValue=lower_[iColumn];
    2100       double upperValue=upper_[iColumn];
     2171      CoinWorkDouble primalValue=solution_[iColumn];
     2172      CoinWorkDouble lowerValue=lower_[iColumn];
     2173      CoinWorkDouble upperValue=upper_[iColumn];
    21012174      // Do dj
    2102       double reducedCost=cost_[iColumn];
     2175      CoinWorkDouble reducedCost=cost_[iColumn];
    21032176      if (lowerBound(iColumn)) {
    21042177        reducedCost+=linearPerturbation_;
     
    21112184             j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
    21122185          int jColumn = columnQuadratic[j];
    2113           double valueJ = solution_[jColumn];
    2114           double elementValue = quadraticElement[j];
     2186          CoinWorkDouble valueJ = solution_[jColumn];
     2187          CoinWorkDouble elementValue = quadraticElement[j];
    21152188          reducedCost += valueJ*elementValue;
    21162189        }
    2117         quadraticNorm = CoinMax(quadraticNorm,fabs(reducedCost));
     2190        quadraticNorm = CoinMax(quadraticNorm,CoinAbs(reducedCost));
    21182191      }
    21192192      dj_[iColumn]=reducedCost;
     
    21322205  for ( iColumn=0;iColumn<numberTotal;iColumn++) {
    21332206    if (!flagged(iColumn)) {
    2134       double primalValue=solution_[iColumn];
    2135       double lowerValue=lower_[iColumn];
    2136       double upperValue=upper_[iColumn];
    2137       double reducedCost=dj_[iColumn];
    2138       double low=0.0;
    2139       double high=0.0;
     2207      CoinWorkDouble primalValue=solution_[iColumn];
     2208      CoinWorkDouble lowerValue=lower_[iColumn];
     2209      CoinWorkDouble upperValue=upper_[iColumn];
     2210      CoinWorkDouble reducedCost=dj_[iColumn];
     2211      CoinWorkDouble low=0.0;
     2212      CoinWorkDouble high=0.0;
    21402213      if (lowerBound(iColumn)) {
    21412214        if (upperBound(iColumn)) {
     
    21482221            high=initialValue;
    21492222          }
    2150           double s = low+extra;
    2151           double ratioZ;
     2223          CoinWorkDouble s = low+extra;
     2224          CoinWorkDouble ratioZ;
    21522225          if (s<zwLarge) {
    21532226            ratioZ=1.0;
    21542227          } else {
    2155             ratioZ=sqrt(zwLarge/s);
     2228            ratioZ=CoinSqrt(zwLarge/s);
    21562229          }
    2157           double t = high+extra;
    2158           double ratioT;
     2230          CoinWorkDouble t = high+extra;
     2231          CoinWorkDouble ratioT;
    21592232          if (t<zwLarge) {
    21602233            ratioT=1.0;
    21612234          } else {
    2162             ratioT=sqrt(zwLarge/t);
     2235            ratioT=CoinSqrt(zwLarge/t);
    21632236          }
    21642237          //modify s and t
     
    21792252            wVec_[iColumn]=CoinMax(-reducedCost , safeObjectiveValue*ratioT);
    21802253          }
    2181           double gammaTerm = gamma2;
     2254          CoinWorkDouble gammaTerm = gamma2;
    21822255          if (primalR_)
    21832256            gammaTerm += primalR_[iColumn];
     
    21882261          low=primalValue-lowerValue;
    21892262          high=0.0;
    2190           double s = low+extra;
    2191           double ratioZ;
     2263          CoinWorkDouble s = low+extra;
     2264          CoinWorkDouble ratioZ;
    21922265          if (s<zwLarge) {
    21932266            ratioZ=1.0;
    21942267          } else {
    2195             ratioZ=sqrt(zwLarge/s);
     2268            ratioZ=CoinSqrt(zwLarge/s);
    21962269          }
    21972270          //modify s
     
    22072280            wVec_[iColumn]=0.0;
    22082281          }
    2209           double gammaTerm = gamma2;
     2282          CoinWorkDouble gammaTerm = gamma2;
    22102283          if (primalR_)
    22112284            gammaTerm += primalR_[iColumn];
     
    22172290          low=0.0;
    22182291          high=upperValue-primalValue;
    2219           double t = high+extra;
    2220           double ratioT;
     2292          CoinWorkDouble t = high+extra;
     2293          CoinWorkDouble ratioT;
    22212294          if (t<zwLarge) {
    22222295            ratioT=1.0;
    22232296          } else {
    2224             ratioT=sqrt(zwLarge/t);
     2297            ratioT=CoinSqrt(zwLarge/t);
    22252298          }
    22262299          //modify t
     
    22362309            wVec_[iColumn]=CoinMax(-reducedCost , safeObjectiveValue*ratioT);
    22372310          }
    2238           double gammaTerm = gamma2;
     2311          CoinWorkDouble gammaTerm = gamma2;
    22392312          if (primalR_)
    22402313            gammaTerm += primalR_[iColumn];
     
    22492322  if (solution_[0]>0.0) {
    22502323    for (int i=0;i<numberTotal;i++)
    2251       printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n",i,fabs(solution_[i]),
    2252              diagonal_[i],fabs(dj_[i]),
     2324      printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n",i,CoinAbs(solution_[i]),
     2325             diagonal_[i],CoinAbs(dj_[i]),
    22532326             lowerSlack_[i],zVec_[i],
    22542327             upperSlack_[i],wVec_[i]);
    22552328  } else {
    22562329    for (int i=0;i<numberTotal;i++)
    2257       printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n",i,fabs(solution_[i]),
    2258              diagonal_[i],fabs(dj_[i]),
     2330      printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n",i,CoinAbs(solution_[i]),
     2331             diagonal_[i],CoinAbs(dj_[i]),
    22592332             upperSlack_[i],wVec_[i],
    22602333             lowerSlack_[i],zVec_[i] );
     
    22662339// complementarityGap.  Computes gap
    22672340//phase 0=as is , 1 = after predictor , 2 after corrector
    2268 double ClpPredictorCorrector::complementarityGap(int & numberComplementarityPairs,
     2341CoinWorkDouble ClpPredictorCorrector::complementarityGap(int & numberComplementarityPairs,
    22692342                                                 int & numberComplementarityItems,
    22702343                                                 const int phase)
    22712344{
    2272   double gap=0.0;
     2345  CoinWorkDouble gap=0.0;
    22732346  //seems to be same coding for phase = 1 or 2
    22742347  numberComplementarityPairs=0;
    22752348  numberComplementarityItems=0;
    22762349  int numberTotal = numberRows_+numberColumns_;
    2277   double toleranceGap=0.0;
    2278   double largestGap=0.0;
    2279   double smallestGap=COIN_DBL_MAX;
     2350  CoinWorkDouble toleranceGap=0.0;
     2351  CoinWorkDouble largestGap=0.0;
     2352  CoinWorkDouble smallestGap=COIN_DBL_MAX;
    22802353  //seems to be same coding for phase = 1 or 2
    22812354  int numberNegativeGaps=0;
    2282   double sumNegativeGap=0.0;
    2283   double largeGap=1.0e2*solutionNorm_;
     2355  CoinWorkDouble sumNegativeGap=0.0;
     2356  CoinWorkDouble largeGap=1.0e2*solutionNorm_;
    22842357  if (largeGap<1.0e10) {
    22852358    largeGap=1.0e10;
    22862359  }
    22872360  largeGap=1.0e30;
    2288   double dualTolerance =  dblParam_[ClpDualTolerance];
    2289   double primalTolerance =  dblParam_[ClpPrimalTolerance];
     2361  CoinWorkDouble dualTolerance =  dblParam_[ClpDualTolerance];
     2362  CoinWorkDouble primalTolerance =  dblParam_[ClpPrimalTolerance];
    22902363  dualTolerance=dualTolerance/scaleFactor_;
    22912364  for (int iColumn=0;iColumn<numberTotal;iColumn++) {
     
    22932366      numberComplementarityPairs++;
    22942367      //can collapse as if no lower bound both zVec and deltaZ 0.0
    2295       double newZ=0.0;
    2296       double newW=0.0;
     2368      CoinWorkDouble newZ=0.0;
     2369      CoinWorkDouble newW=0.0;
    22972370      if (lowerBound(iColumn)) {
    22982371        numberComplementarityItems++;
    2299         double dualValue;
    2300         double primalValue;
     2372        CoinWorkDouble dualValue;
     2373        CoinWorkDouble primalValue;
    23012374        if (!phase) {
    23022375          dualValue=zVec_[iColumn];
    23032376          primalValue=lowerSlack_[iColumn];
    23042377        } else {
    2305           double change;
     2378          CoinWorkDouble change;
    23062379          change =solution_[iColumn]+deltaX_[iColumn]-lowerSlack_[iColumn]-lower_[iColumn];
    23072380          dualValue=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];
     
    23132386          primalValue=largeGap;
    23142387        }
    2315         double gapProduct=dualValue*primalValue;
     2388        CoinWorkDouble gapProduct=dualValue*primalValue;
    23162389        if (gapProduct<0.0) {
    23172390          //cout<<"negative gap component "<<iColumn<<" "<<dualValue<<" "<<
     
    23342407      if (upperBound(iColumn)) {
    23352408        numberComplementarityItems++;
    2336         double dualValue;
    2337         double primalValue;
     2409        CoinWorkDouble dualValue;
     2410        CoinWorkDouble primalValue;
    23382411        if (!phase) {
    23392412          dualValue=wVec_[iColumn];
    23402413          primalValue=upperSlack_[iColumn];
    23412414        } else {
    2342           double change;
     2415          CoinWorkDouble change;
    23432416          change =upper_[iColumn]-solution_[iColumn]-deltaX_[iColumn]-upperSlack_[iColumn];
    23442417          dualValue=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn];
     
    23502423          primalValue=largeGap;
    23512424        }
    2352         double gapProduct=dualValue*primalValue;
     2425        CoinWorkDouble gapProduct=dualValue*primalValue;
    23532426        if (gapProduct<0.0) {
    23542427          //cout<<"negative gap component "<<iColumn<<" "<<dualValue<<" "<<
     
    23742447  if (!phase&&numberNegativeGaps) {
    23752448      handler_->message(CLP_BARRIER_NEGATIVE_GAPS,messages_)
    2376     <<numberNegativeGaps<<sumNegativeGap
     2449        <<numberNegativeGaps<<static_cast<double>(sumNegativeGap)
    23772450    <<CoinMessageEol;
    23782451  }
     
    23932466void ClpPredictorCorrector::setupForSolve(const int phase)
    23942467{
    2395   double extra =eExtra;
     2468  CoinWorkDouble extra =eExtra;
    23962469  int numberTotal = numberRows_ + numberColumns_;
    23972470  int iColumn;
     
    23992472  printf("phase %d in setupForSolve, mu %.18g\n",phase,mu_);
    24002473#endif
    2401   double gamma2 = gamma_*gamma_; // gamma*gamma will be added to diagonal
     2474  CoinWorkDouble gamma2 = gamma_*gamma_; // gamma*gamma will be added to diagonal
     2475  CoinWorkDouble * dualArray = reinterpret_cast<CoinWorkDouble *>(dual_);
    24022476  switch (phase) {
    24032477  case 0:
     
    24052479    if (delta_||dualR_) {
    24062480      // add in regularization
    2407       double delta2 = delta_*delta_;
     2481      CoinWorkDouble delta2 = delta_*delta_;
    24082482      for (int iRow=0;iRow<numberRows_;iRow++) {
    2409         rhsB_[iRow] -= delta2*dual_[iRow];
     2483        rhsB_[iRow] -= delta2*dualArray[iRow];
    24102484        if (dualR_)
    2411           rhsB_[iRow] -= dualR_[iRow]*dual_[iRow];
     2485          rhsB_[iRow] -= dualR_[iRow]*dualArray[iRow];
    24122486      }
    24132487    }
     
    24452519#if 0
    24462520    for (int i=0;i<3;i++) {
    2447       if (!fabs(rhsZ_[i]))
     2521      if (!CoinAbs(rhsZ_[i]))
    24482522        rhsZ_[i]=0.0;
    2449       if (!fabs(rhsW_[i]))
     2523      if (!CoinAbs(rhsW_[i]))
    24502524        rhsW_[i]=0.0;
    2451       if (!fabs(rhsU_[i]))
     2525      if (!CoinAbs(rhsU_[i]))
    24522526        rhsU_[i]=0.0;
    2453       if (!fabs(rhsL_[i]))
     2527      if (!CoinAbs(rhsL_[i]))
    24542528        rhsL_[i]=0.0;
    24552529    }
     
    24992573#if 0
    25002574    for (int i=0;i<numberTotal;i++) {
    2501       if (!fabs(rhsZ_[i]))
     2575      if (!CoinAbs(rhsZ_[i]))
    25022576        rhsZ_[i]=0.0;
    2503       if (!fabs(rhsW_[i]))
     2577      if (!CoinAbs(rhsW_[i]))
    25042578        rhsW_[i]=0.0;
    2505       if (!fabs(rhsU_[i]))
     2579      if (!CoinAbs(rhsU_[i]))
    25062580        rhsU_[i]=0.0;
    2507       if (!fabs(rhsL_[i]))
     2581      if (!CoinAbs(rhsL_[i]))
    25082582        rhsL_[i]=0.0;
    25092583    }
    25102584    if (solution_[0]>0.0) {
    25112585      for (int i=0;i<numberTotal;i++)
    2512         printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n",i,fabs(solution_[i]),
    2513                diagonal_[i],fabs(dj_[i]),
     2586        printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n",i,CoinAbs(solution_[i]),
     2587               diagonal_[i],CoinAbs(dj_[i]),
    25142588               lowerSlack_[i],zVec_[i],
    25152589               upperSlack_[i],wVec_[i]);
    25162590      for (int i=0;i<numberTotal;i++)
    2517         printf("%d %.18g %.18g %.18g %.18g %.18g\n",i,fabs(rhsC_[i]),
     2591        printf("%d %.18g %.18g %.18g %.18g %.18g\n",i,CoinAbs(rhsC_[i]),
    25182592               rhsZ_[i],rhsL_[i],
    25192593               rhsW_[i],rhsU_[i]);
    25202594    } else {
    25212595      for (int i=0;i<numberTotal;i++)
    2522         printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n",i,fabs(solution_[i]),
    2523                diagonal_[i],fabs(dj_[i]),
     2596        printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n",i,CoinAbs(solution_[i]),
     2597               diagonal_[i],CoinAbs(dj_[i]),
    25242598               upperSlack_[i],wVec_[i],
    25252599               lowerSlack_[i],zVec_[i] );
    25262600      for (int i=0;i<numberTotal;i++)
    2527         printf("%d %.18g %.18g %.18g %.18g %.18g\n",i,fabs(rhsC_[i]),
     2601        printf("%d %.18g %.18g %.18g %.18g %.18g\n",i,CoinAbs(rhsC_[i]),
    25282602               rhsW_[i],rhsU_[i],
    25292603               rhsZ_[i],rhsL_[i]);
     
    25492623  case 3:
    25502624    {
    2551       double minBeta = 0.1*mu_;
    2552       double maxBeta = 10.0*mu_;
    2553       double dualStep = CoinMin(1.0,actualDualStep_+0.1);
    2554       double primalStep = CoinMin(1.0,actualPrimalStep_+0.1);
     2625      CoinWorkDouble minBeta = 0.1*mu_;
     2626      CoinWorkDouble maxBeta = 10.0*mu_;
     2627      CoinWorkDouble dualStep = CoinMin(1.0,actualDualStep_+0.1);
     2628      CoinWorkDouble primalStep = CoinMin(1.0,actualPrimalStep_+0.1);
    25552629#ifdef SOME_DEBUG
    25562630      printf("good complementarity range %g to %g\n",minBeta,maxBeta);
     
    25612635        if (!flagged(iColumn)) {
    25622636          if (lowerBound(iColumn)) {
    2563             double change = -rhsL_[iColumn] + deltaX_[iColumn];
    2564             double dualValue=zVec_[iColumn]+dualStep*deltaZ_[iColumn];
    2565             double primalValue=lowerSlack_[iColumn]+primalStep*change;
    2566             double gapProduct=dualValue*primalValue;
     2637            CoinWorkDouble change = -rhsL_[iColumn] + deltaX_[iColumn];
     2638            CoinWorkDouble dualValue=zVec_[iColumn]+dualStep*deltaZ_[iColumn];
     2639            CoinWorkDouble primalValue=lowerSlack_[iColumn]+primalStep*change;
     2640            CoinWorkDouble gapProduct=dualValue*primalValue;
    25672641            if (gapProduct>0.0&&dualValue<0.0)
    25682642              gapProduct = - gapProduct;
     
    25732647                     iColumn,primalValue,dualValue,gapProduct);
    25742648#endif
    2575             double value= 0.0;
     2649            CoinWorkDouble value= 0.0;
    25762650            if (gapProduct<minBeta) {
    25772651              value= 2.0*(minBeta-gapProduct);
     
    25862660          } 
    25872661          if (upperBound(iColumn)) {
    2588             double change = rhsU_[iColumn]-deltaX_[iColumn];
    2589             double dualValue=wVec_[iColumn]+dualStep*deltaW_[iColumn];
    2590             double primalValue=upperSlack_[iColumn]+primalStep*change;
    2591             double gapProduct=dualValue*primalValue;
     2662            CoinWorkDouble change = rhsU_[iColumn]-deltaX_[iColumn];
     2663            CoinWorkDouble dualValue=wVec_[iColumn]+dualStep*deltaW_[iColumn];
     2664            CoinWorkDouble primalValue=upperSlack_[iColumn]+primalStep*change;
     2665            CoinWorkDouble gapProduct=dualValue*primalValue;
    25922666            if (gapProduct>0.0&&dualValue<0.0)
    25932667              gapProduct = - gapProduct;
     
    25982672                     iColumn,primalValue,dualValue,gapProduct);
    25992673#endif
    2600             double value= 0.0;
     2674            CoinWorkDouble value= 0.0;
    26012675            if (gapProduct<minBeta) {
    26022676              value= (minBeta-gapProduct);
     
    26152689  if (cholesky_->type()<20) {
    26162690    for (iColumn=0;iColumn<numberTotal;iColumn++) {
    2617       double value = rhsC_[iColumn];
    2618       double zValue = rhsZ_[iColumn];
    2619       double wValue = rhsW_[iColumn];
     2691      CoinWorkDouble value = rhsC_[iColumn];
     2692      CoinWorkDouble zValue = rhsZ_[iColumn];
     2693      CoinWorkDouble wValue = rhsW_[iColumn];
    26202694#if 0
    26212695#if 1
     
    26442718#else
    26452719      if (lowerBound(iColumn)) {
    2646         double gHat = zValue + zVec_[iColumn]*rhsL_[iColumn];
     2720        CoinWorkDouble gHat = zValue + zVec_[iColumn]*rhsL_[iColumn];
    26472721        value -= gHat/(lowerSlack_[iColumn]+extra);
    26482722      }
    26492723      if (upperBound(iColumn)) {
    2650         double hHat = wValue - wVec_[iColumn]*rhsU_[iColumn];
     2724        CoinWorkDouble hHat = wValue - wVec_[iColumn]*rhsU_[iColumn];
    26512725        value += hHat/(upperSlack_[iColumn]+extra);
    26522726      }
     
    26672741    // KKT
    26682742    for (iColumn=0;iColumn<numberTotal;iColumn++) {
    2669       double value = rhsC_[iColumn];
    2670       double zValue = rhsZ_[iColumn];
    2671       double wValue = rhsW_[iColumn];
     2743      CoinWorkDouble value = rhsC_[iColumn];
     2744      CoinWorkDouble zValue = rhsZ_[iColumn];
     2745      CoinWorkDouble wValue = rhsW_[iColumn];
    26722746      if (lowerBound(iColumn)) {
    2673         double gHat = zValue + zVec_[iColumn]*rhsL_[iColumn];
     2747        CoinWorkDouble gHat = zValue + zVec_[iColumn]*rhsL_[iColumn];
    26742748        value -= gHat/(lowerSlack_[iColumn]+extra);
    26752749      }
    26762750      if (upperBound(iColumn)) {
    2677         double hHat = wValue - wVec_[iColumn]*rhsU_[iColumn];
     2751        CoinWorkDouble hHat = wValue - wVec_[iColumn]*rhsU_[iColumn];
    26782752        value += hHat/(upperSlack_[iColumn]+extra);
    26792753      }
     
    26842758//method: sees if looks plausible change in complementarity
    26852759bool ClpPredictorCorrector::checkGoodMove(const bool doCorrector,
    2686                                           double & bestNextGap,
     2760                                          CoinWorkDouble & bestNextGap,
    26872761                                          bool allowIncreasingGap)
    26882762{
    2689   const double beta3 = 0.99997;
     2763  const CoinWorkDouble beta3 = 0.99997;
    26902764  bool goodMove=false;
    26912765  int nextNumber;
    26922766  int nextNumberItems;
    26932767  int numberTotal = numberRows_+numberColumns_;
    2694   double returnGap=bestNextGap;
    2695   double nextGap=complementarityGap(nextNumber,nextNumberItems,2);
     2768  CoinWorkDouble returnGap=bestNextGap;
     2769  CoinWorkDouble nextGap=complementarityGap(nextNumber,nextNumberItems,2);
    26962770#ifndef NO_RTTI
    26972771  ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
     
    27112785    returnGap=nextGap;
    27122786  }
    2713   double step;
     2787  CoinWorkDouble step;
    27142788  if (actualDualStep_>actualPrimalStep_) {
    27152789    step=actualDualStep_;
     
    27172791    step=actualPrimalStep_;
    27182792  }
    2719   double testValue=1.0-step*(1.0-beta3);
     2793  CoinWorkDouble testValue=1.0-step*(1.0-beta3);
    27202794  //testValue=0.0;
    27212795  testValue*=complementarityGap_;
     
    27292803      //step=actualPrimalStep_;
    27302804    //}
    2731     double gap = bestNextGap;
     2805    CoinWorkDouble gap = bestNextGap;
    27322806    goodMove=checkGoodMove2(step,gap,allowIncreasingGap);
    27332807    if (goodMove)
     
    27602834    while (!goodMove) {
    27612835      pass++;
    2762       double gap = bestNextGap;
     2836      CoinWorkDouble gap = bestNextGap;
    27632837      goodMove=checkGoodMove2(step,gap,allowIncreasingGap);
    27642838      if (goodMove||pass>3) {
     
    27932867  if (goodMove) {
    27942868    //compute delta in objectives
    2795     double deltaObjectivePrimal=0.0;
    2796     double deltaObjectiveDual=
     2869    CoinWorkDouble deltaObjectivePrimal=0.0;
     2870    CoinWorkDouble deltaObjectiveDual=
    27972871      innerProduct(deltaY_,numberRows_,
    27982872                   rhsFixRegion_);
    2799     double error=0.0;
    2800     double * workArray = workArray_;
     2873    CoinWorkDouble error=0.0;
     2874    CoinWorkDouble * workArray = workArray_;
    28012875    CoinZeroN(workArray,numberColumns_);
    28022876    CoinMemcpyN(deltaY_,numberRows_,workArray+numberColumns_);
    28032877    matrix_->transposeTimes(-1.0,deltaY_,workArray);
    2804     //double sumPerturbCost=0.0;
     2878    //CoinWorkDouble sumPerturbCost=0.0;
    28052879    for (int iColumn=0;iColumn<numberTotal;iColumn++) {
    28062880      if (!flagged(iColumn)) {
     
    28132887          deltaObjectiveDual-=deltaW_[iColumn]*upper_[iColumn];
    28142888        }
    2815         double change = fabs(workArray_[iColumn]-deltaZ_[iColumn]+deltaW_[iColumn]);
     2889        CoinWorkDouble change = CoinAbs(workArray_[iColumn]-deltaZ_[iColumn]+deltaW_[iColumn]);
    28162890        error = CoinMax(change,error);
    28172891      }
     
    28192893    }
    28202894    //deltaObjectivePrimal+=sumPerturbCost*linearPerturbation_;
    2821     double testValue;
     2895    CoinWorkDouble testValue;
    28222896    if (error>0.0) {
    28232897      testValue=1.0e1*CoinMax(maximumDualError_,1.0e-12)/error;
     
    28282902    if (testValue<actualDualStep_&&!quadraticObj) {
    28292903      handler_->message(CLP_BARRIER_REDUCING,messages_)
    2830       <<"dual"<<actualDualStep_
    2831       << testValue
     2904        <<"dual"<<static_cast<double>(actualDualStep_)
     2905        << static_cast<double>(testValue)
    28322906      <<CoinMessageEol;
    28332907      actualDualStep_=testValue;
     
    28382912    //check change in AX not too much
    28392913    //??? could be dropped row going infeasible
    2840     double ratio = 1.0e1*CoinMax(maximumRHSError_,1.0e-12)/maximumRHSChange_;
     2914    CoinWorkDouble ratio = 1.0e1*CoinMax(maximumRHSError_,1.0e-12)/maximumRHSChange_;
    28412915    if (ratio<actualPrimalStep_) {
    28422916      handler_->message(CLP_BARRIER_REDUCING,messages_)
    2843       <<"primal"<<actualPrimalStep_
    2844       <<ratio
     2917        <<"primal"<<static_cast<double>(actualPrimalStep_)
     2918        <<static_cast<double>(ratio)
    28452919      <<CoinMessageEol;
    28462920      if (ratio>1.0e-6) {
     
    28572931}
    28582932//:  checks for one step size
    2859 bool ClpPredictorCorrector::checkGoodMove2(double move,
    2860                                            double & bestNextGap,
     2933bool ClpPredictorCorrector::checkGoodMove2(CoinWorkDouble move,
     2934                                           CoinWorkDouble & bestNextGap,
    28612935                                           bool allowIncreasingGap)
    28622936{
    2863   double complementarityMultiplier =1.0/numberComplementarityPairs_;
    2864   const double gamma = 1.0e-8;
    2865   const double gammap = 1.0e-8;
    2866   double gammad = 1.0e-8;
     2937  CoinWorkDouble complementarityMultiplier =1.0/numberComplementarityPairs_;
     2938#if CLP_LONG_CHOLESKY<2
     2939  const CoinWorkDouble gamma = 1.0e-8;
     2940  const CoinWorkDouble gammap = 1.0e-8;
     2941  CoinWorkDouble gammad = 1.0e-8;
     2942#else
     2943  const CoinWorkDouble gamma = 1.0e-8;
     2944  const CoinWorkDouble gammap = 1.0e-8;
     2945  CoinWorkDouble gammad = 1.0e-8;
     2946#endif
    28672947  int nextNumber;
    28682948  int nextNumberItems;
    2869   double nextGap=complementarityGap(nextNumber,nextNumberItems,2);
     2949  CoinWorkDouble nextGap=complementarityGap(nextNumber,nextNumberItems,2);
    28702950  if (nextGap>bestNextGap&&!allowIncreasingGap)
    28712951    return false;
    2872   double lowerBoundGap = gamma*nextGap*complementarityMultiplier;
     2952  CoinWorkDouble lowerBoundGap = gamma*nextGap*complementarityMultiplier;
    28732953  bool goodMove=true;
    28742954  int iColumn;
     
    28762956    if (!flagged(iColumn)) {
    28772957      if (lowerBound(iColumn)) {
    2878         double part1=lowerSlack_[iColumn]+actualPrimalStep_*deltaSL_[iColumn];
    2879         double part2=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];
     2958        CoinWorkDouble part1=lowerSlack_[iColumn]+actualPrimalStep_*deltaSL_[iColumn];
     2959        CoinWorkDouble part2=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];
    28802960        if (part1*part2<lowerBoundGap) {
    28812961          goodMove=false;
     
    28842964      }
    28852965      if (upperBound(iColumn)) {
    2886         double part1=upperSlack_[iColumn]+actualPrimalStep_*deltaSU_[iColumn];
    2887         double part2=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn];
     2966        CoinWorkDouble part1=upperSlack_[iColumn]+actualPrimalStep_*deltaSU_[iColumn];
     2967        CoinWorkDouble part2=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn];
    28882968        if (part1*part2<lowerBoundGap) {
    28892969          goodMove=false;
     
    28932973    }
    28942974  }
    2895    double * nextDj=NULL;
    2896    double maximumDualError = maximumDualError_;
     2975   CoinWorkDouble * nextDj=NULL;
     2976   CoinWorkDouble maximumDualError = maximumDualError_;
    28972977#ifndef NO_RTTI
    28982978  ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
     
    29022982    quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
    29032983#endif
     2984  CoinWorkDouble * dualArray = reinterpret_cast<CoinWorkDouble *>(dual_);
    29042985  if (quadraticObj) {
    29052986    // change gammad
    29062987    gammad=1.0e-4;
    2907     double gamma2 = gamma_*gamma_;
    2908     nextDj = new double [numberColumns_];
    2909     double * nextSolution = new double [numberColumns_];
     2988    CoinWorkDouble gamma2 = gamma_*gamma_;
     2989    nextDj = new CoinWorkDouble [numberColumns_];
     2990    CoinWorkDouble * nextSolution = new CoinWorkDouble [numberColumns_];
    29102991    // put next primal into nextSolution
    29112992    for ( iColumn=0;iColumn<numberColumns_;iColumn++) {
     
    29193000    // do reduced costs
    29203001    CoinMemcpyN(cost_,numberColumns_,nextDj);
    2921     matrix_->transposeTimes(-1.0,dual_,nextDj);
     3002    matrix_->transposeTimes(-1.0,dualArray,nextDj);
    29223003    matrix_->transposeTimes(-actualDualStep_,deltaY_,nextDj);
    29233004    quadraticDjs(nextDj,nextSolution,1.0);
     
    29273008    for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
    29283009      if (!fixedOrFree(iColumn)) {
    2929         double newZ=0.0;
    2930         double newW=0.0;
     3010        CoinWorkDouble newZ=0.0;
     3011        CoinWorkDouble newW=0.0;
    29313012        if (lowerBound(iColumn)) {
    29323013          newZ=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];
     
    29363017        }
    29373018        if (columnQuadraticLength[iColumn]) {
    2938           double gammaTerm = gamma2;
     3019          CoinWorkDouble gammaTerm = gamma2;
    29393020          if (primalR_)
    29403021            gammaTerm += primalR_[iColumn];
    2941           //double dualInfeasibility=
     3022          //CoinWorkDouble dualInfeasibility=
    29423023          //dj_[iColumn]-zVec_[iColumn]+wVec_[iColumn]
    29433024          //+gammaTerm*solution_[iColumn];
    2944           double newInfeasibility=
     3025          CoinWorkDouble newInfeasibility=
    29453026            nextDj[iColumn]-newZ+newW
    29463027            +gammaTerm*(solution_[iColumn]+actualPrimalStep_*deltaX_[iColumn]);
    29473028          maximumDualError = CoinMax(maximumDualError,newInfeasibility);
    2948           //if (fabs(newInfeasibility)>CoinMax(2000.0*maximumDualError_,1.0e-2)) {
     3029          //if (CoinAbs(newInfeasibility)>CoinMax(2000.0*maximumDualError_,1.0e-2)) {
    29493030          //if (dualInfeasibility*newInfeasibility<0.0) {
    29503031          //  printf("%d current %g next %g\n",iColumn,dualInfeasibility,
     
    29623043    solutionNorm_=rhsNorm_;
    29633044  }
    2964   double errorCheck=maximumRHSError_/solutionNorm_;
     3045  CoinWorkDouble errorCheck=maximumRHSError_/solutionNorm_;
    29653046  if (errorCheck<maximumBoundInfeasibility_) {
    29663047    errorCheck=maximumBoundInfeasibility_;
     
    29873068// updateSolution.  Updates solution at end of iteration
    29883069//returns number fixed
    2989 int ClpPredictorCorrector::updateSolution(double nextGap)
     3070int ClpPredictorCorrector::updateSolution(CoinWorkDouble nextGap)
    29903071{
     3072  CoinWorkDouble * dualArray = reinterpret_cast<CoinWorkDouble *>(dual_);
    29913073  int numberTotal = numberRows_+numberColumns_;
    29923074  //update pi
    2993   multiplyAdd(deltaY_,numberRows_,actualDualStep_,dual_,1.0);
     3075  multiplyAdd(deltaY_,numberRows_,actualDualStep_,dualArray,1.0);
    29943076  CoinZeroN(errorRegion_,numberRows_);
    29953077  CoinZeroN(rhsFixRegion_,numberRows_);
    2996   double maximumRhsInfeasibility=0.0;
    2997   double maximumBoundInfeasibility=0.0;
    2998   double maximumDualError=1.0e-12;
    2999   double primalObjectiveValue=0.0;
    3000   double dualObjectiveValue=0.0;
    3001   double solutionNorm=1.0e-12;
     3078  CoinWorkDouble maximumRhsInfeasibility=0.0;
     3079  CoinWorkDouble maximumBoundInfeasibility=0.0;
     3080  CoinWorkDouble maximumDualError=1.0e-12;
     3081  CoinWorkDouble primalObjectiveValue=0.0;
     3082  CoinWorkDouble dualObjectiveValue=0.0;
     3083  CoinWorkDouble solutionNorm=1.0e-12;
    30023084  int numberKilled=0;
    3003   double freeMultiplier=1.0e6;
    3004   double trueNorm =diagonalNorm_/diagonalScaleFactor_;
     3085  CoinWorkDouble freeMultiplier=1.0e6;
     3086  CoinWorkDouble trueNorm =diagonalNorm_/diagonalScaleFactor_;
    30053087  if (freeMultiplier<trueNorm) {
    30063088    freeMultiplier=trueNorm;
     
    30103092  }
    30113093  freeMultiplier=0.5/freeMultiplier;
    3012   double condition = fabs(cholesky_->choleskyCondition());
     3094  CoinWorkDouble condition = CoinAbs(cholesky_->choleskyCondition());
    30133095  bool caution;
    30143096  if ((condition<1.0e10&&trueNorm<1.0e12)||numberIterations_<20) {
     
    30173099    caution=true;
    30183100  }
    3019   double extra=eExtra;
    3020   const double largeFactor=1.0e2;
    3021   double largeGap=largeFactor*solutionNorm_;
     3101  CoinWorkDouble extra=eExtra;
     3102  const CoinWorkDouble largeFactor=1.0e2;
     3103  CoinWorkDouble largeGap=largeFactor*solutionNorm_;
    30223104  if (largeGap<largeFactor) {
    30233105    largeGap=largeFactor;
    30243106  }
    3025   double dualFake=0.0;
    3026   double dualTolerance =  dblParam_[ClpDualTolerance];
     3107  CoinWorkDouble dualFake=0.0;
     3108  CoinWorkDouble dualTolerance =  dblParam_[ClpDualTolerance];
    30273109  dualTolerance=dualTolerance/scaleFactor_;
    30283110  if (dualTolerance<1.0e-12) {
    30293111    dualTolerance=1.0e-12;
    30303112  }
    3031   double offsetObjective=0.0;
    3032   const double killTolerance=primalTolerance();
    3033   double qDiagonal;
     3113  CoinWorkDouble offsetObjective=0.0;
     3114  const CoinWorkDouble killTolerance=primalTolerance();
     3115  CoinWorkDouble qDiagonal;
     3116#if CLP_LONG_CHOLESKY<2
    30343117  if (mu_<1.0) {
    30353118